IGEM:IMPERIAL/2006/Least Square Curve Fit

The Method

This is a brief tutorial on how to do the least square curve fit on Matlab to extract values.

Adopted from previous tutorials by Marcio von Muhlen, Maxim Shusteff, Nate Tedford, and other BE Students at MIT.

In order to generate a fit, you need to tell Matlab the form of the equation you believe the data should follow. Matlab has no way of knowing or telling you if your function is incorrect, so it's up to you to account for this. We first create a Matlab function file. For example, an equation of the form:

<amsmath>

f(x) = \frac{C_1 \times x}{1 + C_2 \times x}

</amsmath>

The function file might look like:

 function yhat = myfun(beta, x)
 c1 = beta(1);
 c2 = beta(2);
 yhat = (c1*x ./ (1 + c2.*x);

Since this is a function, it must be named myfun.m. The coefficients to be fitted are here contained in the input variable beta. Make sure to have this function in your working directory.

Create a matrix fitData which contains the information you want. You could use:

Some sample data can be found here.

 >>fitData = dlmread('fitdata.txt', '\t', 1, 0);

Assign the X column to the variable X and the f(X) column to the variable Y, then plot the data.

 >>X = fitData(:,1);
 >>Y = fitData(:,2);
 >>semilogx(X, Y, 'bd')

We need to provide Matlab with initial guesses for our coefficient values. At this point it's a good idea to plot our function to see if the data makes sense with regards to it. We can start by guessing 1 for both coefficients; in a real experiment, your prior knowledge of what you are measuring would determine these guesses.

 >>hold on
 >>semilogx(X, (X ./(1 + 1.*X))

Now let's have Matlab figure out the best fit (least squares sense) of the data to this function, then plot the data together.

 >>Cfs = lsqcurvefit(@myfun, [1 1], X, Y)
 >>hold off
 >>semilogx(X, Y, 'bd', X, (Cfs(1)*X ./ (1 + Cfs(2)*X)), 'k-')

Here we plotted two sets of dta in the same command. See the help file for semilogx to understand how this works. Notice that although our guess wasn't very good, Matlab had no trouble finding good coefficients. This won't always be the case with more complicated data functions.

Additional resources: MIT's Matlab Answer Page

Another Matlab tutorial, with additional material on looping and solving ODE's in Matlab.

J37016 Curve Fitting

The function that I have decided to fit is:

[math]\displaystyle{ f(x) = \frac{C_1 \times x}{C_2 + 0.93X} }[/math]

This is adapted from the equation

[math]\displaystyle{ \frac{d[GFP]}{dt} = \frac{V_{max}[AHL][LuxR]}{[AHL][LuxR] + \frac{K_D}{K_{D \alpha}}} - \delta_{GFP}[GFP] \lt br\gt }[/math]

We can now solve this differential equation for [GFP] to obtain a relation which we can use to fit the curve we obtain from experiment, since we obtained a graph of fluorescence (directly related to the concentration of GFP) vs. concentration.

We first let:

[math]\displaystyle{ \gamma = \frac{V_{max}[AHL][LuxR]}{[AHL][LuxR] + \frac{K_D}{K_{D \alpha}}} }[/math]

Now we solve the differential equation

[math]\displaystyle{ \frac{d[GFP]}{dt} + \delta_{GFP}[GFP] = \gamma }[/math]

From the initial conditions t = 0, [GFP] = 0, we can come up with the solution:

[math]\displaystyle{ [GFP] = \frac{\gamma}{\delta_{GFP}} \times (1 - exp(-\delta_{GFP}t)) }[/math]

I'm not sure this is correct, since we are plotting fluorescence vs. [AHL], so please update it if I'm wrong and I can see what I need to do...

Comments

VR:Well, the model you refer to has 3 variables ([AHL],[LuxR],[GFP]) and 2 parameters (Vmax, Kd). Your fitting function needs to have the same degrees of freedom. Your fitting function should be like:

[math]\displaystyle{ GFP = \frac{C_1 \times X \times Y}{C_2 + XY} }[/math]

Out of the fitting operation you want to get C_1 and C_2.Then you have to figure out how to constrain your function with the experimental data we have. Remember that we operated at steady state and with the assumption that the degradation of GFP equals degradation of LuxR (mostly due to dilution from growth).

Johnsy: I will try again with the model that you suggest. I thought that we were assuming that LuxR was a constant? How can we extract values from the graph of LuxR, since we only have a graph of GFP vs. [AHL]? Also, what I have done to constrain the parameters was to "ignore" everything before expression occurs. So in most of our cases, this would correspond to 10 nM (I think) before GFP production began to increase significantly.