Сравнение SSA и явных стохастических решателей с скачками Тау

В этом примере показано, как создать и симулировать модель с помощью решателя SSA и решателя Explicit Tau-Leaping.

Будут построены и стохастически моделированы следующие реакции димеризации распада:

Реакция 1: s1 - > null, с постоянной скоростью реакции, c1 = 1,0
Реакция 2:2 s1 < - > s2, с постоянными скоростью реакции, вперед: c2f = 0,002 reverse: c2r = 0,5
Реакция 3: s2 - > s3, с постоянной скоростью реакции, c3 = 0,04
Начальные условия: s1 = 100000 молекул, s2 = s3 = 0

Этот пример использует параметры и условия, как описано в Daniel T. Gillespie, 2001, «Приблизительная ускоренная стохастическая симуляция химически взаимодействующих систем», Journal of Chemical Physics, vol. 115, no. 4, pp. 1716-1733.

Создайте модель затухания-димеризации

model = sbiomodel('Decaying-Dimerizing Reaction Set');

Ввод реакций

r1 = addreaction(model, 's1 -> null');
r2 = addreaction(model, '2 s1 <-> s2');
r3 = addreaction(model, 's2 -> s3');

Установите реакции в MassAction

kl1 = addkineticlaw(r1, 'MassAction');
kl2 = addkineticlaw(r2, 'MassAction');
kl3 = addkineticlaw(r3, 'MassAction');

Добавьте константы скорости для каждой реакции

p1  = addparameter(kl1, 'c1',  'Value', 1.0);
p2f = addparameter(kl2, 'c2f', 'Value', 0.002);
p2r = addparameter(kl2, 'c2r', 'Value', 0.5);
p3  = addparameter(kl3, 'c3',  'Value', 0.04);

Установите константы кинетического закона для каждого кинетического закона.

kl1.ParameterVariableNames = {'c1'};
kl2.ParameterVariableNames = {'c2f', 'c2r'};
kl3.ParameterVariableNames = {'c3'};

Задайте начальное количество каждого вида

model.species(1).InitialAmount = 100000;    % s1
model.species(2).InitialAmount = 0;         % s2
model.species(3).InitialAmount = 0;         % s3

Отображение завершенных объектов модели

model

model = 
   SimBiology Model - Decaying-Dimerizing Reaction Set 

   Model Components:
     Compartments:      1
     Events:            0
     Parameters:        4
     Reactions:         3
     Rules:             0
     Species:           3
     Observables:       0

Отображение объектов реакции

model.Reactions

ans = 
   SimBiology Reaction Array

   Index:    Reaction:  
   1         s1 -> null 
   2         2 s1 <-> s2
   3         s2 -> s3

Отображение видовых объектов

model.Species

ans = 
   SimBiology Species Array

   Index:    Compartment:    Name:    Value:    Units:
   1         unnamed         s1       100000          
   2         unnamed         s2       0               
   3         unnamed         s3       0

Получите активную конфигурацию модели для модели.

cs = getconfigset(model,'active');

Симулируйте модель с помощью решателя SSA Stochastic и постройте график

tfinal = 30, регистрация каждой 10-й точки данных.

cs.SolverType = 'ssa';
cs.StopTime = 30;
solver = cs.SolverOptions;
solver.LogDecimation = 10;
cs.CompileOptions.DimensionalAnalysis = false;
[t_ssa, x_ssa] = sbiosimulate(model);

h1 = subplot(2,1,1);
plot(h1, t_ssa, x_ssa(:,1),'.');
h2 = subplot(2,1,2);
plot(h2, t_ssa, x_ssa(:,2:3),'.');
grid(h1,'on'); 
grid(h2,'on');
title(h1,'Decay Dimerizing Reactions');
ylabel(h1,'Amount of S1');
ylabel(h2,'Amount of S2 & S3');
xlabel(h2,'Time');
legend(h2, 'S2', 'S3');

Figure contains 2 axes. Axes 1 with title Decay Dimerizing Reactions contains an object of type line. Axes 2 contains 2 objects of type line. These objects represent S2, S3.

Симулируйте модель с помощью явного решателя Tau-Leaping и постройте график на той же фигуре

Не закрывая окно рисунка, постройте график результатов от использования Явного Tau-Leaping Solver.

tfinal = 30, регистрация каждой 10-й точки данных. Допустимый допуск ошибок для решателя, 0,03.

cs.StopTime = 30;
cs.SolverType = 'explTau';
solver = cs.SolverOptions;
solver.LogDecimation = 10;
[t_etl, x_etl] = sbiosimulate(model);

hold(h1,'on'); 
hold(h2,'on');
plot(h1, t_etl, x_etl(:,1),'o');
plot(h2, t_etl, x_etl(:,2:3),'o');
legend(h2, 'S2 (SSA)', 'S3 (SSA)', 'S2 (Exp. Tau)', 'S3 (Exp. Tau)');
hold(h1,'off'); 
hold(h2,'off');

Figure contains 2 axes. Axes 1 with title Decay Dimerizing Reactions contains 2 objects of type line. Axes 2 contains 4 objects of type line. These objects represent S2 (SSA), S3 (SSA), S2 (Exp. Tau), S3 (Exp. Tau).

Сравнение количества шагов для SSA и явных алгоритмов скачка Тау

fprintf('Approximate Number of SSA steps: %d\n', (length(t_ssa) * 10));

Approximate Number of SSA steps: 616010

fprintf('Approximate Number of Explicit Tau-Leaping steps: %d\n', ...
        (length(t_etl) * 10));

Approximate Number of Explicit Tau-Leaping steps: 620

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