Введение

Когда поведение модели является стохастическим по своей природе, одна запущенная симуляция не обеспечивает достаточно понимания модели. Нужно выполнить ансамбль запусков. Запуски ансамбля производят большие объемы данных, которые требуют систематического анализа.

Этот пример иллюстрирует, как сделать запуски ансамбля с помощью SimBiology и как анализировать сгенерированные данные.

Модель загрузки

Мы будем использовать модель G Protein, которая была создана с помощью опубликованных данных от И и др. (2003). Загрузите модель G Protein дикого типа и посмотрите на ее разновидности и реакции.

sbioloadproject gprotein_norules m1
m1.Species

ans = 
   SimBiology Species Array

   Index:    Compartment:    Name:    Value:       Units:
   1         unnamed         G        7000               
   2         unnamed         Gd       3000               
   3         unnamed         Ga       0                  
   4         unnamed         RL       0                  
   5         unnamed         L        6.022e+17          
   6         unnamed         R        10000              
   7         unnamed         Gbg      3000               

m1.Reactions

ans = 
   SimBiology Reaction Array

   Index:    Reaction:              
   1         L + R <-> RL           
   2         R <-> null             
   3         RL -> null             
   4         Gd + Gbg -> G          
   5         G + RL -> Ga + Gbg + RL
   6         Ga -> Gd               

Выполните запуски ансамбля

Запуски ансамбля могут быть сделаны, только если вы используете стохастические решатели. С детерминированными решателями это не целесообразно делать запуски ансамбля, поскольку вы будете всегда получать точно идентичные результаты. Из-за стохастичности, демонстрационные графики на этой странице не могут точно совпадать с графиками, которые вы получаете от выполнения этого примера сами.

Измените тип решателя в SSA и выполните 10 запусков модели. Чтобы ускорить симуляцию, давайте регистрировать только каждое 100-е событие реакции.

numRuns = 10;
configsetObj = getconfigset(m1,'active');
configsetObj.SolverType = 'ssa';
configsetObj.SolverOptions.LogDecimation = 100;
simdata = sbioensemblerun(m1, numRuns);

График необработанных данных

Отобразите на графике необработанные данные, соответствующие разновидностям по имени G, Ga и RL, и аннотируйте график соответственно.

figure;
speciesNames = {'G', 'Ga', 'RL'};
[t, x] = selectbyname(simdata, speciesNames);
hold on;
for i = 1:numRuns
    plot(t{i}, x{i});
end

grid on;
xlabel('Time in seconds');
ylabel('Species amount');
title('Raw Ensemble Data');
legend(speciesNames);

Статистика ансамбля

Давайте вычислим статистику ансамбля для этих разновидностей.

[timeVals, meanVals, varianceVals] = sbioensemblestats(simdata, speciesNames);

Давайте построим среднее и стандартное отклонение как функцию времени.

figure;
plot(timeVals, meanVals);
grid on;
xlabel('Time in seconds');
ylabel('Mean');
title('Variation of Mean');
legend(speciesNames, 'Location', 'NorthEastOutside');

figure;
plot(timeVals, sqrt(varianceVals));
grid on;
xlabel('Time in seconds');
ylabel('Standard Deviation');
title('Variation of Standard Deviation');
legend(speciesNames, 'Location', 'NorthEastOutside');

Из этих графиков кажется, что G и Ga имеют точно идентичные стандартные отклонения, хотя изменение их средних значений отличается. Давайте смотреть, можем ли мы узнать, почему это так. Первый шаг должен был бы выяснить, существует ли отношение между ними. Давайте посмотрим на реакции в модели еще раз, чтобы видеть, очевидно ли что-нибудь.

m1.Reactions

ans = 
   SimBiology Reaction Array

   Index:    Reaction:              
   1         L + R <-> RL           
   2         R <-> null             
   3         RL -> null             
   4         Gd + Gbg -> G          
   5         G + RL -> Ga + Gbg + RL
   6         Ga -> Gd               

Сохранение половины

От реакций, отношения между G и Ga не совсем ясно. Чтобы анализировать далее, мы собираемся использовать sbioconsmoiety видеть, существует ли какое-либо сохранение половины.

sbioconsmoiety(m1, 'semipos', 'p')

ans = 2x1 cell
    {'G + Gbg'    }
    {'G + Gd + Ga'}

От этого мы видим, что 'G+ Gd + Ga' всегда сохраняется. Но это не полностью объясняет, почему отклонения G и Ga идентичны. Что относительно Gd? Почему его отклонение не влияет на G или Ga? Чтобы изучить его далее, давайте вычислим статистику ансамбля для Gd и давайте построим их изменение.

[timeValsGd, meanValsGd, varianceValsGd] = sbioensemblestats(simdata, 'Gd');
figure;
plot(timeValsGd, meanValsGd, '-', ...
    timeValsGd, sqrt(varianceValsGd), 'r:');
axis([-50 600 -500 3000]);
xlabel('Time in seconds');
title('Mean and Standard Deviation of Gd');
legend('Mean', 'Standard Deviation')

Объяснение идентичных отклонений

Из графика это видно, который Gd начинает с ненулевым значением во время = 0, но и его среднее значение и нуль подхода отклонения очень резко и остаются там. Таким образом, когда Gd остается около нуля, уравнение сохранения половины уменьшает до

Это означает, когда G повышается, Ga спускается равной суммой и наоборот. Это объясняет, почему отклонения G и Ga почти идентичны. Если вы посмотрите на данные в матрице varianceVals, вы будете видеть, что эти два отклонения очень близки, но не точно равны. Это происходит из-за присутствия Gd, который является очень близко к нулю, но не точно нулевой.

Графики ансамбля: 2D графики распределения

Один способ визуализировать стохастические данные ансамбля состоит в том, чтобы построить гистограммы концентраций разновидностей в точках определенного времени. Каждая гистограмма показывает распределение концентраций для конкретной разновидности по целому ансамблю запусков. Эти гистограммы могут быть сгенерированы с помощью команды SimBiology sbioensembleplot.

Давайте создадим гистограммы для разновидностей G, Ga и RL во время t = 10. Обратите внимание на то, что в этом примере, мы сгенерировали данные ансамбля, не задавая опцию интерполяции для sbioensemblerun. Временные векторы для каждого запуска в ансамбле поэтому отличаются друг от друга. sbioensembleplot интерполирует данные моделирования, чтобы найти суммы разновидностей для каждого запуска в точное время t = 10.

sbioensembleplot(simdata, speciesNames, 10);

Графики ансамбля: 3D горные графики

Ясно, что, чтобы получить полное понимание того, как распределение изменяется со временем, необходимо будет сделать эти графики распределения в каждом временном интервале интереса. Кроме того, графики распределения должны быть замечены наряду с графиками отклонения и средним значением. Было бы хорошо соединить всю эту информацию во всего одном графике.

Ну, 3D графики данных ансамбля делают точно это. С этими 3D графиками можно просмотреть, как средний и отклонение изменяются как функция времени. Кроме того, вместо того, чтобы иметь необходимость построить распределение разновидностей на каждом возможном временном шаге, в одном представлении вы видите, как подходящее нормальное распределение, с тем же средним значением и отклонением как фактические данные, изменяется со временем. 3D график ансамбля превосходен для того, чтобы получить представление о том, как средний, отклонение и распределение варьируются как функция времени.

Давайте см. 3D график ансамбля разновидностей Ga.

sbioensembleplot(simdata, 'Ga');

Этот пример показал, как стохастические данные ансамбля моделей SimBiology могут анализироваться с помощью различных инструментов в SimBiology.