1
+ %% Initial Setup
2
+
3
+ close all
4
+ clear all
5
+ addpath(genpath(' ../../src' ));
6
+
7
+ M = [1 2 3 6 ]; % Number of intermediate states
8
+
9
+ kOn = 0.1 ;
10
+ kOff = kOn * 1000 / 2 ;
11
+
12
+ tau = 50 ;
13
+ delayPeriods = 10 ;
14
+ tSpan = linspace(0 , tau * delayPeriods , 1 + delayPeriods );
15
+
16
+ Model = SSIT ;
17
+ Model.species = {' Dempty' ; ' Dfull' ; ' X' };
18
+
19
+ initialConditionX = 44 ; % Result from actual SSA runs
20
+ Model.initialCondition = [1 ; 0 ; initialConditionX ];
21
+
22
+ parameters = ({...
23
+ ' kPlus' , 100 ; ... % (Delayed) rate of protein production
24
+ ' kMinus' , 1 ; ... % Rate of protein degradation
25
+ ' kOn' , kOn ; ... % Rate of operator site emptying / activation
26
+ ' kOff' , kOff ; ... % Rate of operator site filling / deactivation
27
+ ' tau' , tau ... % Delay in protein production
28
+ });
29
+ Model.parameters = parameters ;
30
+
31
+ Model.propensityFunctions = {
32
+ ' kOn*Dfull' , ... % Operator site emptying
33
+ ' kOff*Dempty*X' , ... % Operator site filling
34
+ ' kPlus*Dempty' , ... % (Delayed) protein production
35
+ ' kMinus*X' ... % Protein degradation
36
+ };
37
+
38
+ Model.stoichiometry = [...
39
+ 1 , - 1 , 0 , 0 ; ... % Dempty
40
+ - 1 , 1 , 0 , 0 ; ... % Dfull
41
+ 0 , 0 , 1 , - 1 ... % X
42
+ ];
43
+
44
+ % The only delayed reaction in the system is the third, with a delay equal
45
+ % to tau. It has a single reactant, one Dempty, so we must include it in
46
+ % the delayed-reaction scheduling stochiometries:
47
+
48
+ Model.delayedReactions = [3 , Model.parameters{5 ,2 }];
49
+
50
+ Model.delayedReactionSchedulingS = [...
51
+ 0 , 0 , - 1 , 0 ; ... % Dempty
52
+ 0 , 0 , 0 , 0 ; ... % Dfull
53
+ 0 , 0 , 0 , 0 ... % X
54
+ ];
55
+
56
+ datasetSizes = [1 , 2 , 5 , 10 , 15 , 20 , 30 ] * 1e3 ;
57
+ datasetFractions = datasetSizes / 3e4 ;
58
+ for dsCntr = 1 : length(datasetSizes )
59
+ datasetNames{dsCntr } = [' ssa_' , num2str(datasetSizes(dsCntr )), ' .csv' ];
60
+ end
61
+
62
+ fspTimes = zeros(length(M ), 1 );
63
+ fitTimes = zeros(length(M ), length(datasetSizes ));
64
+
65
+ %% %% Protein intermediates
66
+
67
+ paramsToFit = [1 , 4 , 5 ]; % kPlus, kOff, tau
68
+
69
+ fitParamValues = zeros(length(M ), length(datasetSizes ), length(paramsToFit ));
70
+
71
+ for Mcntr = 1 : length(M )
72
+ curM = M(Mcntr );
73
+ piModel = getNFModelWithProteinIntermediates(...
74
+ curM , initialConditionX , parameters );
75
+ piModel.tSpan = tSpan ;
76
+ piModel.solutionScheme = ' FSP' ;
77
+ tic ;
78
+ [piFSPsoln , piModel .fspOptions .bounds ] = piModel .solve ;
79
+ fspTimes(Mcntr ) = toc ;
80
+
81
+ for dsCntr = 1 : length(datasetSizes )
82
+ curModel = piModel ; % Clone the model solved via FSP
83
+
84
+ % Load the subset of the SSA data
85
+
86
+ curModel = curModel .loadData(datasetNames{dsCntr }, ...
87
+ {' Dempty' ,' exp1_s1' ,' Dfull' ,' exp1_s2' ,' X' ,' exp1_s3' });
88
+
89
+ % Only a subset of the parameters will be made "free"
90
+
91
+ curModel.fittingOptions.modelVarsToFit = paramsToFit ;
92
+
93
+ fitOptions = optimset(' Display' , ' none' , ' MaxIter' , 1000 );
94
+ tic ;
95
+ fitParameters = curModel .maximizeLikelihood([], fitOptions );
96
+ fitTimes(Mcntr ,dsCntr ) = toc ;
97
+ curModel .parameters(curModel .fittingOptions .modelVarsToFit , 2 ) = ...
98
+ num2cell(fitParameters );
99
+ for paramCntr = 1 : length(paramsToFit )
100
+ fitParamValues(Mcntr , dsCntr , paramCntr ) = ...
101
+ curModel .parameters(paramsToFit(paramCntr ), 2 );
102
+ end
103
+ end % for dsCntr = 1:length(datasetSizes) ...
104
+
105
+ % piModel.makePlot(piFSPsoln,'meansAndDevs',[],[]) % Make plot of mean vs. time.
106
+ % piModel.makePlot(piFSPsoln,'marginals',[],[]) % Make plot of marginals
107
+ end % for Mcntr = 1:length(M)
108
+
109
+ disp(fspTimes )
110
+ save piFSPtimes.mat fspTimes
111
+ disp(fitTimes )
112
+ save piFitTimes.mat fitTimes
113
+ disp(fitParamValues )
114
+ save piFitParamValues.mat fitParamValues
0 commit comments