IGEM:Melbourne/2008/BioClock/David Version/clock code

function xdot = Stephen_clock(time,x)

format long e; global pulse g_p;

% synopsis: % xdot = (time, x) % x0 = % % can be used with odeN functions as follows: % % x0 = ; % [t,x] = ode23s(@, [0 100], ); % plot (t,x); % % where 100 is the end time % % When is used without any arguments it returns a vector of % the initial concentrations of the 14 floating species. % Otherwise should be called with two arguments: % time and x. time is the current time. x is the vector of the % concentrations of the 14 floating species. % When these parameters are supplied returns a vector of % the rates of change of the concentrations of the 14 floating species. % % the following table shows the mapping between the vector % index of a floating species and the species name. % % NOTE for compartmental models % matlab translator generates code that when simulated in matlab, % produces results which have the units of species amounts. Users % should divide the results for each species with the volume of the % compartment it resides in, in order to obtain concentrations. % % Indx Name % x(1) EnvZ_dephos % x(2) OmpR_dephos % x(3) X % x(4) GFP % x(5) RiboKey1 % x(6) Y % x(7) RFP % x(8) Off_X % x(9) RiboKey2 % x(10) GFP_sink % x(11) RFP_sink % x(12) X_sink % x(13) Y_sink % x(14) EnvZ_dephos_sink % x(15) OmpR_dephos_sink % x(16) Intermediate1 % x(17) RiboKey1_sink % x(18) Intermediate2 % x(19) Intermediate3 % x(20) X_sink_617291

xdot = zeros(20, 1); g_p = ones(30,1);

% List of Compartments vol__compartment = 1;

% Global Parameters g_p(1) = 0.1; % j1_k g_p(2) = 0.1; % j10_k g_p(3) = 0.1; % j14_k g_p(4) = 0.1; % j16_k g_p(5) = 5; % j2_vmax g_p(6) = 0.2; % j2_km1 g_p(7) = 0.3; % j2_ki g_p(8) = 0.5; % j5_vmax g_p(9) = 0.3; % j5_km1 g_p(10) = 0.1; % j5_ki g_p(11) = 0.5; % j12_vmax g_p(12) = 0.3; % j12_km1 g_p(13) = 1; % j12_ki g_p(14) = 0.5; % j19_k g_p(15) = 0.1; % j20_k g_p(16) = 2; % j21_k g_p(17) = 0.5; % j8_vmax g_p(18) = 0.3; % j8_km1 g_p(19) = 0.1; % j8_ki g_p(20) = 0.5; % j9_vmax g_p(21) = 0.3; % j9_km1 g_p(22) = 0.1; % j9_ki g_p(23) = 0.1; % j17_k g_p(24) = 0.1; % j22_k g_p(25) = 1; % j18_k g_p(26) = 0.1; % j23_vmax g_p(27) = 1e-005; % j23_km1 g_p(28) = 0.1; % j24_k g_p(29) = 0.1; % j25_k g_p(30) = 5000; % j26_k

% Boundary Conditions g_p31 = 1; % envz_phos [Concentration] %g_p32 = 0; % red_light [Concentration] g_p33 = 1; % ompr_phos [Concentration] g_p34 = 0; % off_y [Concentration]

 if pulse == 1
      g_p32 = getSG0(time); %red_light similar pulse to Sg0 in the other model
      dlmwrite('time_v2.dat',time,'-append');
      dlmwrite('sg0_v2.dat',g_p32,'-append');
  else
      g_p32 = 0;		% red_light [Concentration]
      
 end

% Other Parameters

  if (nargin == 0)
   % set initial conditions
  xdot(1) = 0*vol__compartment;		% EnvZ_dephos [Concentration]
  xdot(2) = 0*vol__compartment;		% OmpR_dephos [Concentration]
  xdot(3) = 0*vol__compartment;		% X [Concentration]
  xdot(4) = 0*vol__compartment;		% GFP [Concentration]
  xdot(5) = 0*vol__compartment;		% RiboKey1 [Concentration]
  xdot(6) = 0*vol__compartment;		% Y [Concentration]
  xdot(7) = 0*vol__compartment;		% RFP [Concentration]
  xdot(8) = 0*vol__compartment;		% Off_X [Concentration]
  xdot(9) = 0*vol__compartment;		% RiboKey2 [Concentration]
  xdot(10) = 0*vol__compartment;		% GFP_sink [Concentration]
  xdot(11) = 0*vol__compartment;		% RFP_sink [Concentration]
  xdot(12) = 0*vol__compartment;		% X_sink [Concentration]
  xdot(13) = 0*vol__compartment;		% Y_sink [Concentration]
  xdot(14) = 0*vol__compartment;		% EnvZ_dephos_sink [Concentration]
  xdot(15) = 0*vol__compartment;		% OmpR_dephos_sink [Concentration]
  xdot(16) = 0*vol__compartment;		% Intermediate1 [Concentration]
  xdot(17) = 0*vol__compartment;		% RiboKey1_sink [Concentration]
  xdot(18) = 0*vol__compartment;		% Intermediate2 [Concentration]
  xdot(19) = 0*vol__compartment;		% Intermediate3 [Concentration]
  xdot(20) = 0*vol__compartment;		% X_sink_617291 [Concentration]

else

   % calculate rates of change
  R0 = g_p(1)*(x(1)/vol__compartment)*(g_p33/vol__compartment);
  R1 = g_p(5)*(x(2)/vol__compartment)/(g_p(6)*(1+(g_p33/vol__compartment)/g_p(7))+(x(2)/vol__compartment));
  R2 = g_p(8)*(x(3)/vol__compartment)/(g_p(9)*(1+(x(8)/vol__compartment)/g_p(10))+(x(3)/vol__compartment));
  R3 = g_p(17)*(x(6)/vol__compartment)/(g_p(18)*(1+(g_p34/vol__compartment)/g_p(19))+(x(6)/vol__compartment));
  R4 = g_p(20)*(x(6)/vol__compartment)/(g_p(21)*(1+(g_p34/vol__compartment)/g_p(22))+(x(6)/vol__compartment));
  R5 = g_p(2)*(x(6)/vol__compartment);
  R6 = g_p(11)*(x(3)/vol__compartment)/(g_p(12)*(1+(x(8)/vol__compartment)/g_p(13))+(x(3)/vol__compartment));
  R7 = g_p(3)*(g_p31/vol__compartment)*(g_p32/vol__compartment);
  R8 = g_p(4)*(x(8)/vol__compartment);
  R9 = g_p(25)*(x(4)/vol__compartment);
  R10 = g_p(14)*(x(7)/vol__compartment);
  R11 = g_p(15)*(x(3)/vol__compartment);
  R12 = g_p(16)*(x(6)/vol__compartment);
  R13 = g_p(23)*(x(1)/vol__compartment);
  R14 = g_p(24)*(x(2)/vol__compartment);
  R15 = g_p(26)*(x(3)/vol__compartment)/(g_p(27)+(x(3)/vol__compartment));
  R16 = g_p(28)*(x(5)/vol__compartment);
  R17 = g_p(29)*(x(16)/vol__compartment)*(x(2)/vol__compartment);
  R18 = g_p(30)*(x(16)/vol__compartment)*(x(3)/vol__compartment);

  xdot = [
     - R0 + R7 - R13
     + R0 - R1 - R14 - R17
     + R1 + R2 - R2 - R6 - R11 - R15 - R18
     + R6 - R9
     + 50*R15 - R16
     + R3 - R3 - R4 - R5 - R12 + 20*R17
     + R4 - R10
     + R8 - R8
     + R5
     + R9
     + R10
     + R11
     + R12
     + R13
     + R14
     + R16 - R17 - R18
       0
       0
       0
     + R18
  ];
  end;