Team:Paris/Modeling/

From 2008.igem.org

Estimation of the parameters

function Popt=findparam

function Popt=findparam

 

 

clear, close all

 

%Parameters

Tmax=100; %Time of the simulation

Npoints=40; %Number of points available. We assume for the tests that they are regularly distant.

x=1:Tmax/Npoints:Tmax;%Definition of the inducer vector.

 

%The following function provides a Hill function, given a X vector and

%some parameter that characterizes it.

function y= Hill(x,K,Vmax,n)

l=length(x);

for i=1:l

if K<0 s(i)=K; else s(i)=x(i); end

y(i)=(Vmax*(s(i)^n))/((x(i)^n)+(K^n));

end

end

 

%Then, we want to introduce some noise in the data obtained, in order to

%simulate biological results.

function y = NoisyHill(x,K,Vmax,n)

y=Hill(x,K,Vmax,n);

for i=1:length(x)

y(i)=random('normal',y(i),Vmax/20);

end

end

 

%Finally, we create a 'noisy Hill' vector that will be considered as

%biological data.

function set = genset

K=30;

Vmax=2;

n=3;

set=NoisyHill([1:Tmax/Npoints:Tmax],K,Vmax,n);

end

set=genset();

 

%The heart of the problem is to minimize the distance between a

%deterministic Hill function and the noisy one obtained. 'Toopt' is

%just the criteria to minimize.

function y = toopt(x,P)

y=Hill(x,10*P(1),P(2),P(3))-set;

end

%These lines use a Matlab function to carry out the optimization process

options=optimset('LevenbergMarquardt','on','TolX',1e-10,'MaxFunEvals',1e7,'TolFun',1e-5);

Popt=lsqnonlin(@(P)toopt(x,P),[0.1,1,1],[0,0,0],[5,10,10],options);

Popt=Popt(:);

 

%In order to see the results, we plot both the noisy curve and the

%estimated one.

plot(set)

hold on

plot(Hill(x,10*Popt(1),Popt(2),Popt(3) ),'r');

legend('Experimental Hill function','Estimation of the experimental hill function');

gtext('NPoints=40, Standard Deviation=1/10');

xlabel('Concentration of X');

ylabel('Rate of Y production');

end