Model for DTC matlab course

dydt.m

function res = dydt(t,y,S, w)
  indices2={[1,5], [11], [11], [2,5], [12], [12], [5,13], [12,6], [3,9], [14], [14], [4,9], [15], [15], [7,7], [12], [2,6], [13], [10,6], [16], [10,13], [17], [17], [8], [7], [6], [4], [3]};
  % w=ones(28,1);
  % w=randn(28,1);
  v = w;
  for j=1:28
    for i=indices2{j}
      v(j) = v(j) * y(i);
    end
  end
  res = S*v;
end

model.m

S=zeros(17,28);
indices1={[11], [1,5], [1,5,3], [12], [2,5], [2,5,4], [12,6], [5,13], [14], [3,9], [7,3,9], [15], [4,9], [8,4,9], [6], [7,7], [13], [2,6], [16], [10,6], [17], [10,13], [2,16], [], [], [], [], []};
indices2={[1,5], [11], [11], [2,5], [12], [12], [5,13], [12,6], [3,9], [14], [14], [4,9], [15], [15], [7,7], [12], [2,6], [13], [10,6], [16], [10,13], [17], [17], [8], [7], [6], [4], [3]};
for j=1:28
  for i = indices1{j}
    S(i,j) = S(i,j) + 1;
  end
end
for j=1:28
  for i = indices2{j}
    S(i,j) = S(i,j) + -1;
  end
end
% y0 = ones(17,1);
% y0 = randi([1,10], 17,1);
% y0(7) = 0; y0(3) = 0; y0(4) = 0; y0(6) = 0; y0(7) = 0; y0(8) = 0; y0(13) = 0; y0(14) = 0; y0(15) = 0; y0(16) = 0; y0(17) = 0;
y0 = zeros(17,1);
y0(1) = 0.1;
y0(2) = 0.1;
y0(5) = 0.1;
y0(9) = 0.1;
y0(10) = 0.1;
% w=1*rand(28,1);
w(4) = 0.5; % transcription binding affinity
w(5) = 0.2; % reversion transcription binding affinity
w(3) = 15; % transcription rate
w(12) = 0.1; % translation binding affinity
w(13) = 0.2; % reverse translation binding affinity
w(14) = 1; % translation rate
w(21) = 0.1; % inducer binding affinity
w(20) = 0.5; % reverse inducer binding affinity
w(23) = 1; % inducer makes repressor unbind the promoter
w(25) = 0.002; % protein degradation
w(28) = 0.01; % mRNA degradation
w(17) = 0.1; % repressor binding affinity
w(18) = 0.05; % reverse repressor binding affinity
w(15) = 0.1; % dimerization rate
w(16) = 0.05; % reverse dimerization rate
w(1) = w(4); w(2) = w(5); w(6) = w(3); w(9) = w(12); w(10) = w(13); w(11) = w(14); w(18) = w(20); w(19) = w(21); w(26) = w(25); w(24) = w(25); w(27) = w(28);
% w = w + 0.5*rand(28,1);
% w(1) = 0.1;
[t,y] = ode45(@(t,y) dydt(t,y,S,w), [0 1300], y0);
% plot(t,y)
plot(t/60,y(:,8))
% hold on
for mm = 0:0.1:1
  w(21) =mm; %4
  [t,y] = ode45(@(t,y) dydt(t,y,S,w), [0 50], y0);
  plot(t,y(:,7))
  hold on
end
% plot(t,y)

% while y(end,7) > 1 || y(end/2,7) <
%   w=rand(28,1);
%   [t,y] = ode45(@(t,y) dydt(t,y,S,w), [0 120], y0);
% end
% plot(t,y(:,7))

% 0.5752
% 0.8676
% 0.9695
% 0.5564
% 0.3430
% 0.7485
% 0.9355
% 0.1917
% 0.3283
% 0.1171
% 0.3652
% 0.0661
% 0.2062
% 0.1926
% 0.4291
% 0.1069
% 0.1304
% 0.0855
% 0.0328
% 0.4484
% 0.6831
% 0.9545
% 0.3369
% 0.5863
% 0.9830
% 0.5024
% 0.6987
% 0.6699