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;