Введение

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

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

Загрузка модели

Мы будем использовать модель G-белка, которая была построена с использованием опубликованных данных Yi et al. (2003). Загрузите модель G-белка дикого типа и посмотрите на его виды и реакции.

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               

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

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

Измените тип решателя на 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 почти идентичны. Если вы посмотрите на данные в matrix 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.