В этом примере показано, как создать и симулировать модель с помощью решателя стохастики 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, «Точная стохастическая симуляция связанных химических реакций», Журнал физической химии, том 81, № 25, стр. 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');