Accelerating the pace of engineering and science

# Documentation

## Stochastic Simulation of Radioactive Decay

This example shows how to build and simulate a model using the SSA stochastic solver.

The following model will be constructed and stochastically simulated:

• Reaction 1: x -> z with a first-order reaction rate, c = 0.5.

• Initial conditions: x = 1000 molecules, z = 0.

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.

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

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

plot(t,X);
legend('x', 'z');
ylabel('Number of molecules');
xlabel('Time (seconds)');
```

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 = getconfigset(model,'active');
cs.SolverType = 'sundials';
cs.StopTime = 20;
[t,X] = sbiosimulate(model);
plot(t,X,'red');
hold off;
```