2

I want to run stochastic simulations of a model which has a rate constant that is time dependent. I am not aware of any stochastic simulator that supports time-dependent rate constants.

So far I have tried to use Python (stochpy and gillespy2) and Matlab (Simbiology).

I am trying to do this for the simple model written below (code runs in Matlab (Simbiology)). My question is: in this model how do I make k1 time dependent, for example making k1 changing linearly with time during the stochastic simulation?

% Create a new model
model = sbiomodel('GeneExpressionModel');

% Add reactions for gene expression
reaction1 = addreaction(model, 'Gene -> Gene + mRNA');
reaction2 = addreaction(model, 'mRNA -> mRNA + Protein');
reaction3 = addreaction(model, 'mRNA -> null');
reaction4 = addreaction(model, 'Protein -> null');

% Set reaction rate constants
kineticLaw1 = addkineticlaw(reaction1, 'MassAction');
kineticLaw2 = addkineticlaw(reaction2, 'MassAction');
kineticLaw3 = addkineticlaw(reaction3, 'MassAction');
kineticLaw4 = addkineticlaw(reaction4, 'MassAction');

% Add parameters for reaction rate constants
p1 = addparameter(kineticLaw1, 'k1', 'Value', 0.1);
p2 = addparameter(kineticLaw2, 'k2', 'Value', 0.05);
p3 = addparameter(kineticLaw3, 'k3', 'Value', 0.01);
p4 = addparameter(kineticLaw4, 'k4', 'Value', 0.002);

% Set the rate equation for each reaction
kineticLaw1.ParameterVariableNames = {p1.Name};
kineticLaw2.ParameterVariableNames = {p2.Name};
kineticLaw3.ParameterVariableNames = {p3.Name};
kineticLaw4.ParameterVariableNames = {p4.Name};

% Set initial conditions
model.Species(1).InitialAmount = 1;
model.Species(2).InitialAmount = 0;
model.Species(3).InitialAmount = 0;

% Configure the stochastic solver
configSet = getconfigset(model, 'active');
set(configSet, 'SolverType', 'ssa'); % Set the solver to stochastic (Gillespie algorithm)
set(configSet, 'StopTime', 100); % Set the simulation end time

% Run the stochastic simulation
simulationData = sbiosimulate(model);

% Plot the simulation results
figure;
plot(simulationData.Time, simulationData.Data);
xlabel('Time');
ylabel('Species Count');
legend('Gene', 'mRNA', 'Protein');
title('Stochastic Simulation of Gene Expression');
Cristina Palma
  • 21
  • 1

1 Answers1

0

This is not possible in SimBiology or any software that implements the Gillespie algorithm. This algorithm requires that the rate constant is, well, constant.

How does your rate constant vary with time? If it's a discrete change, could you do separate simulations for each time period that has a particular constant value? If the change is not discrete, then you will likely need to use some other mathematical approach for modeling this system. However, that's not my area of expertise, so I can't recommend anything specific.

agoldsipe
  • 181
  • 4
  • Thanks agodsipe. My change is continuous but maybe can be approximated by making it discrete with minimal error. I will try that approach. I also don't know any other mathematical approach. I have heard about 'hybrid models' (deterministic + stochastic approach) . Maybe I will explore that. Thanks. – Cristina Palma May 01 '23 at 18:02