www.gusucode.com > simbio 案例源码程序 matlab代码 > simbio/lotkavolterra.m
%% Stochastic Simulation of the Lotka-Volterra Reactions % % 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 + y1 -> 2 y1 + x, with rate constant, c1 = 10. % * Reaction 2: y1 + y2 -> 2 y2, with rate constant, c2 = 0.01. % * Reaction 3: y2 -> z, with rate constant, c3 = 10. % * Initial conditions: x=1 (constant), y1=y2=1000, z=0. % * Note: Species 'x' in Reaction 1 is represented on both sides of the % reaction to model the assumption that the amount of x is % constant. % % These reactions can be interpreted as a simple predator-prey model if % one considers that the prey population (y1) increases in the presence % of food (x) (Reaction 1), that the predator population (y2) increases % as they eat prey (Reaction 2), and that predators (y2) % die of natural causes (Reaction 3). % % 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. %% Register Units for the Model sbioaddtolibrary(sbiounit('rabbit','molecule',1)); sbioaddtolibrary(sbiounit('coyote','molecule',1)); sbioaddtolibrary(sbiounit('food','molecule',1)); sbioaddtolibrary(sbiounit('amountDimension','molecule',1)); %% Create the Lotka-Volterra Model model = sbiomodel('Lotka-Volterra Model'); c = addcompartment(model,'C'); c.CapacityUnits = 'liter'; %% Add Reaction 1 to the Model Object r1 = addreaction(model,'x + y1 -> 2 y1 + x') % Set the Kinetic Law for Reaction 1. kl1 = addkineticlaw(r1, 'MassAction'); % Add rate constant parameter, c1, to reaction with value = 10 p1 = addparameter(kl1, 'c1', 'Value', 10); kl1.ParameterVariableNames = {'c1'}; % Add units to c1 p1.ValueUnits = '1/(second*rabbit)'; % Set initial amounts for species in Reaction 1 r1.Reactants(1).InitialAmount = 1; % x r1.Reactants(2).InitialAmount = 1000; % y1 % Set the initial amount units for species in Reaction 1 r1.Reactants(1).InitialAmountUnits = 'food'; % x r1.Reactants(2).InitialAmountUnits = 'rabbit'; % y1 %% Add Reaction 2 to the Model Object r2 = addreaction(model,'y1 + y2 -> 2 y2') % Set the kinetic law for Reaction 2. kl2 = addkineticlaw(r2, 'MassAction'); % Add rate constant parameter, c2, to kinetic law with value = 0.01 p2 = addparameter(kl2, 'c2', 'Value', 0.01); kl2.ParameterVariableNames = {'c2'}; % Add units to c2 p2.ValueUnits = '1/(second*coyote)'; % Set initial amounts for new species in Reaction 2 r2.Products(1).InitialAmount = 1000; % y2 % Set the initial amount units for new species in Reaction 2 r2.Products(1).InitialAmountUnits = 'coyote'; % y2 %% Add Reaction 3 to the Model Object r3 = addreaction(model,'y2 -> z') % Add "bogus" units to trash variable 'z' r3.Products(1).InitialAmountUnits = 'amountDimension'; % Set the kinetic law for Reaction 3. kl3 = addkineticlaw(r3, 'MassAction'); % Add rate constant parameter, c3, to reaction with value = 10 p3 = addparameter(kl3, 'c3', 'Value', 10); kl3.ParameterVariableNames = {'c3'}; % Add units to c3 p3.ValueUnits = '1/second'; %% Display the Completed Model Objects model %% Display the Reaction Objects model.Reactions %% Display the Species Objects model.Species %% Simulate with the Stochastic (SSA) Solver & Plot cs = getconfigset(model,'active'); cs.SolverType = 'ssa'; cs.StopTime = 30; cs.SolverOptions.LogDecimation = 200; cs.CompileOptions.UnitConversion = true; [t,X] = sbiosimulate(model); plot(t, X(:,2), t, X(:,3)); legend('Y1', 'Y2'); title('Lotka-Volterra Reaction - State History'); ylabel('Number of predator-prey'); grid on; %% Show Phase Portrait of Y1 to Y2 plot(X(:,2),X(:,3)); title('Lotka-Volterra Reaction - Y1 vs. Y2'); xlabel('Number of Y1 rabbits'); ylabel('Number of Y2 coyotes'); % Clean up units. sbioremovefromlibrary('unit', 'rabbit'); sbioremovefromlibrary('unit', 'coyote'); sbioremovefromlibrary('unit', 'food'); sbioremovefromlibrary('unit', 'amountDimension');