Trying to model the dynamics of the non-regulated construct and the auto-regulated construct using the Gillespie algorithm, do you know why it doesn't reach 200 and stay?
1 view (last 30 days)
Show older comments
t_max= 1200; % in s
k=0.01; % in s^-1
g_nr=2; % in molecules/s
t_nr=0; % in s
x_nr=5; % in molecules
step_counter=1;
% Simulation of the non-regulated construct
t=0;
prod_rate=g_nr;
deg_rate=k*x_nr;
while t<t_max
wait_time=-log(rand)/((prod_rate)+(deg_rate*x_nr(step_counter)));
prob_prod=(prod_rate)/((prod_rate)+(deg_rate*x_nr(step_counter))); % Propensity of production
t=t+wait_time; % Update current time
step_counter=step_counter+1; % Update the number of steps (reactions) associated with the experiment.
t_nr(step_counter)=t; % Add the current time to the time log
if rand < prob_prod % Defines whether production takes place based on Monte Carlo method.
x_nr(step_counter)=x_nr(step_counter-1)+1; % Implements production
else
x_nr(step_counter)=x_nr(step_counter-1)-1; % Implements degradation
end
if t>t_max
t_nr(end)=[];
x_nr(end)=[];
end
end
close all
figure,
plot(t_nr,'r','LineWidth',2)
hold on
plot(x_nr,'b','LineWidth',2)
legend('t_nr','x_nr')
x_nr should oscillate around a value of 200 within ~500 s after the start of the experiment, but my code only goes up to approximately 90. The production and degradation rate as well as the probability of production are all given as equations that I'm sure I got right, so I'm not sure what's causing it to only go up to 90 molecules.
0 Comments
Answers (1)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!