www.gusucode.com > simbio 案例源码程序 matlab代码 > simbio/ensemblerun_data_analysis.m
%% Analysis of Stochastic Ensemble Data in SimBiology(R) % % This example shows how to make ensemble runs and how to analyze the % generated data in SimBiology(R). % % Copyright 2004-2014 The MathWorks, Inc. %% Introduction % % When the behavior of a model is stochastic in nature, a single simulation % run does not provide enough insight into the model. One has to perform an % ensemble of runs. Ensemble runs produce large amounts of data that % require systematic analysis. % % This example illustrates how to make ensemble runs using SimBiology and % how to analyze the generated data. % %% Load Model % % We will use the G Protein model that was built using published data from % Yi et al. (2003). Load the wild-type G Protein model and look at its % species and reactions. sbioloadproject gprotein_norules m1 m1.Species m1.Reactions %% Perform Ensemble Runs % % Ensemble runs can be done only if you are using stochastic solvers. With % deterministic solvers, it does not make sense to do ensemble runs since % you will always get exactly identical results. Due to stochasticity, % the sample plots on this page might not exactly match the plots you get % from running this example yourself. %% % Change the solver type to SSA and perform 10 runs of the model. To % speed up the simulation, let's log only every 100th reaction event. % numRuns = 10; configsetObj = getconfigset(m1,'active'); configsetObj.SolverType = 'ssa'; configsetObj.SolverOptions.LogDecimation = 100; simdata = sbioensemblerun(m1, numRuns); %% Plot of Raw Data % Plot the raw data corresponding to the species named G, Ga and RL and % annotate the plot appropriately. 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); %% Ensemble Statistics % Let's compute ensemble statistics for these species. [timeVals, meanVals, varianceVals] = sbioensemblestats(simdata, speciesNames); %% % Let's plot the mean and standard deviation as a function of time. 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'); %% % From these plots it appears that G and Ga have exactly identical standard % deviations, though the variation of their means is different. Let's see % if we can find out why that is so. The first step would be to figure out % if there is any relationship between them. Let's look at the reactions in % the model once again to see if anything is obvious. m1.Reactions %% Moiety Conservation % From the reactions, the relationship between G and Ga is not quite clear. % To analyze further we are going to use |sbioconsmoiety| to see if % there is any moiety conservation. sbioconsmoiety(m1, 'semipos', 'p') %% % From this we can see that 'G + Gd + Ga' is always conserved. But that % does not completely explain why the variances of G and Ga are identical. % What about Gd? Why doesn't its variance affect G or Ga? To look into it % further, let's compute the ensemble statistics for Gd and plot their % variation. [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') %% Explanation of Identical Variances % From the plot, it can be seen that Gd starts out with a non-zero value at % time = 0, but both its mean and variance approach zero very sharply and % stay there. Thus, when Gd stays near zero, the moiety conservation % equation reduces to % % $$G + Ga \approx constant$$ % % This means as G goes up Ga goes down by an equal amount and vice versa. % This explains why the variances of G and Ga are almost identical. If you % look at the data in matrix varianceVals, you'll see that the two % variances are very close but not exactly equal. This is due to the % presence of Gd which is very close to zero but not exactly zero. % %% Ensemble Plots: 2D Distribution Plots % % One way to visualize stochastic ensemble data is to plot histograms of % species concentrations at particular time points. Each histogram shows % the distribution of concentrations for a particular species over the % entire ensemble of runs. These histograms may be generated using the % SimBiology command |sbioensembleplot|. % % Let's create histograms for the species G, Ga, and RL at time t = 10. % Note that in this example, we generated the ensemble data without % specifying an interpolation option for |sbioensemblerun|. The time % vectors for each run within the ensemble are therefore different from % each other. |sbioensembleplot| interpolates the simulation data to find % species amounts for every run at the precise time t = 10. % sbioensembleplot(simdata, speciesNames, 10); %% Ensemble Plots: 3D Mountain Plots % It is clear that to get a full understanding of how the distribution % changes with time, you will need to make these distribution plots at % every time interval of interest. Moreover, the distribution plots need to % be seen along with the mean and variance plots. It would be nice to put % all of this information together in just one plot. % % Well, the 3D ensemble data plots do exactly that. With these 3D plots you % can view how mean and variance change as a function of time. In addition, % instead of having to plot the distribution of species at every possible % time step, in one view you can see how a fitted normal distribution, with % the same mean and variance as the actual data, changes with time. The 3D % ensemble plot is excellent for getting an overview of how mean, variance % and distribution vary as a function of time. % % Let's see a 3D ensemble plot of species Ga. % sbioensembleplot(simdata, 'Ga'); %% % This example showed how stochastic ensemble data of SimBiology models can % be analyzed using various tools in SimBiology. %