SimBiology 3.1

Stochastic Simulation of Radioactive Decay

Stochastically simulates the following model:

This model can also be used to represent irreversible isomerization

This example uses parameters and conditions as described in Daniel T. Gillespie, 1977, "Exact Stochastic Simulation of Coupled Chemical Reactions," The Journal of Physical Chemistry, vol. 81, no. 25, pp. 2340-2361.

Contents

Read the Radioactive Decay Model Saved in SBML Format

SBML = Systems Biology Markup Language, www.sbml.org

model  = sbmlimport('radiodecay.xml')

   SimBiology Model - RadioactiveDecay

   Model Components:
     Compartments:      1
     Events:            0
     Parameters:        1
     Reactions:         1
     Rules:             0
     Species:           2

View Species Objects of the Model

model.Species

   SimBiology Species Array

   Index:    Compartment:    Name:    InitialAmount:    InitialAmountUnits:
   1         unnamed         x        1000              molecule
   2         unnamed         z        0                 molecule

View Reaction Objects of the Model

model.Reactions

   SimBiology Reaction Array

   Index:    Reaction:
   1         x -> z

View Parameter Objects for the Kinetic Law

model.Reactions(1).KineticLaw(1).Parameters

   SimBiology Parameter Array

   Index:    Name:    Value:    ValueUnits:
   1         c        0.5       1/second

Update the Reaction to use MassAction Kinetic Law for Stochastic Solvers.

model.Reactions(1).KineticLaw(1).KineticLawName = 'MassAction';
model.Reactions(1).KineticLaw(1).ParameterVariableNames = {'c'};

Simulate the Model Using the Stochastic (SSA) Solver & Plot

StopTime == 'inf' used to instruct the stochastic solver to run until there were no more reactions to execute (no more reactants), but this behavior changed, so now we set an explicit StopTime (g469950).

cs = model.getconfigset('active');
cs.SolverType = 'ssa';
cs.StopTime = 14.0;
cs.CompileOptions.DimensionalAnalysis = false;
[t,X] = sbiosimulate(model);

plot(t,X);
legend('x', 'z');
title('Stochastic Radioactive Decay Simulation');
ylabel('Number of molecules');

Repeat the Simulation to Show Run-to-Run Variability

title('Multiple Stochastic Radioactive Decay Simulations');
hold on;
for loop = 1:20
    [t,X] = sbiosimulate(model);
    plot(t,X);    % Just plot number of reactant molecules
    drawnow;
end

Overlay the Reaction's ODE Solution in Red

cs = model.getconfigset('active');
cs.SolverType = 'ode15s';
cs.StopTime = 20;
[t,X] = sbiosimulate(model);
plot(t,X,'red');
hold off;