Team:Paris/Modeling/Others prog

complexes

function complexes = complexes(x,y,K,n)
% amount of complexes n*x + y = x_-_y in function of x, y
 
% x, y = initial amounts of binding molecules
% K = dissociation constant
% n = stoechiometric coefficient of x
 
syms z; % to treat z as a symbolic variable
Eqn = ( ( x - n*z )^n )*( y - z ) - K*z; % to define the equation in whose
                                          % the amount of complexes is a root  
 
S = eval(solve( Eqn, 'z' )); % the vector of the roots
 
Index = find ((0 < S) & (n*S < x) & (S < y) & (imag(S) == 0)); % determine  
                                                                % which
                                                                % root is  
                                                                % acceptable
 
complexes = S(Index(1)); % the result
 
end

hill

function compl = hill(x,K,n)
% complexation nx + Y = x_-_Y when x is in large excess :
% gives the ratio (x_-_Y)/(Ytotal) in function of x
 
% x = molecule that binds to Y
% K^n = dissociation constant
% n = stoechiometric coefficient of x ; also known as 'cooperativity'
 
compl = x^n/(K^n + x^n);
 
end

Global

function labo=Core_Sytem_Simulation (t, init, prom1, prom2, retro)
% Simulation of our Core_System
 
% prom1 = 'pTet' or 'pFlhDC'
% prom2 = 'pFlgA' or 'pFlgB'
% retro = 'EnvZ' or 'OmpR' or quantity of aTc
 
% t = discret time scale
% init = vector of the initial state
%   [FlhDC;
%    FliA;
%    FP1;
%    FP2;
%    FP3;
%    TetR / OmpR / EnvZ] regards to the studied version
 
global gamma beta16 K13 n13 K12 n12 K15 n15 beta17 K6 n6 ...
     beta22 beta18 K1 n1 beta23 K7 n7 beta24 K2 n2 ...
     beta25 K8 n8 gamma37 beta26 K3 n3 beta27 K9 n9 ...
     gamma38 beta28 K4 n4 beta29 K10 n10 beta30 K5 n5 ...
     beta31 K11 n11 gamma39 EnvZ_b OmpR_b K14 n14;  
% A program must load all parameters !
 
if prom2 == 'pFlgA'
     F = @f6
elseif prom2 == 'pFlgA'
         F = @f7
     end
else error('Wrong promoter for FP2')
end
 
if prom1 == 'pTet'
     try aTc = eval(retro)
     catch error('Wrong aTc Value')
     end
       
     function ydot = deriv(t,y) % dy/dt
                                % pTet : y(6) = TetR ; retro = aTc
         ydot = [f1(y(6), aTc);
                 f4(y(1), y(2));
                 f5(y(1), y(2)) - gamma36;
                 F(y(1), y(2)) - gamma37;
                 f8(y(1), y(2)) - gamma38;
                 f8(y(1), y(2));] - gamma*ones(6,1)
     end
elseif prom1 == 'pFlhDC'
     if retro == 'OmpR'
         function ydot = deriv(t,y)  % dy/dt
                                     % pFlhDC : y(6) = OmpR
             ydot = [f3(y(2), y(6));
                     f4(y(1), y(2));
                     f5(y(1), y(2)) - gamma36;
                     F(y(1), y(2)) - gamma37;
                     f8(y(1), y(2)) - gamma38;
                     f8(y(1), y(2));] - gamma*ones(6,1)
         end
     elseif retro == 'EnvZ'
         function ydot = deriv(t,y)  % dy/dt
                                     % pFlhDC : y(6) = EnvZ
             ydot = [f3bis(y(6), y(2));
                     f4(y(1), y(2));
                     f5(y(1), y(2)) - gamma36;
                     F(y(1), y(2)) - gamma37;
                     f8(y(1), y(2)) - gamma38;
                     f8(y(1), y(2));] - gamma*ones(6,1)
         end
     else
         error('gene for negative feed-back')
     end
else
     error('Wrong promoter for flhDC')
     end
end
 
[t,labo]=ode45(deriv,t,init);
 
% The 'full modularity' of this Approach allow other variations,  
% like  changing the place of the negative feed-back
% (changing f8 by f7, f6, f5)...