diff --git a/Examples/example_ExtrinsicNoiseModels.m b/Examples/example_ExtrinsicNoiseModels.m index 909691e..de488df 100644 --- a/Examples/example_ExtrinsicNoiseModels.m +++ b/Examples/example_ExtrinsicNoiseModels.m @@ -29,4 +29,5 @@ 2+0.4*randn,... 1+0.2*randn]; extrinsicModel = extrinsicSSIT(Model1,parDistributions,20); -Model1.makePlot(extrinsicModel.averagedResults,'marginals',[],[],[2,3]) \ No newline at end of file +Model1.makePlot(extrinsicModel.averagedResults,'marginals',[],[],[2,3]) +legend('Intrinsic Only','Intrinsic + Extrinsic','Location','EastOutside') \ No newline at end of file diff --git a/Examples/example_FIMCalculation.m b/Examples/example_FIMCalculation.m index 3fb610a..054ac45 100644 --- a/Examples/example_FIMCalculation.m +++ b/Examples/example_FIMCalculation.m @@ -13,21 +13,35 @@ %% (2) Solve FSP for model F2.solutionScheme = 'FSP'; % Set solution scheme to FSP. [FSPsoln,F2.fspOptions.bounds] = F2.solve; % Solve the FSP analysis +% Plot the results from the FSP analysis +fig1 = figure(1);clf; set(fig1,'Name','Marginal Distributions, mRNA'); +fig2 = figure(2);clf; set(fig2,'Name','Marginal Distributions, Protein'); +F2.makePlot(FSPsoln,'marginals',[1:4:21],false,[fig1,fig2],{'r-','linewidth',2}) %% (3) Solve FSP Sensitivity F2.solutionScheme = 'fspSens'; % Set solutions scheme to FSP Sensitivity [sensSoln,bounds] = F2.solve(FSPsoln.stateSpace); % Solve the sensitivity problem +% Plot the results from the sensitivity analysis +fig3 = figure(3);clf; set(fig3,'Name','Marginal Sensitivity, mRNA'); +fig4 = figure(4);clf; set(fig4,'Name','Marginal Sensitivity, Protein'); +F2.makePlot(sensSoln,'marginals',[],false,[fig3,fig4],{'b','linewidth',2}) %% (4) Compute FIM using FSP Sensitivity Results fimResults = F2.computeFIM(sensSoln.sens); % Compute the FIM for full observations and no distortion. cellCounts = 10*ones(size(F2.tSpan)); % Number of cells in each experiment. [fimTotal,mleCovEstimate,fimMetrics] = F2.evaluateExperiment(fimResults,cellCounts) +fig5 = figure(5);clf; set(fig5,'Name','Fim-Predicted Uncertainty Ellipses'); +F2.plotMHResults([],fimTotal,'lin',[],fig5) +legend('FIM - Full Observation') %% (5) Compute FIM for Partial Observations F2.pdoOptions.PDO=[]; F2.pdoOptions.unobservedSpecies = 'x1'; [fimResults_partialObs] = F2.computeFIM(sensSoln.sens); % Compute the FIM for full observations and no distortion. [fimTotal_partialObs,mleCovEstimate_partialObs,fimMetrics_partialObs] = F2.evaluateExperiment(fimResults_partialObs,cellCounts) +fig6 = figure(6);clf; set(fig6,'Name','Fim-Predicted Uncertainty Ellipses'); +F2.plotMHResults([],[fimTotal,fimTotal_partialObs],'lin',[],fig6) +legend('FIM - Full Observation','FIM - Protein Only') %% (6) Compute FIM for Distorted Observation (Probabilistic Distortion Operator) F2.pdoOptions.unobservedSpecies = 'x1'; @@ -36,17 +50,20 @@ pdoOptions.props.CaptureProbabilityS1 = 0; % Use zero for unobserved species. pdoOptions.props.CaptureProbabilityS2 = 0.9; -% call method to generate the PDO. +% Call method to generate the PDO. F2.pdoOptions.PDO = F2.generatePDO(pdoOptions,[],FSPsoln.fsp); -% plot the PDO +% Plot the PDO N = size(F2.pdoOptions.PDO.conditionalPmfs{1}); -figure; +fig7 = figure(7); set(fig7,'Name','Probabilistic Distortion Operator for Protein'); contourf([0:N(1)-1],[0:N(2)-1],log10(F2.pdoOptions.PDO.conditionalPmfs{1})); % Here we wanted the first PDO for 'x2' because 'x1' was unobserved. xlabel('Actual');ylabel('Observable');colorbar; -% solve FIM using the specified PDO +% Solve FIM using the specified PDO [fimResults_BinomialPDO] = F2.computeFIM(sensSoln.sens); % Compute the FIM for full observations and no distortion. [fimTotal_BinomialPDO,mleCovEstimate_BinomialPDO,fimMetrics_BinomialPDO] =... - F2.evaluateExperiment(fimResults_BinomialPDO,cellCounts) \ No newline at end of file + F2.evaluateExperiment(fimResults_BinomialPDO,cellCounts) +fig8 = figure(8);clf; set(fig8,'Name','Fim-Predicted Uncertainty Ellipses'); +F2.plotMHResults([],[fimTotal,fimTotal_partialObs,fimTotal_BinomialPDO],'lin',[],fig8) +legend('FIM - Full Observation','FIM - Protein Only','FIM - Protein with Error') diff --git a/Examples/example_HybridModel.m b/Examples/example_HybridModel.m index 689f91d..66bc330 100644 --- a/Examples/example_HybridModel.m +++ b/Examples/example_HybridModel.m @@ -27,6 +27,7 @@ Model2.hybridOptions.upstreamODEs = {'rna'}; [fspSoln2, Model2.fspOptions.bounds] = Model2.solve; Model2.makePlot(fspSoln2,'marginals',[],[],3) +legend('Full','QSSA for (n)','Location','eastoutside') %% Example 2 - 5-species MAPK induction Model % In this example, we consider a model of MAPK translocation to the nucleus @@ -71,4 +72,4 @@ Model5.hybridOptions.upstreamODEs = {'mapkCyt','mapkNuc'}; [fspSoln5, Model5.fspOptions.bounds] = Model5.solve; Model5.makePlot(fspSoln5,'marginals',[],[],[11,12,15]) - +legend('Full','QSSA for (gene,MAKP)','QSSA for (MAKP)','Location','eastoutside') diff --git a/Examples/example_IterativeBayesianExperimentDesign.m b/Examples/example_IterativeBayesianExperimentDesign.m index 799f206..1500042 100644 --- a/Examples/example_IterativeBayesianExperimentDesign.m +++ b/Examples/example_IterativeBayesianExperimentDesign.m @@ -34,6 +34,8 @@ fitParameters = [1:4]; end +experimentDesignCriteria = + saveFileName = ['IterativeExperimentResults_',example]; %% Generate Model Propensity Functions and Solve True Model diff --git a/Examples/example_ModelReductionTool.m b/Examples/example_ModelReductionTool.m index 387321e..1f4287c 100644 --- a/Examples/example_ModelReductionTool.m +++ b/Examples/example_ModelReductionTool.m @@ -135,7 +135,10 @@ % Make Figures to compare the results. Here, we will plot the original % model in blue and the reduced model in red lines. Model1.makePlot(fspSoln,'meansAndDevs',[],[],1,{'Color',[0,0,1]}) -Model1.makePlot(fspSoln,'marginals',[],[],[2,3],{'Color',[0,0,1]}) - Model2.makePlot(fspSolnRed,'meansAndDevs',[],[],1,{'Color',[1,0,0]}) -Model2.makePlot(fspSolnRed,'marginals',[],[],[2,3],{'Color',[1,0,0]}) \ No newline at end of file +figure(1);legend('Full','Reduced','Location','southeast') + +Model1.makePlot(fspSoln,'marginals',[],[],[2,3],{'Color',[0,0,1]}) +Model2.makePlot(fspSolnRed,'marginals',[],[],[2,3],{'Color',[1,0,0]}) +figure(2);legend('Full','Reduced','Location','eastoutside') +figure(3);legend('Full','Reduced','Location','eastoutside') \ No newline at end of file diff --git a/Examples/html/example_EscapeTimes.html b/Examples/html/example_EscapeTimes.html index 54e4027..f6a3baf 100644 --- a/Examples/html/example_EscapeTimes.html +++ b/Examples/html/example_EscapeTimes.html @@ -6,7 +6,7 @@ example_EscapeTimes

Contents

example_EscapeTimes

In this script, we demonstrate how to create and solve a fiorst passage time problem.

close all
 clear all
-addpath('../CommandLine');
+addpath(genpath('../src'));
 

Example 1 - a simple transcription/translation model

First create a full model (e.g., for mRNA and protein)

Model1 = SSIT;
 Model1.species = {'rna','protein'};
 Model1.initialCondition = [0;0];
@@ -112,7 +112,7 @@
 % time problem.
 close all
 clear all
-addpath('../CommandLine');
+addpath(genpath('../src'));
 %% Example 1 - a simple transcription/translation model
 % First create a full model (e.g., for mRNA and protein)
 Model1 = SSIT;
diff --git a/Examples/html/example_EscapeTimes.png b/Examples/html/example_EscapeTimes.png
index b6b4153..e591f61 100644
Binary files a/Examples/html/example_EscapeTimes.png and b/Examples/html/example_EscapeTimes.png differ
diff --git a/Examples/html/example_EscapeTimes_01.png b/Examples/html/example_EscapeTimes_01.png
index 6b5657c..d5420f7 100644
Binary files a/Examples/html/example_EscapeTimes_01.png and b/Examples/html/example_EscapeTimes_01.png differ
diff --git a/Examples/html/example_EscapeTimes_02.png b/Examples/html/example_EscapeTimes_02.png
index b4416fe..2ed9015 100644
Binary files a/Examples/html/example_EscapeTimes_02.png and b/Examples/html/example_EscapeTimes_02.png differ
diff --git a/Examples/html/example_EscapeTimes_03.png b/Examples/html/example_EscapeTimes_03.png
index 949879d..ae21492 100644
Binary files a/Examples/html/example_EscapeTimes_03.png and b/Examples/html/example_EscapeTimes_03.png differ
diff --git a/Examples/html/example_EscapeTimes_04.png b/Examples/html/example_EscapeTimes_04.png
index f4a0748..099dd71 100644
Binary files a/Examples/html/example_EscapeTimes_04.png and b/Examples/html/example_EscapeTimes_04.png differ
diff --git a/Examples/html/example_ExtrinsicNoiseModels.html b/Examples/html/example_ExtrinsicNoiseModels.html
index faea29f..ceae636 100644
--- a/Examples/html/example_ExtrinsicNoiseModels.html
+++ b/Examples/html/example_ExtrinsicNoiseModels.html
@@ -6,7 +6,7 @@
    example_ExtrinsicNoiseModels

Contents

example_ExtrinsicNoiseModels

In this example, we show how to sample an FSM model over intrinsic noise in its various parameters.

close all
 clear all
-addpath('../CommandLine');
+addpath(genpath('../src'));
 

Example 1 - transcription and translation

First create a full model (e.g., for mRNA and protein)

Model1 = SSIT();
 Model1.species = {'rna','protein'};
 Model1.initialCondition = [0;0];
@@ -80,6 +80,7 @@
 Model1.fspOptions.initApproxSS = false;
 Model1.tSpan = linspace(0,5,10);
 [fspSoln1,Model1.fspOptions.bounds] = Model1.solve;
+Model1 = Model1.formPropensitiesGeneral('Model1');
 Model1.makePlot(fspSoln1,'marginals',[],[],[2,3])
 

Generate, solve and plot results for an extrinsic noise version

Specify the rules for the extrinsic noise. This must be a function that returns a vector of parameters, that are in the same order as the parameters provided. You can choose a different distribution for the extrinsic noise in each parameter.

parDistributions = @()[10+2*randn,...
     0.5+0.1*randn,...
@@ -87,6 +88,7 @@
     1+0.2*randn];
 extrinsicModel = extrinsicSSIT(Model1,parDistributions,20);
 Model1.makePlot(extrinsicModel.averagedResults,'marginals',[],[],[2,3])
+legend('Intrinsic Only','Intrinsic + Extrinsic','Location','EastOutside')
 
\ No newline at end of file diff --git a/Examples/html/example_ExtrinsicNoiseModels.png b/Examples/html/example_ExtrinsicNoiseModels.png index 3bf93a6..d94aaf6 100644 Binary files a/Examples/html/example_ExtrinsicNoiseModels.png and b/Examples/html/example_ExtrinsicNoiseModels.png differ diff --git a/Examples/html/example_ExtrinsicNoiseModels_01.png b/Examples/html/example_ExtrinsicNoiseModels_01.png index c8a7c6d..7c4efb0 100644 Binary files a/Examples/html/example_ExtrinsicNoiseModels_01.png and b/Examples/html/example_ExtrinsicNoiseModels_01.png differ diff --git a/Examples/html/example_ExtrinsicNoiseModels_02.png b/Examples/html/example_ExtrinsicNoiseModels_02.png index 2dc5fde..429803d 100644 Binary files a/Examples/html/example_ExtrinsicNoiseModels_02.png and b/Examples/html/example_ExtrinsicNoiseModels_02.png differ diff --git a/Examples/html/example_ExtrinsicNoiseModels_03.png b/Examples/html/example_ExtrinsicNoiseModels_03.png index f0249a6..100ce8c 100644 Binary files a/Examples/html/example_ExtrinsicNoiseModels_03.png and b/Examples/html/example_ExtrinsicNoiseModels_03.png differ diff --git a/Examples/html/example_ExtrinsicNoiseModels_04.png b/Examples/html/example_ExtrinsicNoiseModels_04.png index 0f7f53f..704b465 100644 Binary files a/Examples/html/example_ExtrinsicNoiseModels_04.png and b/Examples/html/example_ExtrinsicNoiseModels_04.png differ diff --git a/Examples/html/example_FIMCalculation.html b/Examples/html/example_FIMCalculation.html index b4ebc4a..a77ed26 100644 --- a/Examples/html/example_FIMCalculation.html +++ b/Examples/html/example_FIMCalculation.html @@ -6,7 +6,7 @@ example_FIMCalculation

Contents

fimExample

In this script, we show how to set up and solve the FSP-FIM matrix with partial observations and probabilistic distortion.

clear all
 close all
+addpath(genpath('../src'));
 

(1) Set up Model

ModelChoice = 'CentralDogma';  % Two species problem (mRNa and protein)
 F2 = SSIT(ModelChoice);
 F2 = F2.formPropensitiesGeneral('FIMExample');
 

(2) Solve FSP for model

F2.solutionScheme = 'FSP';    % Set solution scheme to FSP.
 [FSPsoln,F2.fspOptions.bounds] = F2.solve;  % Solve the FSP analysis
-

(3) Solve FSP Sensitivity

F2.solutionScheme = 'fspSens'; % Set solutions scheme to FSP Sensitivity
+% Plot the results from the FSP analysis
+fig1 = figure(1);clf; set(fig1,'Name','Marginal Distributions, mRNA');
+fig2 = figure(2);clf; set(fig2,'Name','Marginal Distributions, Protein');
+F2.makePlot(FSPsoln,'marginals',[1:4:21],false,[fig1,fig2],{'r-','linewidth',2})
+

(3) Solve FSP Sensitivity

F2.solutionScheme = 'fspSens'; % Set solutions scheme to FSP Sensitivity
 [sensSoln,bounds] = F2.solve(FSPsoln.stateSpace);  % Solve the sensitivity problem
-
Error with Analytical Sensitivity Calculations - Switching to Finite Difference Method
-

(4) Compute FIM using FSP Sensitivity Results

fimResults = F2.computeFIM(sensSoln.sens); % Compute the FIM for full observations and no distortion.
+% Plot the results from the sensitivity analysis
+fig3 = figure(3);clf; set(fig3,'Name','Marginal Sensitivity, mRNA');
+fig4 = figure(4);clf; set(fig4,'Name','Marginal Sensitivity, Protein');
+F2.makePlot(sensSoln,'marginals',[],false,[fig3,fig4],{'b','linewidth',2})
+

(4) Compute FIM using FSP Sensitivity Results

fimResults = F2.computeFIM(sensSoln.sens); % Compute the FIM for full observations and no distortion.
 cellCounts = 10*ones(size(F2.tSpan));  % Number of cells in each experiment.
 [fimTotal,mleCovEstimate,fimMetrics] = F2.evaluateExperiment(fimResults,cellCounts)
-
fimTotal =
+fig5 = figure(5);clf; set(fig5,'Name','Fim-Predicted Uncertainty Ellipses');
+F2.plotMHResults([],fimTotal,'lin',[],fig5)
+legend('FIM - Full Observation')
+
+fimTotal =
+
   1×1 cell array
+
     {4×4 double}
+
+
 mleCovEstimate =
+
   1×1 cell array
+
     {4×4 double}
+
+
 fimMetrics = 
+
   struct with fields:
 
-          det: 15056909185.6849
-        trace: 31032.4672821342
-    minEigVal: 1.01963637716716
-

(5) Compute FIM for Partial Observations

F2.pdoOptions.PDO=[];
+          det: 15056904911.6332
+        trace: 31032.4672181693
+    minEigVal: 1.01963610988029
+
+

(5) Compute FIM for Partial Observations

F2.pdoOptions.PDO=[];
 F2.pdoOptions.unobservedSpecies = 'x1';
 [fimResults_partialObs] = F2.computeFIM(sensSoln.sens); % Compute the FIM for full observations and no distortion.
 [fimTotal_partialObs,mleCovEstimate_partialObs,fimMetrics_partialObs] = F2.evaluateExperiment(fimResults_partialObs,cellCounts)
-
fimTotal_partialObs =
+fig6 = figure(6);clf; set(fig6,'Name','Fim-Predicted Uncertainty Ellipses');
+F2.plotMHResults([],[fimTotal,fimTotal_partialObs],'lin',[],fig6)
+legend('FIM - Full Observation','FIM - Protein Only')
+
+fimTotal_partialObs =
+
   1×1 cell array
+
     {4×4 double}
+
+
 mleCovEstimate_partialObs =
+
   1×1 cell array
+
     {4×4 double}
+
+
 fimMetrics_partialObs = 
+
   struct with fields:
 
-          det: 15334575.7734958
-        trace: 28026.2035489583
-    minEigVal: 0.265831411968639
-

(6) Compute FIM for Distorted Observation (Probabilistic Distortion Operator)

F2.pdoOptions.unobservedSpecies = 'x1';
+          det: 15334574.1344478
+        trace: 28026.2035042016
+    minEigVal: 0.265831373533071
+
+

(6) Compute FIM for Distorted Observation (Probabilistic Distortion Operator)

F2.pdoOptions.unobservedSpecies = 'x1';
 pdoOptions.type = 'Binomial';
 % Need to define loss parameter for each species S1, S2,...
 pdoOptions.props.CaptureProbabilityS1 = 0;  % Use zero for unobserved species.
 pdoOptions.props.CaptureProbabilityS2 = 0.9;
 
-% call method to generate the PDO.
+% Call method to generate the PDO.
 F2.pdoOptions.PDO = F2.generatePDO(pdoOptions,[],FSPsoln.fsp);
 
-% plot the PDO
+% Plot the PDO
 N = size(F2.pdoOptions.PDO.conditionalPmfs{1});
-figure;
+fig7 = figure(7); set(fig7,'Name','Probabilistic Distortion Operator for Protein');
 contourf([0:N(1)-1],[0:N(2)-1],log10(F2.pdoOptions.PDO.conditionalPmfs{1}));
 % Here we wanted the first PDO for 'x2' because 'x1' was unobserved.
 xlabel('Actual');ylabel('Observable');colorbar;
 
-% solve FIM using the specified PDO
+% Solve FIM using the specified PDO
 [fimResults_BinomialPDO] = F2.computeFIM(sensSoln.sens); % Compute the FIM for full observations and no distortion.
 [fimTotal_BinomialPDO,mleCovEstimate_BinomialPDO,fimMetrics_BinomialPDO] =...
     F2.evaluateExperiment(fimResults_BinomialPDO,cellCounts)
-
fimTotal_BinomialPDO =
+fig8 = figure(8);clf; set(fig8,'Name','Fim-Predicted Uncertainty Ellipses');
+F2.plotMHResults([],[fimTotal,fimTotal_partialObs,fimTotal_BinomialPDO],'lin',[],fig8)
+legend('FIM - Full Observation','FIM - Protein Only','FIM - Protein with Error')
+
+fimTotal_BinomialPDO =
+
   1×1 cell array
+
     {4×4 double}
+
+
 mleCovEstimate_BinomialPDO =
+
   1×1 cell array
+
     {4×4 double}
+
+
 fimMetrics_BinomialPDO = 
+
   struct with fields:
 
-          det: 11753388.9719856
-        trace: 26697.1703866449
-    minEigVal: 0.240145247197538
-
\ No newline at end of file diff --git a/Examples/html/example_FIMCalculation.png b/Examples/html/example_FIMCalculation.png index 91240b2..bffcd66 100644 Binary files a/Examples/html/example_FIMCalculation.png and b/Examples/html/example_FIMCalculation.png differ diff --git a/Examples/html/example_FIMCalculation_01.png b/Examples/html/example_FIMCalculation_01.png index 711bce6..c7f97c8 100644 Binary files a/Examples/html/example_FIMCalculation_01.png and b/Examples/html/example_FIMCalculation_01.png differ diff --git a/Examples/html/example_FIMCalculation_02.png b/Examples/html/example_FIMCalculation_02.png new file mode 100644 index 0000000..45877ea Binary files /dev/null and b/Examples/html/example_FIMCalculation_02.png differ diff --git a/Examples/html/example_FIMCalculation_03.png b/Examples/html/example_FIMCalculation_03.png new file mode 100644 index 0000000..73465c2 Binary files /dev/null and b/Examples/html/example_FIMCalculation_03.png differ diff --git a/Examples/html/example_FIMCalculation_04.png b/Examples/html/example_FIMCalculation_04.png new file mode 100644 index 0000000..891edb4 Binary files /dev/null and b/Examples/html/example_FIMCalculation_04.png differ diff --git a/Examples/html/example_FIMCalculation_05.png b/Examples/html/example_FIMCalculation_05.png new file mode 100644 index 0000000..1042345 Binary files /dev/null and b/Examples/html/example_FIMCalculation_05.png differ diff --git a/Examples/html/example_FIMCalculation_06.png b/Examples/html/example_FIMCalculation_06.png new file mode 100644 index 0000000..29d9a96 Binary files /dev/null and b/Examples/html/example_FIMCalculation_06.png differ diff --git a/Examples/html/example_FIMCalculation_07.png b/Examples/html/example_FIMCalculation_07.png new file mode 100644 index 0000000..424ecae Binary files /dev/null and b/Examples/html/example_FIMCalculation_07.png differ diff --git a/Examples/html/example_FIMCalculation_08.png b/Examples/html/example_FIMCalculation_08.png new file mode 100644 index 0000000..7b19881 Binary files /dev/null and b/Examples/html/example_FIMCalculation_08.png differ diff --git a/Examples/html/example_Fit_STL1_mRNA_Data.html b/Examples/html/example_Fit_STL1_mRNA_Data.html index fa3c6ae..b4e3323 100644 --- a/Examples/html/example_Fit_STL1_mRNA_Data.html +++ b/Examples/html/example_Fit_STL1_mRNA_Data.html @@ -6,7 +6,7 @@ example_Fit_STL1_mRNA_Data

Contents

Example script to show SSIT for the Vanderbilt Q_BIO Group

In this script, we are going to show how to create, solve and fit a CME model to some single-cell smFISH data. For this example, we will use some data collected in Dr. Gregor Neuert's laboratory at Vanderbilt.

close all
 clear all
-addpath('../CommandLine')
+addpath(genpath('../src'));
 

Create SSIT Model

First, we are going to create an FSP model for a bursting gene expression model. This model will consist of 3 species: OFF Gene, ON Gene, and mRNA. There will be four reactions: activation, inactivation, transcription and degradation. The activation rate will be assumed to be time varying and controlled by a MAPK signal.

Model = SSIT;    % Create SSIT instance and call it 'Model'.
 
 % Set species names for bursting gene expression model:
@@ -110,7 +110,7 @@
 
 % Set the code to start at steady state at t=0;
 Model.fspOptions.initApproxSS =true;
-Model = Model.formPropensitiesGeneral('STL1Model');
+Model = Model.formPropensitiesGeneral('STL1Model',true);
 [FSPsoln,Model.fspOptions.bounds] = Model.solve;  % Solve the FSP analysis
 
 % Next we make plots of the marginal distributions at time points 3, 5, 7,
@@ -260,37 +260,37 @@
 % again - it is possible you can still find a better model to explain your
 % data.  This can take several rounds of iteration before convergence.  I
 % recommend creating a while loop to make it automated.
-
n=0; acc=0.08. TMPmh_24.mat
-n=100; acc=0.125. TMPmh_24.mat
-n=200; acc=0.093333. TMPmh_24.mat
-n=300; acc=0.09. TMPmh_24.mat
-n=400; acc=0.084. TMPmh_24.mat
-n=500; acc=0.071667. TMPmh_24.mat
-n=600; acc=0.077143. TMPmh_24.mat
-n=700; acc=0.075. TMPmh_24.mat
-n=800; acc=0.068889. TMPmh_24.mat
-n=900; acc=0.067. TMPmh_24.mat
-n=1000; acc=0.063636. TMPmh_24.mat
-n=1100; acc=0.066667. TMPmh_24.mat
-n=1200; acc=0.065385. TMPmh_24.mat
-n=1300; acc=0.065. TMPmh_24.mat
-n=1400; acc=0.064667. TMPmh_24.mat
-n=1500; acc=0.066875. TMPmh_24.mat
-n=1600; acc=0.069412. TMPmh_24.mat
-n=1700; acc=0.070556. TMPmh_24.mat
-n=1800; acc=0.071053. TMPmh_24.mat
-n=1900; acc=0.0715. TMPmh_24.mat
-n=2000; acc=0.071905. TMPmh_24.mat
-n=2100; acc=0.075455. TMPmh_24.mat
-n=2200; acc=0.075217. TMPmh_24.mat
-n=2300; acc=0.077083. TMPmh_24.mat
-n=2400; acc=0.0768. TMPmh_24.mat
-n=2500; acc=0.076154. TMPmh_24.mat
-n=2600; acc=0.074074. TMPmh_24.mat
-n=2700; acc=0.073571. TMPmh_24.mat
-n=2800; acc=0.072759. TMPmh_24.mat
-n=2900; acc=0.072. TMPmh_24.mat
-n=3000; acc=0.07129. TMPmh_24.mat
+
n=0; acc=0.13. TMPmh_2.mat
+n=100; acc=0.135. TMPmh_2.mat
+n=200; acc=0.15. TMPmh_2.mat
+n=300; acc=0.1375. TMPmh_2.mat
+n=400; acc=0.128. TMPmh_2.mat
+n=500; acc=0.12333. TMPmh_2.mat
+n=600; acc=0.12286. TMPmh_2.mat
+n=700; acc=0.12625. TMPmh_2.mat
+n=800; acc=0.12556. TMPmh_2.mat
+n=900; acc=0.122. TMPmh_2.mat
+n=1000; acc=0.12. TMPmh_2.mat
+n=1100; acc=0.12. TMPmh_2.mat
+n=1200; acc=0.11615. TMPmh_2.mat
+n=1300; acc=0.11857. TMPmh_2.mat
+n=1400; acc=0.11467. TMPmh_2.mat
+n=1500; acc=0.10938. TMPmh_2.mat
+n=1600; acc=0.10941. TMPmh_2.mat
+n=1700; acc=0.10722. TMPmh_2.mat
+n=1800; acc=0.10474. TMPmh_2.mat
+n=1900; acc=0.1065. TMPmh_2.mat
+n=2000; acc=0.10476. TMPmh_2.mat
+n=2100; acc=0.10545. TMPmh_2.mat
+n=2200; acc=0.10435. TMPmh_2.mat
+n=2300; acc=0.10208. TMPmh_2.mat
+n=2400; acc=0.1028. TMPmh_2.mat
+n=2500; acc=0.10269. TMPmh_2.mat
+n=2600; acc=0.10111. TMPmh_2.mat
+n=2700; acc=0.10107. TMPmh_2.mat
+n=2800; acc=0.098621. TMPmh_2.mat
+n=2900; acc=0.096667. TMPmh_2.mat
+n=3000; acc=0.095161. TMPmh_2.mat
 

Iterating between MLE and MH.

Let's run a few rounds of MLE and MH to see if we can get better convergence.

Model.parameters(:,2) = num2cell(pars);
 for i=1:3
     % Maximize likelihood
@@ -319,109 +319,110 @@
  the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 
  and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04 
 
-n=0; acc=0.31. TMPmh_25.mat
-n=100; acc=0.395. TMPmh_25.mat
-n=200; acc=0.39. TMPmh_25.mat
-n=300; acc=0.39. TMPmh_25.mat
-n=400; acc=0.4. TMPmh_25.mat
-n=500; acc=0.42. TMPmh_25.mat
-n=600; acc=0.42571. TMPmh_25.mat
-n=700; acc=0.43375. TMPmh_25.mat
-n=800; acc=0.43333. TMPmh_25.mat
-n=900; acc=0.435. TMPmh_25.mat
-n=1000; acc=0.43091. TMPmh_25.mat
-n=1100; acc=0.42917. TMPmh_25.mat
-n=1200; acc=0.43308. TMPmh_25.mat
-n=1300; acc=0.43286. TMPmh_25.mat
-n=1400; acc=0.426. TMPmh_25.mat
-n=1500; acc=0.42438. TMPmh_25.mat
-n=1600; acc=0.42941. TMPmh_25.mat
-n=1700; acc=0.43111. TMPmh_25.mat
-n=1800; acc=0.42895. TMPmh_25.mat
-n=1900; acc=0.427. TMPmh_25.mat
-n=2000; acc=0.42762. TMPmh_25.mat
-n=2100; acc=0.42545. TMPmh_25.mat
-n=2200; acc=0.4287. TMPmh_25.mat
-n=2300; acc=0.42792. TMPmh_25.mat
-n=2400; acc=0.4312. TMPmh_25.mat
-n=2500; acc=0.43077. TMPmh_25.mat
-n=2600; acc=0.42852. TMPmh_25.mat
-n=2700; acc=0.42929. TMPmh_25.mat
-n=2800; acc=0.43138. TMPmh_25.mat
-n=2900; acc=0.43333. TMPmh_25.mat
-n=3000; acc=0.43484. TMPmh_25.mat
+n=0; acc=0.42. TMPmh_3.mat
+n=100; acc=0.385. TMPmh_3.mat
+n=200; acc=0.42667. TMPmh_3.mat
+n=300; acc=0.425. TMPmh_3.mat
+n=400; acc=0.422. TMPmh_3.mat
+n=500; acc=0.43167. TMPmh_3.mat
+n=600; acc=0.43571. TMPmh_3.mat
+n=700; acc=0.43125. TMPmh_3.mat
+n=800; acc=0.42778. TMPmh_3.mat
+n=900; acc=0.422. TMPmh_3.mat
+n=1000; acc=0.42545. TMPmh_3.mat
+n=1100; acc=0.42083. TMPmh_3.mat
+n=1200; acc=0.41692. TMPmh_3.mat
+n=1300; acc=0.41714. TMPmh_3.mat
+n=1400; acc=0.41867. TMPmh_3.mat
+n=1500; acc=0.41625. TMPmh_3.mat
+n=1600; acc=0.41647. TMPmh_3.mat
+n=1700; acc=0.41889. TMPmh_3.mat
+n=1800; acc=0.41789. TMPmh_3.mat
+n=1900; acc=0.4205. TMPmh_3.mat
+n=2000; acc=0.42143. TMPmh_3.mat
+n=2100; acc=0.41591. TMPmh_3.mat
+n=2200; acc=0.4187. TMPmh_3.mat
+n=2300; acc=0.42208. TMPmh_3.mat
+n=2400; acc=0.4216. TMPmh_3.mat
+n=2500; acc=0.42192. TMPmh_3.mat
+n=2600; acc=0.42222. TMPmh_3.mat
+n=2700; acc=0.41929. TMPmh_3.mat
+n=2800; acc=0.42241. TMPmh_3.mat
+n=2900; acc=0.42533. TMPmh_3.mat
+n=3000; acc=0.42742. TMPmh_3.mat
  
 Optimization terminated:
  the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 
  and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04 
 
-n=0; acc=0.33. TMPmh_26.mat
-n=100; acc=0.375. TMPmh_26.mat
-n=200; acc=0.39667. TMPmh_26.mat
-n=300; acc=0.4. TMPmh_26.mat
-n=400; acc=0.408. TMPmh_26.mat
-n=500; acc=0.41. TMPmh_26.mat
-n=600; acc=0.40429. TMPmh_26.mat
-n=700; acc=0.415. TMPmh_26.mat
-n=800; acc=0.41889. TMPmh_26.mat
-n=900; acc=0.425. TMPmh_26.mat
-n=1000; acc=0.42273. TMPmh_26.mat
-n=1100; acc=0.4175. TMPmh_26.mat
-n=1200; acc=0.41385. TMPmh_26.mat
-n=1300; acc=0.41571. TMPmh_26.mat
-n=1400; acc=0.41933. TMPmh_26.mat
-n=1500; acc=0.42562. TMPmh_26.mat
-n=1600; acc=0.42235. TMPmh_26.mat
-n=1700; acc=0.42111. TMPmh_26.mat
-n=1800; acc=0.42211. TMPmh_26.mat
-n=1900; acc=0.424. TMPmh_26.mat
-n=2000; acc=0.42286. TMPmh_26.mat
-n=2100; acc=0.42182. TMPmh_26.mat
-n=2200; acc=0.42304. TMPmh_26.mat
-n=2300; acc=0.42125. TMPmh_26.mat
-n=2400; acc=0.4208. TMPmh_26.mat
-n=2500; acc=0.42038. TMPmh_26.mat
-n=2600; acc=0.42333. TMPmh_26.mat
-n=2700; acc=0.4225. TMPmh_26.mat
-n=2800; acc=0.42345. TMPmh_26.mat
-n=2900; acc=0.421. TMPmh_26.mat
-n=3000; acc=0.41839. TMPmh_26.mat
+n=0; acc=0.44. TMPmh_4.mat
+n=100; acc=0.425. TMPmh_4.mat
+n=200; acc=0.4. TMPmh_4.mat
+n=300; acc=0.415. TMPmh_4.mat
+n=400; acc=0.418. TMPmh_4.mat
+n=500; acc=0.42833. TMPmh_4.mat
+n=600; acc=0.41143. TMPmh_4.mat
+n=700; acc=0.39875. TMPmh_4.mat
+n=800; acc=0.39333. TMPmh_4.mat
+n=900; acc=0.4. TMPmh_4.mat
+n=1000; acc=0.41091. TMPmh_4.mat
+n=1100; acc=0.4125. TMPmh_4.mat
+n=1200; acc=0.41769. TMPmh_4.mat
+n=1300; acc=0.42143. TMPmh_4.mat
+n=1400; acc=0.42267. TMPmh_4.mat
+n=1500; acc=0.41937. TMPmh_4.mat
+n=1600; acc=0.42647. TMPmh_4.mat
+n=1700; acc=0.42667. TMPmh_4.mat
+n=1800; acc=0.42632. TMPmh_4.mat
+n=1900; acc=0.428. TMPmh_4.mat
+n=2000; acc=0.42857. TMPmh_4.mat
+n=2100; acc=0.42545. TMPmh_4.mat
+n=2200; acc=0.42522. TMPmh_4.mat
+n=2300; acc=0.42. TMPmh_4.mat
+n=2400; acc=0.4244. TMPmh_4.mat
+n=2500; acc=0.42385. TMPmh_4.mat
+n=2600; acc=0.42259. TMPmh_4.mat
+n=2700; acc=0.4225. TMPmh_4.mat
+n=2800; acc=0.42345. TMPmh_4.mat
+n=2900; acc=0.42533. TMPmh_4.mat
+n=3000; acc=0.42613. TMPmh_4.mat
  
 Optimization terminated:
  the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 
  and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04 
 
-n=0; acc=0.44. TMPmh_27.mat
-n=100; acc=0.415. TMPmh_27.mat
-n=200; acc=0.43667. TMPmh_27.mat
-n=300; acc=0.4075. TMPmh_27.mat
-n=400; acc=0.404. TMPmh_27.mat
-n=500; acc=0.40167. TMPmh_27.mat
-n=600; acc=0.40286. TMPmh_27.mat
-n=700; acc=0.38875. TMPmh_27.mat
-n=800; acc=0.39222. TMPmh_27.mat
-n=900; acc=0.398. TMPmh_27.mat
-n=1000; acc=0.39455. TMPmh_27.mat
-n=1100; acc=0.39417. TMPmh_27.mat
-n=1200; acc=0.39385. TMPmh_27.mat
-n=1300; acc=0.395. TMPmh_27.mat
-n=1400; acc=0.396. TMPmh_27.mat
-n=1500; acc=0.39937. TMPmh_27.mat
-n=1600; acc=0.39412. TMPmh_27.mat
-n=1700; acc=0.39778. TMPmh_27.mat
-n=1800; acc=0.40105. TMPmh_27.mat
-n=1900; acc=0.403. TMPmh_27.mat
-n=2000; acc=0.41048. TMPmh_27.mat
-n=2100; acc=0.41045. TMPmh_27.mat
-n=2200; acc=0.40435. TMPmh_27.mat
-n=2300; acc=0.40542. TMPmh_27.mat
-n=2400; acc=0.404. TMPmh_27.mat
-n=2500; acc=0.40077. TMPmh_27.mat
-n=2600; acc=0.4037. TMPmh_27.mat
-n=2700; acc=0.4075. TMPmh_27.mat
-n=2800; acc=0.40862. TMPmh_27.mat
-n=2900; acc=0.40967. TMPmh_27.mat
-n=3000; acc=0.41323. TMPmh_27.mat
+Parallel pool using the 'Processes' profile is shutting down.
+n=0; acc=0.42. TMPmh_5.mat
+n=100; acc=0.41. TMPmh_5.mat
+n=200; acc=0.44333. TMPmh_5.mat
+n=300; acc=0.4475. TMPmh_5.mat
+n=400; acc=0.452. TMPmh_5.mat
+n=500; acc=0.44. TMPmh_5.mat
+n=600; acc=0.44143. TMPmh_5.mat
+n=700; acc=0.4475. TMPmh_5.mat
+n=800; acc=0.44111. TMPmh_5.mat
+n=900; acc=0.439. TMPmh_5.mat
+n=1000; acc=0.43818. TMPmh_5.mat
+n=1100; acc=0.44. TMPmh_5.mat
+n=1200; acc=0.43923. TMPmh_5.mat
+n=1300; acc=0.435. TMPmh_5.mat
+n=1400; acc=0.43333. TMPmh_5.mat
+n=1500; acc=0.43875. TMPmh_5.mat
+n=1600; acc=0.43529. TMPmh_5.mat
+n=1700; acc=0.43167. TMPmh_5.mat
+n=1800; acc=0.43421. TMPmh_5.mat
+n=1900; acc=0.4395. TMPmh_5.mat
+n=2000; acc=0.43429. TMPmh_5.mat
+n=2100; acc=0.43409. TMPmh_5.mat
+n=2200; acc=0.43304. TMPmh_5.mat
+n=2300; acc=0.43042. TMPmh_5.mat
+n=2400; acc=0.432. TMPmh_5.mat
+n=2500; acc=0.43038. TMPmh_5.mat
+n=2600; acc=0.4337. TMPmh_5.mat
+n=2700; acc=0.43393. TMPmh_5.mat
+n=2800; acc=0.43621. TMPmh_5.mat
+n=2900; acc=0.43467. TMPmh_5.mat
+n=3000; acc=0.43484. TMPmh_5.mat
 

Evaluating the MH results

Here we will generate three plots. The first one will show the Likelihood function as we search over parameter space. In order to get a good estimate of the parameter uncertainty, we want this to quickly reach ts maximum value and then to fluctuate around that value for a significant amount of time. If you see that it is still incereasing, you know that the fit has not yet converged.

figure
 subplot(1,3,1)
 plot(chainResults.mhValue)
@@ -469,7 +470,7 @@
 
 Neff =
 
-           41.266228308708
+          43.1954111381028
 
 
example_FittingAndDesigningExperiments

Contents

example_FittingAndDesigningExperiments

In this script, we show how the SSIT can be used to identify a time-inhomogeneous model for the activation of Dusp1 mRNA expression under Dexamethasome stimulation of Glucocorticoid Receptors.

clear all
 clc
 close all
-addpath('../CommandLine');
+addpath(genpath('../src'));
 

Define SSIT Model

Here we set up a simple model where there is an upstream transcription factor (GR) that activates a gene. Once active, the gene can transcribe nuclear RNA, which can later decay or leave the nucleus.

Model = SSIT;
 Model.species = {'activeGene';'rna'};
 Model.initialCondition = [0;0];
@@ -83,6 +83,7 @@
 

Solve the model using the FSP

Next, we can solve the model using the FSP. In this example, we show how to run the code twice. First call finds the FSP projection needed to solve the problem, and the second call solves using that projection.

Model.solutionScheme = 'FSP';
 Model.fspOptions.fspTol = 1e-4;
 Model.fspOptions.bounds(3:4) = [2,400];
+Model = Model.formPropensitiesGeneral('Model');
 [fspSoln,Model.fspOptions.bounds] = Model.solve;
 

Load and Fit smFISH Data

Next we load experimental data from a CSV file and associate the species 'rna' with the collumn 'RNA_nuc'.

Model = Model.loadData('../ExampleData/DUSP1_Dex_100nM_Rep1_Rep2.csv',{'rna','RNA_nuc'});
 Model.fittingOptions.modelVarsToFit = 1:7;
@@ -106,7 +107,7 @@
 Model.fspOptions.fspTol = 1e-6;
 [sensSoln] = Model.solve(fspSoln.stateSpace);
 Model.makePlot(sensSoln,'marginals')
-

Compute Fisher Information (in log-parameter space)

Model.pdoOptions.unobservedSpecies = {'activeGene'};
+

Compute Fisher Information (in log-parameter space)

Model.pdoOptions.unobservedSpecies = {'activeGene'};
 fimResults = Model.computeFIM(sensSoln.sens,'log');
 FIMlog =  Model.evaluateExperiment(fimResults,Model.dataSet.nCells);
 

Metropolis Hastings to Quantify Parameter Uncertainty

Next, we determine the uncertainty in a subset of the estimated parameters. In this case, we assume that the GR parameters ([5:7]) are correct, and we only care about the uncertainity in the other parameters ([1:4]).

Model.solutionScheme = 'FSP';
@@ -122,54 +123,54 @@
 % calculated using the first parameters set, this will cause the MH to
 % converge more slowly.  You will likely want to re-run the fit again
 % multiple times to have a better chance of convergence.
-
n=0; acc=0.5. TMPmh_5.mat
-n=100; acc=0.485. TMPmh_5.mat
-n=200; acc=0.47667. TMPmh_5.mat
-n=300; acc=0.495. TMPmh_5.mat
-n=400; acc=0.504. TMPmh_5.mat
-n=500; acc=0.50167. TMPmh_5.mat
-n=600; acc=0.49143. TMPmh_5.mat
-n=700; acc=0.485. TMPmh_5.mat
-n=800; acc=0.48444. TMPmh_5.mat
-n=900; acc=0.48. TMPmh_5.mat
-n=1000; acc=0.48273. TMPmh_5.mat
-n=1100; acc=0.48083. TMPmh_5.mat
-n=1200; acc=0.48. TMPmh_5.mat
-n=1300; acc=0.47643. TMPmh_5.mat
-n=1400; acc=0.47467. TMPmh_5.mat
-n=1500; acc=0.475. TMPmh_5.mat
-n=1600; acc=0.47353. TMPmh_5.mat
-n=1700; acc=0.47611. TMPmh_5.mat
-n=1800; acc=0.47053. TMPmh_5.mat
-n=1900; acc=0.474. TMPmh_5.mat
-n=2000; acc=0.48. TMPmh_5.mat
-

Solve for Fisher Information Matrix at all Time Points

Using the sensitivity calculation, we next computer the Fisher Information Matrix. Because we only can observe one species 'rna'), we must include that the species 'activeGene' is not observed when formulating the PDO.

Model.pdoOptions.unobservedSpecies = {'activeGene'};
+
n=0; acc=0.51. TMPmh_1.mat
+n=100; acc=0.52. TMPmh_1.mat
+n=200; acc=0.5. TMPmh_1.mat
+n=300; acc=0.4775. TMPmh_1.mat
+n=400; acc=0.474. TMPmh_1.mat
+n=500; acc=0.47833. TMPmh_1.mat
+n=600; acc=0.48571. TMPmh_1.mat
+n=700; acc=0.4925. TMPmh_1.mat
+n=800; acc=0.5. TMPmh_1.mat
+n=900; acc=0.496. TMPmh_1.mat
+n=1000; acc=0.48818. TMPmh_1.mat
+n=1100; acc=0.49417. TMPmh_1.mat
+n=1200; acc=0.49615. TMPmh_1.mat
+n=1300; acc=0.49143. TMPmh_1.mat
+n=1400; acc=0.48933. TMPmh_1.mat
+n=1500; acc=0.48875. TMPmh_1.mat
+n=1600; acc=0.48882. TMPmh_1.mat
+n=1700; acc=0.48667. TMPmh_1.mat
+n=1800; acc=0.48316. TMPmh_1.mat
+n=1900; acc=0.4805. TMPmh_1.mat
+n=2000; acc=0.48476. TMPmh_1.mat
+

Solve for Fisher Information Matrix at all Time Points

Using the sensitivity calculation, we next computer the Fisher Information Matrix. Because we only can observe one species 'rna'), we must include that the species 'activeGene' is not observed when formulating the PDO.

Model.pdoOptions.unobservedSpecies = {'activeGene'};
 fims = Model.computeFIM(sensSoln.sens);
 FIM = Model.evaluateExperiment(fims,Model.dataSet.nCells);
 Model.plotMHResults(mhResults,FIM);
-

Optimize Experiment Design (Same Number of Cells and Timepoints)

Now that we have calculate the FIM and PDO, we can optimize an experiment design to match the total number of cellular measurements as conducted in the original experiment.

nTotal = sum(Model.dataSet.nCells);
+

Optimize Experiment Design (Same Number of Cells and Timepoints)

Now that we have calculate the FIM and PDO, we can optimize an experiment design to match the total number of cellular measurements as conducted in the original experiment.

nTotal = sum(Model.dataSet.nCells);
 nCellsOpt = Model.optimizeCellCounts(fims,nTotal,'TR[1:4]');
 fimOpt = Model.evaluateExperiment(fims,nCellsOpt);
 Model.plotMHResults(mhResults,[FIM,fimOpt]);
-

Calibrate PDO from Multi-Modal Experimental Data

Calibration the PDO from empirical data. Here, the number of spots has been measured using different assays in data columns 'nTotal' for the 'true' data set and in the columns 'nSpots0' for a different label or 'intens1' for the integrated intensity. We calibrate two different PDOs for this case. In both cases, we assume an 'AffinePoiss' PDO where the obervation probability is a Poisson distribution where the mean value is affine linearly related to the true value: P(y|x) = Poiss(a0 + a1*x);

ModelPDOSpots = Model.calibratePDO('../ExampleData/pdoCalibrationData.csv',...
+

Calibrate PDO from Multi-Modal Experimental Data

Calibration the PDO from empirical data. Here, the number of spots has been measured using different assays in data columns 'nTotal' for the 'true' data set and in the columns 'nSpots0' for a different label or 'intens1' for the integrated intensity. We calibrate two different PDOs for this case. In both cases, we assume an 'AffinePoiss' PDO where the obervation probability is a Poisson distribution where the mean value is affine linearly related to the true value: P(y|x) = Poiss(a0 + a1*x);

ModelPDOSpots = Model.calibratePDO('../ExampleData/pdoCalibrationData.csv',...
     {'rna'},{'nTotal'},{'nSpots0'},'AffinePoiss',true);
 
 ModelPDOIntens = Model.calibratePDO('../ExampleData/pdoCalibrationData.csv',...
     {'rna'},{'nTotal'},{'intens1'},'AffinePoiss',true,[1,1000,1]);
-

Compute the FIM for the calculated PDO

Now that we have the new PDO, we can calculate the FIM for the different distorted observation experiments. In each case, we compare the results to the Metropolis Hastings results in which the data was assumed to be undistorted. First, for the alternate label distortion.

fimsPDOSpot = ModelPDOSpots.computeFIM(sensSoln.sens);
+

Compute the FIM for the calculated PDO

Now that we have the new PDO, we can calculate the FIM for the different distorted observation experiments. In each case, we compare the results to the Metropolis Hastings results in which the data was assumed to be undistorted. First, for the alternate label distortion:

fimsPDOSpot = ModelPDOSpots.computeFIM(sensSoln.sens);
 fimPDOSpots = ModelPDOSpots.evaluateExperiment(fimsPDOSpot,nCellsOpt);
 Model.plotMHResults(mhResults,[FIM,fimOpt,fimPDOSpots]);
 
-% Next, for the intensity integration distortion.
+% Next, for the intensity integration distortion:
 fimsPDOIntens = ModelPDOIntens.computeFIM(sensSoln.sens);
 fimPDOIntens = ModelPDOIntens.evaluateExperiment(fimsPDOIntens,nCellsOpt);
 Model.plotMHResults(mhResults,[FIM,fimOpt,fimPDOSpots,fimPDOIntens]);
 
-% Finally for an extended experiment with a larger number of cells.
+% Finally, for an extended experiment with a larger number of cells;
 fimsPDOIntens = ModelPDOIntens.computeFIM(sensSoln.sens);
 fimPDOIntens2x = ModelPDOIntens.evaluateExperiment(fimsPDOIntens,2.218*nCellsOpt);
 Model.plotMHResults(mhResults,[FIM,fimOpt,fimPDOSpots,fimPDOIntens,fimPDOIntens2x]);
-
example_HybridModel

Contents

example_HybridModel

In this script, we demonstrate how to adjust a model to treat some species (i.e., upstream reactions) using an ODE formulation, while having other species (i.e., downstream species) evolving in a discrete stochastic manner.

close all
 clear all
-addpath('../CommandLine');
+addpath(genpath('../src'));
 

Example 1 - transcription and translation

First create a full model (e.g., for mRNA and protein)

Model1 = SSIT;
 Model1.species = {'rna','protein'};
 Model1.initialCondition = [0;0];
@@ -86,6 +86,7 @@
 Model2.hybridOptions.upstreamODEs = {'rna'};
 [fspSoln2, Model2.fspOptions.bounds] = Model2.solve;
 Model2.makePlot(fspSoln2,'marginals',[],[],3)
+legend('Full','QSSA for (n)','Location','eastoutside')
 

Example 2 - 5-species MAPK induction Model

In this example, we consider a model of MAPK translocation to the nucleus followed by binding to a gene and then transcription activation. We will assume that there are two alleles with one each starting in the active and inactive state.

Model3 = SSIT;
 Model3.species = {'geneInactive','geneActive','mapkCyt','mapkNuc','rna'};
 Model3.initialCondition = [1;1;50;0;0];
@@ -114,14 +115,13 @@
 Model5.hybridOptions.upstreamODEs = {'mapkCyt','mapkNuc'};
 [fspSoln5, Model5.fspOptions.bounds] = Model5.solve;
 Model5.makePlot(fspSoln5,'marginals',[],[],[11,12,15])
-
Warning: Reaction 3 changes both
-ODE and stochastic species.
-Removing effect on upstream
-species. 
-Warning: Reaction 4 changes both
-ODE and stochastic species.
-Removing effect on upstream
-species. 
+legend('Full','QSSA for (gene,MAKP)','QSSA for (MAKP)','Location','eastoutside')
+
Warning: Reaction 4 changes both ODE and stochastic species. Removing effect on upstream species.
+> In <a href="matlab:matlab.internal.language.introspective.errorDocCallback('parallel_function>make_general_channel/channel_general', '/Applications/MATLAB_R2023a.app/toolbox/matlab/lang/parallel_function.m', 837)" style="font-weight:bold">parallel_function>make_general_channel/channel_general</a> (<a href="matlab: opentoline('/Applications/MATLAB_R2023a.app/toolbox/matlab/lang/parallel_function.m',837,0)">line 837</a>)
+In <a href="matlab:matlab.internal.language.introspective.errorDocCallback('parallel.internal.parfor.cppRemoteParallelFunction', '/Applications/MATLAB_R2023a.app/toolbox/parallel/cluster/+parallel/+internal/+parfor/cppRemoteParallelFunction.m', 53)" style="font-weight:bold">parallel.internal.parfor.cppRemoteParallelFunction</a> (<a href="matlab: opentoline('/Applications/MATLAB_R2023a.app/toolbox/parallel/cluster/+parallel/+internal/+parfor/cppRemoteParallelFunction.m',53,0)">line 53</a>)
+Warning: Reaction 3 changes both ODE and stochastic species. Removing effect on upstream species.
+> In <a href="matlab:matlab.internal.language.introspective.errorDocCallback('parallel_function>make_general_channel/channel_general', '/Applications/MATLAB_R2023a.app/toolbox/matlab/lang/parallel_function.m', 837)" style="font-weight:bold">parallel_function>make_general_channel/channel_general</a> (<a href="matlab: opentoline('/Applications/MATLAB_R2023a.app/toolbox/matlab/lang/parallel_function.m',837,0)">line 837</a>)
+In <a href="matlab:matlab.internal.language.introspective.errorDocCallback('parallel.internal.parfor.cppRemoteParallelFunction', '/Applications/MATLAB_R2023a.app/toolbox/parallel/cluster/+parallel/+internal/+parfor/cppRemoteParallelFunction.m', 53)" style="font-weight:bold">parallel.internal.parfor.cppRemoteParallelFunction</a> (<a href="matlab: opentoline('/Applications/MATLAB_R2023a.app/toolbox/parallel/cluster/+parallel/+internal/+parfor/cppRemoteParallelFunction.m',53,0)">line 53</a>)
 
\ No newline at end of file diff --git a/Examples/html/example_HybridModel.png b/Examples/html/example_HybridModel.png index d725aaf..26ec9c0 100644 Binary files a/Examples/html/example_HybridModel.png and b/Examples/html/example_HybridModel.png differ diff --git a/Examples/html/example_HybridModel_01.png b/Examples/html/example_HybridModel_01.png index 138a633..3fa8b17 100644 Binary files a/Examples/html/example_HybridModel_01.png and b/Examples/html/example_HybridModel_01.png differ diff --git a/Examples/html/example_HybridModel_02.png b/Examples/html/example_HybridModel_02.png index b1c1c54..8a241a2 100644 Binary files a/Examples/html/example_HybridModel_02.png and b/Examples/html/example_HybridModel_02.png differ diff --git a/Examples/html/example_HybridModel_03.png b/Examples/html/example_HybridModel_03.png index 0bbeb76..9a983e8 100644 Binary files a/Examples/html/example_HybridModel_03.png and b/Examples/html/example_HybridModel_03.png differ diff --git a/Examples/html/example_HybridModel_04.png b/Examples/html/example_HybridModel_04.png index 663f242..d17cafc 100644 Binary files a/Examples/html/example_HybridModel_04.png and b/Examples/html/example_HybridModel_04.png differ diff --git a/Examples/html/example_HybridModel_05.png b/Examples/html/example_HybridModel_05.png index b288751..2927d6d 100644 Binary files a/Examples/html/example_HybridModel_05.png and b/Examples/html/example_HybridModel_05.png differ diff --git a/Examples/html/example_HybridModel_06.png b/Examples/html/example_HybridModel_06.png index 92d015c..9d51521 100644 Binary files a/Examples/html/example_HybridModel_06.png and b/Examples/html/example_HybridModel_06.png differ diff --git a/Examples/html/example_HybridModel_07.png b/Examples/html/example_HybridModel_07.png index d08ccc4..77149d9 100644 Binary files a/Examples/html/example_HybridModel_07.png and b/Examples/html/example_HybridModel_07.png differ diff --git a/Examples/html/example_HybridModel_08.png b/Examples/html/example_HybridModel_08.png index 4962f88..c1ecd1e 100644 Binary files a/Examples/html/example_HybridModel_08.png and b/Examples/html/example_HybridModel_08.png differ diff --git a/Examples/html/example_HybridModel_09.png b/Examples/html/example_HybridModel_09.png index 6f56d48..aa7855d 100644 Binary files a/Examples/html/example_HybridModel_09.png and b/Examples/html/example_HybridModel_09.png differ diff --git a/Examples/html/example_HybridModel_10.png b/Examples/html/example_HybridModel_10.png index 99a859b..aa61390 100644 Binary files a/Examples/html/example_HybridModel_10.png and b/Examples/html/example_HybridModel_10.png differ diff --git a/Examples/html/example_HybridModel_11.png b/Examples/html/example_HybridModel_11.png index 0c92d21..4de6c01 100644 Binary files a/Examples/html/example_HybridModel_11.png and b/Examples/html/example_HybridModel_11.png differ diff --git a/Examples/html/example_HybridModel_12.png b/Examples/html/example_HybridModel_12.png index c4bcae9..57d88f1 100644 Binary files a/Examples/html/example_HybridModel_12.png and b/Examples/html/example_HybridModel_12.png differ diff --git a/Examples/html/example_ModelReductionTool.html b/Examples/html/example_ModelReductionTool.html index a093f3f..b6f858d 100644 --- a/Examples/html/example_ModelReductionTool.html +++ b/Examples/html/example_ModelReductionTool.html @@ -6,7 +6,7 @@ example_ModelReductionTool

Contents

Using the SSIT to fit Multiple Models and Data sets with Shared Parameters

In this script, we show how multiple SSIT models and data sets can be fit simultaneously. This is most useful in situations where: 1) the analysis considers different experimental conditions (e.g., different time points, different inducer concentrations, different genetic mutations). 2) replica to replica variations are expected that would result in slightly different parameter combinations

close all
 clear all
-addpath('../CommandLine');
+addpath(genpath('../src'));
 

Define SSIT Model

SSIT models are defined as usual:

Model1 = SSIT;
 Model1.species = {'onGene';'rna'};
 Model1.initialCondition = [0;0];
@@ -87,6 +87,8 @@
 
 % We generate functions for model propensities
 Model1 = Model1.formPropensitiesGeneral('Model1FSP');
+
Starting parallel pool (parpool) using the 'Processes' profile ...
+Connected to parallel pool with 6 workers.
 

Create Second Model and associate to its own data

Model2 = Model1;
 Model2 = Model2.loadData('../ExampleData/DUSP1_Dex_100nM_Rep1_Rep2.csv',{'rna','RNA_nuc'},...
     {'Rep_num','2'}); % This would load the data assign onGene and rna and condition on Rep_num = 1;
@@ -208,7 +210,7 @@
 %   slightly different parameter combinations
 close all
 clear all
-addpath('../CommandLine');
+addpath(genpath('../src'));
 
 %% Define SSIT Model
 % SSIT models are defined as usual:
diff --git a/Examples/html/example_MultiModelTool.png b/Examples/html/example_MultiModelTool.png
index ed844e0..6ff23ed 100644
Binary files a/Examples/html/example_MultiModelTool.png and b/Examples/html/example_MultiModelTool.png differ
diff --git a/Examples/html/example_MultiModelTool_01.png b/Examples/html/example_MultiModelTool_01.png
index acdd44a..ad5c416 100644
Binary files a/Examples/html/example_MultiModelTool_01.png and b/Examples/html/example_MultiModelTool_01.png differ
diff --git a/Examples/html/example_MultiModelTool_02.png b/Examples/html/example_MultiModelTool_02.png
index 4177d1f..6d2f81b 100644
Binary files a/Examples/html/example_MultiModelTool_02.png and b/Examples/html/example_MultiModelTool_02.png differ
diff --git a/Examples/html/example_MultiModelTool_03.png b/Examples/html/example_MultiModelTool_03.png
index 70f14ec..d7d6a70 100644
Binary files a/Examples/html/example_MultiModelTool_03.png and b/Examples/html/example_MultiModelTool_03.png differ
diff --git a/Examples/html/example_MultiModelTool_04.png b/Examples/html/example_MultiModelTool_04.png
index 8812cc1..0d926cb 100644
Binary files a/Examples/html/example_MultiModelTool_04.png and b/Examples/html/example_MultiModelTool_04.png differ
diff --git a/Examples/html/example_MultiModelTool_05.png b/Examples/html/example_MultiModelTool_05.png
index abdfbe4..135fc1c 100644
Binary files a/Examples/html/example_MultiModelTool_05.png and b/Examples/html/example_MultiModelTool_05.png differ
diff --git a/Examples/html/example_MultiModelTool_06.png b/Examples/html/example_MultiModelTool_06.png
index 2a67128..0501551 100644
Binary files a/Examples/html/example_MultiModelTool_06.png and b/Examples/html/example_MultiModelTool_06.png differ
diff --git a/Examples/html/example_MultiModelTool_07.png b/Examples/html/example_MultiModelTool_07.png
index b1808d1..2b686e9 100644
Binary files a/Examples/html/example_MultiModelTool_07.png and b/Examples/html/example_MultiModelTool_07.png differ
diff --git a/Examples/html/example_MultiModelTool_08.png b/Examples/html/example_MultiModelTool_08.png
index 3a3849f..f7803ae 100644
Binary files a/Examples/html/example_MultiModelTool_08.png and b/Examples/html/example_MultiModelTool_08.png differ
diff --git a/Examples/html/example_MultiModelTool_09.png b/Examples/html/example_MultiModelTool_09.png
index e3ee06f..1d575d1 100644
Binary files a/Examples/html/example_MultiModelTool_09.png and b/Examples/html/example_MultiModelTool_09.png differ
diff --git a/Examples/html/example_MultiModelTool_10.png b/Examples/html/example_MultiModelTool_10.png
index 914ad38..d14424b 100644
Binary files a/Examples/html/example_MultiModelTool_10.png and b/Examples/html/example_MultiModelTool_10.png differ
diff --git a/Examples/html/example_MultiModelTool_11.png b/Examples/html/example_MultiModelTool_11.png
index 6cea4b3..eac2e45 100644
Binary files a/Examples/html/example_MultiModelTool_11.png and b/Examples/html/example_MultiModelTool_11.png differ
diff --git a/Examples/html/example_MultiModelTool_12.png b/Examples/html/example_MultiModelTool_12.png
index ae39395..78effe3 100644
Binary files a/Examples/html/example_MultiModelTool_12.png and b/Examples/html/example_MultiModelTool_12.png differ
diff --git a/Examples/html/example_MultiModelTool_13.png b/Examples/html/example_MultiModelTool_13.png
index 73eaa4c..1f1c666 100644
Binary files a/Examples/html/example_MultiModelTool_13.png and b/Examples/html/example_MultiModelTool_13.png differ
diff --git a/Examples/html/example_MultiModelTool_14.png b/Examples/html/example_MultiModelTool_14.png
index 700b522..f65265a 100644
Binary files a/Examples/html/example_MultiModelTool_14.png and b/Examples/html/example_MultiModelTool_14.png differ
diff --git a/Examples/html/example_MultiModelTool_15.png b/Examples/html/example_MultiModelTool_15.png
index fc60324..10d9f48 100644
Binary files a/Examples/html/example_MultiModelTool_15.png and b/Examples/html/example_MultiModelTool_15.png differ
diff --git a/Examples/html/example_MultiModelTool_16.png b/Examples/html/example_MultiModelTool_16.png
index 7297a9f..d134ce9 100644
Binary files a/Examples/html/example_MultiModelTool_16.png and b/Examples/html/example_MultiModelTool_16.png differ
diff --git a/Examples/html/example_MultiModelTool_17.png b/Examples/html/example_MultiModelTool_17.png
index 646a196..3342c14 100644
Binary files a/Examples/html/example_MultiModelTool_17.png and b/Examples/html/example_MultiModelTool_17.png differ
diff --git a/Examples/html/example_MultiModelTool_18.png b/Examples/html/example_MultiModelTool_18.png
index c9acf95..10b37be 100644
Binary files a/Examples/html/example_MultiModelTool_18.png and b/Examples/html/example_MultiModelTool_18.png differ
diff --git a/Examples/html/example_MultiModelTool_19.png b/Examples/html/example_MultiModelTool_19.png
index 3361851..fc4fa11 100644
Binary files a/Examples/html/example_MultiModelTool_19.png and b/Examples/html/example_MultiModelTool_19.png differ
diff --git a/Examples/html/example_MultiModelTool_20.png b/Examples/html/example_MultiModelTool_20.png
index 3021a89..cf66713 100644
Binary files a/Examples/html/example_MultiModelTool_20.png and b/Examples/html/example_MultiModelTool_20.png differ
diff --git a/Examples/html/example_MultiModelTool_21.png b/Examples/html/example_MultiModelTool_21.png
index 0cffbfa..43f34e0 100644
Binary files a/Examples/html/example_MultiModelTool_21.png and b/Examples/html/example_MultiModelTool_21.png differ
diff --git a/Examples/html/example_MultiModelTool_22.png b/Examples/html/example_MultiModelTool_22.png
index 2222b8c..0e6cc0c 100644
Binary files a/Examples/html/example_MultiModelTool_22.png and b/Examples/html/example_MultiModelTool_22.png differ
diff --git a/Examples/html/example_MultiModelTool_23.png b/Examples/html/example_MultiModelTool_23.png
index b899118..ed2e216 100644
Binary files a/Examples/html/example_MultiModelTool_23.png and b/Examples/html/example_MultiModelTool_23.png differ
diff --git a/Examples/html/example_MultiModelTool_24.png b/Examples/html/example_MultiModelTool_24.png
index 4ce9adf..68f3c5d 100644
Binary files a/Examples/html/example_MultiModelTool_24.png and b/Examples/html/example_MultiModelTool_24.png differ
diff --git a/Examples/html/example_MultiModelTool_25.png b/Examples/html/example_MultiModelTool_25.png
index 08d09fa..554391c 100644
Binary files a/Examples/html/example_MultiModelTool_25.png and b/Examples/html/example_MultiModelTool_25.png differ
diff --git a/Examples/html/example_MultiModelTool_26.png b/Examples/html/example_MultiModelTool_26.png
index 4c2b678..5935cb3 100644
Binary files a/Examples/html/example_MultiModelTool_26.png and b/Examples/html/example_MultiModelTool_26.png differ
diff --git a/Examples/html/example_MultiModelTool_27.png b/Examples/html/example_MultiModelTool_27.png
index 3a6fd1e..3fb63a0 100644
Binary files a/Examples/html/example_MultiModelTool_27.png and b/Examples/html/example_MultiModelTool_27.png differ
diff --git a/Examples/html/example_MultiModelTool_28.png b/Examples/html/example_MultiModelTool_28.png
index 73ae4c9..484d3c6 100644
Binary files a/Examples/html/example_MultiModelTool_28.png and b/Examples/html/example_MultiModelTool_28.png differ
diff --git a/Examples/html/example_MultiModelTool_29.png b/Examples/html/example_MultiModelTool_29.png
index 797406c..b60230f 100644
Binary files a/Examples/html/example_MultiModelTool_29.png and b/Examples/html/example_MultiModelTool_29.png differ
diff --git a/Examples/html/example_MultiModelTool_30.png b/Examples/html/example_MultiModelTool_30.png
index dac231b..8a11130 100644
Binary files a/Examples/html/example_MultiModelTool_30.png and b/Examples/html/example_MultiModelTool_30.png differ
diff --git a/Examples/html/example_MultiModelTool_31.png b/Examples/html/example_MultiModelTool_31.png
index e708baa..fdcf3e5 100644
Binary files a/Examples/html/example_MultiModelTool_31.png and b/Examples/html/example_MultiModelTool_31.png differ
diff --git a/Examples/html/example_MultiModelTool_32.png b/Examples/html/example_MultiModelTool_32.png
index b165534..53f1a48 100644
Binary files a/Examples/html/example_MultiModelTool_32.png and b/Examples/html/example_MultiModelTool_32.png differ
diff --git a/Examples/html/example_MultiModelTool_33.png b/Examples/html/example_MultiModelTool_33.png
index 1e1e62d..50e6439 100644
Binary files a/Examples/html/example_MultiModelTool_33.png and b/Examples/html/example_MultiModelTool_33.png differ
diff --git a/Examples/html/example_MultiModelTool_34.png b/Examples/html/example_MultiModelTool_34.png
index 2e1b4df..a5b99da 100644
Binary files a/Examples/html/example_MultiModelTool_34.png and b/Examples/html/example_MultiModelTool_34.png differ
diff --git a/Examples/html/example_MultiModelTool_35.png b/Examples/html/example_MultiModelTool_35.png
index b12b193..b35294e 100644
Binary files a/Examples/html/example_MultiModelTool_35.png and b/Examples/html/example_MultiModelTool_35.png differ
diff --git a/Examples/html/example_MultiModelTool_36.png b/Examples/html/example_MultiModelTool_36.png
index 242bb97..9dc66e3 100644
Binary files a/Examples/html/example_MultiModelTool_36.png and b/Examples/html/example_MultiModelTool_36.png differ
diff --git a/Examples/html/example_MultiModelTool_37.png b/Examples/html/example_MultiModelTool_37.png
index 140f54d..94dc2b5 100644
Binary files a/Examples/html/example_MultiModelTool_37.png and b/Examples/html/example_MultiModelTool_37.png differ
diff --git a/Examples/html/example_MultiModelTool_38.png b/Examples/html/example_MultiModelTool_38.png
index 75ecbac..9f991ad 100644
Binary files a/Examples/html/example_MultiModelTool_38.png and b/Examples/html/example_MultiModelTool_38.png differ
diff --git a/Examples/html/example_MultiModelTool_39.png b/Examples/html/example_MultiModelTool_39.png
index 933cc61..b2084ed 100644
Binary files a/Examples/html/example_MultiModelTool_39.png and b/Examples/html/example_MultiModelTool_39.png differ
diff --git a/Examples/html/example_MultiModelTool_40.png b/Examples/html/example_MultiModelTool_40.png
index 66da5c1..a11dbfa 100644
Binary files a/Examples/html/example_MultiModelTool_40.png and b/Examples/html/example_MultiModelTool_40.png differ
diff --git a/Examples/html/example_MultiModelTool_41.png b/Examples/html/example_MultiModelTool_41.png
index 59b829e..4261900 100644
Binary files a/Examples/html/example_MultiModelTool_41.png and b/Examples/html/example_MultiModelTool_41.png differ
diff --git a/Examples/html/example_MultiModelTool_42.png b/Examples/html/example_MultiModelTool_42.png
index 15e7e39..f02247f 100644
Binary files a/Examples/html/example_MultiModelTool_42.png and b/Examples/html/example_MultiModelTool_42.png differ
diff --git a/Examples/html/example_MultiModelTool_43.png b/Examples/html/example_MultiModelTool_43.png
index 96e6987..a4a4163 100644
Binary files a/Examples/html/example_MultiModelTool_43.png and b/Examples/html/example_MultiModelTool_43.png differ
diff --git a/Examples/html/example_MultiModelTool_44.png b/Examples/html/example_MultiModelTool_44.png
index a90cacf..d658981 100644
Binary files a/Examples/html/example_MultiModelTool_44.png and b/Examples/html/example_MultiModelTool_44.png differ
diff --git a/Examples/html/example_SSITBasics.html b/Examples/html/example_SSITBasics.html
index 6530a2c..1e3cf3f 100644
--- a/Examples/html/example_SSITBasics.html
+++ b/Examples/html/example_SSITBasics.html
@@ -6,7 +6,7 @@
    example_SSITBasics

Contents

clear all
 close all
+addpath(genpath('../src'));
 

(1) Creating and Modifying SSIT Model

(1A) Choosing a Pre-Made Model

ModelChoice = 'BirthDeath';    % One species problem
 % ModelChoice = 'ToggleSwitch';  % Two species problem (non-linear toggle switch)
 % ModelChoice = 'CentralDogmaTV';  % Two species problem (mRNa and protein) with time varying transcription rate
@@ -84,6 +85,8 @@
 F1.initialTime = -1;
 F1.solutionScheme = 'FSP';    % Set solutions scheme to FSP.
 [FSPsoln,F1.fspOptions.bounds] = F1.solve;  % Solve the FSP analysis
+
Starting parallel pool (parpool) using the 'Processes' profile ...
+Connected to parallel pool with 6 workers.
 

(2A.1) Solve FSP Model again using the bounds from the last solution

If we start with the bounds computed in the first analysis, the solution is often much faster.

[FSPsoln] = F1.solve(FSPsoln.stateSpace);  % Solve the FSP analysis
 

(2A.2) Make plots of FSP solution

F1.makePlot(FSPsoln,'meansAndDevs',[],[],1,{'linewidth',3,'color',[0,1,1]}) % Make plot of mean vs. time.
 F1.makePlot(FSPsoln,'marginals',[],[],2,{'linewidth',3,'color',[0,0,1]}) % Make plot of mean vs. time.
@@ -97,15 +100,14 @@
 

(3) Sensitivity Analysis using FSP

F4 = F1;
 F4.solutionScheme = 'fspSens'; % Set solutions scheme to FSP Sensitivity
 [sensSoln,bounds] = F4.solve(FSPsoln.stateSpace);  % Solve the sensitivity problem
-
Error with Analytical Sensitivity Calculations - Switching to Finite Difference Method
 

(3A) Make plot of sensitivities

F4.makePlot(sensSoln,'marginals',[],[],4,{'linewidth',3,'color',[0,0,1]}) % Plot marginal sensitivities