www.gusucode.com > simbio 案例源码程序 matlab代码 > simbio/radiodecay.m

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

% Copyright 2004-2014 The MathWorks, Inc.


%% Read the Radioactive Decay Model Saved in SBML Format
% SBML = Systems Biology Markup Language, <http://www.sbml.org www.sbml.org>

model  = sbmlimport('radiodecay.xml')

%% View Species Objects of the Model
model.Species


%% View Reaction Objects of the Model
model.Reactions


%% View Parameter Objects for the Kinetic Law
model.Reactions(1).KineticLaw(1).Parameters

%% 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');
title('Stochastic Radioactive Decay Simulation');
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;