В этом примере показано, как построить и смоделировать модель с помощью стохастического решателя SSA.
Будет построена и стохастически смоделирована следующая модель:
Реакция 1: x + y1 - > 2 y1 + x, с постоянной скоростью, c1 = 10.
Реакция 2: y1 + y2 - > 2 y2, с постоянной скоростью, c2 = 0,01.
Реакция 3: y2 - > z, с постоянной скоростью, c3 = 10.
Исходные условия: x = 1 (константа), y1 = y2 = 1000, z = 0.
Примечание: Вид 'x' в реакции 1 представлен с обеих сторон реакции для моделирования предположения, что количество x является постоянным.
Эти реакции можно интерпретировать как простую модель хищника-добычи, если считать, что популяция добычи (y1) увеличивается в присутствии пищи (x) (реакция 1), что популяция хищника (y2) увеличивается по мере того, как они едят добычу (реакция 2), и что хищники (y2) умирают по естественным причинам (реакция 3).
В этом примере используются параметры и условия, описанные в Daniel T. Gillespie, 1977, «Точное стохастическое моделирование связанных химических реакций», The Journal of Physical Chemistry, vol. 81, no. 25, pp. 2340-2361.
sbioaddtolibrary(sbiounit('rabbit','molecule',1)); sbioaddtolibrary(sbiounit('coyote','molecule',1)); sbioaddtolibrary(sbiounit('food','molecule',1)); sbioaddtolibrary(sbiounit('amountDimension','molecule',1));
model = sbiomodel('Lotka-Volterra Model'); c = addcompartment(model,'C'); c.CapacityUnits = 'meter^3';
r1 = addreaction(model,'x + y1 -> 2 y1 + x')r1 = SimBiology Reaction Array Index: Reaction: 1 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
r2 = addreaction(model,'y1 + y2 -> 2 y2')r2 = SimBiology Reaction Array Index: Reaction: 1 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
r3 = addreaction(model,'y2 -> z')r3 = SimBiology Reaction Array Index: Reaction: 1 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';
model
model =
SimBiology Model - Lotka-Volterra Model
Model Components:
Compartments: 1
Events: 0
Parameters: 3
Reactions: 3
Rules: 0
Species: 4
Observables: 0
model.Reactions
ans = SimBiology Reaction Array Index: Reaction: 1 x + y1 -> 2 y1 + x 2 y1 + y2 -> 2 y2 3 y2 -> z
model.Species
ans = SimBiology Species Array Index: Compartment: Name: Value: Units: 1 C x 1 food 2 C y1 1000 rabbit 3 C y2 1000 coyote 4 C z 0 amountDimension
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;
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');