Стохастическая симуляция реакций Лотки - Вольтерры

В этом примере показано, как создать и симулировать модель с помощью стохастического решателя 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');

Документация