Введение

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

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

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

Мы будем использовать модель G Protein, которая была построена с использованием опубликованных данных 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               

Выполнение ансамблевых прогонов

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

Измените тип решателя на 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 идентичны. А как же Гд? Почему его дисперсия не влияет на 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 почти идентичны. Если вы посмотрите на данные в матричном варьировании Vals, вы увидите, что два отклонения очень близки, но не точно равны. Это связано с присутствием Gd, которое очень близко к нулю, но не совсем равно нулю.

Ансамблевые участки: 2D Распределительные участки

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

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

sbioensembleplot(simdata, speciesNames, 10);

Ансамблевые участки: 3D Горные участки

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

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

Посмотрим 3D ансамблевый сюжет вида Га.

sbioensembleplot(simdata, 'Ga');

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