Clear Filters
Clear Filters

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)
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.

Answers (1)

Jeffrey Clark
Jeffrey Clark on 30 Jun 2022
@Jonathan Louie, I believe your deg_rate=k*x_nr should just be deg_rate=k.

Categories

Find more on Scripts in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!