Сравнение SSA и явных Tau-прыгающих стохастических решателей

Скрипт Open Live Script

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

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

Реакция 1: s1-> пустой указатель, с постоянной скоростью реакции, c1 = 1.0
Реакция 2: 2 s1 <-> s2, с константами скорости реакции, передают: c2f = 0,002 реверса: c2r = 0.5
Реакция 3: s2-> s3, с постоянной скоростью реакции, c3 = 0.04
Начальные условия: s1 = 100 000 молекул, s2 = s3 = 0

Этот пример использует параметры и условия как описано в Дэниеле Т. Гиллеспи, 2001, "Аппроксимированная ускоренная стохастическая симуляция химически реагирующих систем", Журнал Химической Физики, издания 115, № 4, стр 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 и график

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');

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

Не закрывая окно рисунка, постройте результаты использования Явного Tau-прыгающего Решателя.

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');

Сравнение количества шагов для SSA и явных Tau-прыгающих алгоритмов

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

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