Стохастическая симуляция реакций Лотки - Вольтерры
В этом примере показано, как создать и симулировать модель с помощью стохастического решателя 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).
Этот пример использует параметры и условия как описано в Дэниеле Т. Гиллеспи, 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';
Добавьте реакцию 1 в объект модели
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
Добавьте реакцию 2 в объект модели
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
Добавьте реакцию 3 в объект модели
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
Симулируйте со стохастической Solver & Plot (SSA)
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 Y1 к 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');