Skip to content

Commit 60cf301

Browse files
alexpopingaalexpopinga
committed
Adding comments and testing prior and FIM effects
1 parent e649345 commit 60cf301

File tree

1 file changed

+216
-44
lines changed

1 file changed

+216
-44
lines changed

‎WorkSpace/EricModel/a0_Fit_GR_and_DUSP1_models.m

+216-44
Original file line numberDiff line numberDiff line change
@@ -13,9 +13,10 @@
1313
savedWorkspace = 'workspaceJuly24';
1414
addpath('tmpPropensityFunctions');
1515

16-
%% STEP 0 -- Preliminaries.
16+
%% GR model
17+
%% STEP 0 -- Preliminaries: Load previously computed data or define GR model from scratch
1718
if loadPrevious
18-
% STEP 0.A -- Option to load previous workspace from file -- just for testing.
19+
% STEP 0.A. -- Option to load previous workspace from file -- just for testing.
1920
varNames = {'ModelGR',
2021
'GRfitCases'
2122
'log10PriorMean'
@@ -39,27 +40,41 @@
3940
GRpars
4041
%GRpars = GRpars';
4142
else
43+
%% STEP 0.B. -- Create Base Model for GR Only
44+
% Here, I set up the model for the GR translocation dynamics.
4245
fitOptions = optimset('Display','iter','MaxIter',300);
43-
%% STEP 0.B.1. -- Create Base Model for GR Only
44-
% Here, I set up the model for the GR translocation dynamics.
46+
47+
% Create blank SSIT model.
4548
ModelGR = SSIT;
49+
50+
% Set species names.
4651
ModelGR.species = {'cytGR';'nucGR'};
52+
53+
% Set initial condition.
4754
ModelGR.initialCondition = [20;1];
55+
56+
% Define propensity functions and input signals:
4857
ModelGR.propensityFunctions = {'(kcn0 + (t>0)*kcn1*IDex/(MDex+IDex)) * cytGR';'knc*nucGR';...
4958
'kg1';'gg1*cytGR';'gg2*nucGR'};
59+
60+
% Define stoichiometry.
5061
ModelGR.stoichiometry = [-1,1,1,-1,0;...
5162
1,-1,0,0,-1];
63+
64+
% Specify parameter guesses.
5265
ModelGR.parameters = ({'koff',0.1;'kon',0.1;'kr',1;'gr',0.02;...
5366
'kcn0',0.005;'kcn1',0.02;'gDex',0.003;'knc',0.01;'kg1',14e-5;...
5467
'gg1',1e-5;'gg2',1e-6;'MDex',5;'Dex0',100});
68+
69+
5570
log10PriorMean = [-1 -1 0 -2,...
5671
-1 -3 -2 -1 -2 -2 -2 0.5, 2];
5772
log10PriorStd = 2*ones(1,13);
5873
% the log prior will be applied to the fit to multiple models as an
5974
% additional constraint.
60-
6175
ModelGR.fittingOptions.logPrior = [];
6276
% So it is left out of the prior, since we only want it to be calculated once.
77+
6378
ModelGR.fspOptions.initApproxSS = true;
6479
ModelGR.fittingOptions.modelVarsToFit = (5:12);
6580

@@ -70,10 +85,12 @@
7085

7186
[FSPGrSoln,ModelGR.fspOptions.bounds] = ModelGR.solve;
7287
[FSPGrSoln,ModelGR.fspOptions.bounds] = ModelGR.solve(FSPGrSoln.stateSpace);
73-
%% STEP 0.B.2. -- Load previously fit parameter values (optional)
88+
89+
% STEP 0.B.2. -- Define GR parameters
7490
fitIters = 3;
75-
GRpars = cell2mat(ModelGR.parameters(5:12,2))';
76-
%% STEP 0.B.3. -- Associate GR Data with Different Instances of Model (10,100nm Dex)
91+
GRpars = cell2mat(ModelGR.parameters(5:12,2))';
92+
93+
% STEP 0.B.3. -- Associate GR Data with Different Instances of Model (10,100nm Dex)
7794
GRfitCases = {'1','1',101,'GR Fit (1nM Dex)';...
7895
'10','10',102,'GR Fit (10nM Dex)';...
7996
'100','100',103,'GR Fit (100nM Dex)'};
@@ -90,7 +107,8 @@
90107
% the entire class of models. In this case, these refer to
91108
% global parameters 5:15 (GR parameters).
92109
end
93-
%% STEP 0.B.4. -- Make Guesses for the FSP bounds
110+
111+
% STEP 0.B.4. -- Make Guesses for the FSP bounds
94112
% This is sometimes necessary when using an uninduced steady state as the
95113
% initial condition. You need to guess a reasonable statespace or the
96114
% computation of the SS can be inaccurate.
@@ -101,52 +119,38 @@
101119
% custom.
102120
end
103121
end
104-
%%
105-
%% STEP 1 -- GR Model
106-
%% STEP 1.A. -- Combine all three GR models and fit using a single parameter set.
122+
123+
%% STEP 1 -- Fit GR Models
124+
% STEP 1.A. -- Specify dataset time points.
107125
for i = 1:3
108126
ModelGRfit{i}.tSpan = ModelGRfit{i}.dataSet.times;
109127
end
110-
%%
111-
%% STEP 1.A. -- Run MH on GR Models.
112-
%GRpars = GRpars';
113-
%MHFitOptions.thin=1;
114-
%MHFitOptions.numberOfSamples=1000;
115-
%MHFitOptions.burnIn=0;
116-
%MHFitOptions.progress=true;
117-
%MHFitOptions.numChains = 1;
118-
%MHFitOptions.useFIMforMetHast = true;
119-
%[~,~,MHResultsGR] = ModelGRfit{1}.maximizeLikelihood(...
120-
% GRpars, MHFitOptions, 'MetropolisHastings');
121-
%%
122-
%figNew = figure;
123-
%ModelGR.plotMHResults(MHResultsGR,[],'log',[],figNew)
124-
%for i = 1:7
125-
% for j = i+1:7
126-
% subplot(7,7,(i-1)*7+j-1)
127-
% CH = get(gca,'Children');
128-
% CH(1).Color=[1,0,1]; %
129-
% CH(1).LineWidth = 3;
130-
% end
131-
%end
132-
%%
133-
logPriorGR = @(x)-sum((log10(x)-log10PriorMean(5:12)').^2./(2*log10PriorStd(5:12)'.^2));
128+
129+
% STEP 1.B. -- Specify log prior (NOTE: must transpose due to Matlab update that
130+
% no longer correctly assumes format when adding single value vector to
131+
% column vector).
132+
%logPriorGR = @(x)-sum((log10(x)-log10PriorMean(5:12)').^2./(2*log10PriorStd(5:12)'.^2));
133+
134+
% STEP 1.C. -- Combine all three GR models and fit using a single parameter set.
134135
for jj = 1:fitIters
135-
combinedGRModel = SSITMultiModel(ModelGRfit,ModelGRparameterMap,logPriorGR);
136+
%combinedGRModel = SSITMultiModel(ModelGRfit,ModelGRparameterMap,logPriorGR);
137+
combinedGRModel = SSITMultiModel(ModelGRfit,ModelGRparameterMap);
136138
combinedGRModel = combinedGRModel.initializeStateSpaces(boundGuesses);
137139
combinedGRModel = combinedGRModel.updateModels(GRpars,false);
138140
GRpars = combinedGRModel.maximizeLikelihood(...
139141
GRpars, fitOptions);
140142
save('EricModel_MMDex','GRpars')
141143
end
142-
%% STEP 1.B. -- Compute FIM for GR parameters.
143-
%combinedGRModel = combinedGRModel.computeFIMs([],'log');
144-
%fimGR_withPrior = combinedGRModel.FIM.totalFIM+... % the FIM in log space.
145-
% diag(1./(log10PriorStd(ModelGR.fittingOptions.modelVarsToFit)*log(10)).^2); % Add prior in log space.
146-
%% STEP 1.C. -- Run MH on GR Models.
144+
145+
%% STEP 1.D. -- Compute FIM for GR parameters.
146+
combinedGRModel = combinedGRModel.computeFIMs([],'log');
147+
fimGR_withPrior = combinedGRModel.FIM.totalFIM+... % the FIM in log space.
148+
diag(1./(log10PriorStd(ModelGR.fittingOptions.modelVarsToFit)*log(10)).^2); % Add prior in log space.
149+
150+
%% STEP 1.C. -- Run MH on GR Models.
147151
%GRpars = GRpars';
148-
MHFitOptions.thin=3;
149-
MHFitOptions.numberOfSamples=3000;
152+
MHFitOptions.thin=1;
153+
MHFitOptions.numberOfSamples=1000;
150154
MHFitOptions.burnIn=0;
151155
MHFitOptions.progress=true;
152156
MHFitOptions.numChains = 1;
@@ -167,6 +171,174 @@
167171
end
168172
end
169173

174+
%% STEP 1.D. -- Make Plots of GR Fit Results
175+
makeGRPlots(combinedGRModel,GRpars)
176+
177+
%% STEP 2 -- Extend Model to Include DUSP1 Activation, Production, and Degradation
178+
if loadPrevious
179+
varNamesDUSP1 = {'ModelGRDusp100nM'
180+
'GRfitCases'
181+
'log10PriorMean'
182+
'log10PriorStd'
183+
'duspLogPrior'
184+
'DUSP1pars'
185+
'Dusp1FitCases'};
186+
load(savedWorkspace,varNamesDUSP1{:})
187+
fitOptions = optimset('Display','iter','MaxIter',10);
188+
fitIters = 1;
189+
try
190+
ModelGRDusp100nM.propensitiesGeneral{1}.stateDependentFactor(0);
191+
catch
192+
ModelGRDusp100nM = ModelGRDusp100nM.formPropensitiesGeneral('DUSP1_Model');
193+
end
194+
else
195+
%% STEP 2.A.1. -- Add DUSP1 to the existing GR model.
196+
% Copy parameters from the 100nM Dex stim case in GR.
197+
fitIters = 3;
198+
ModelGRDusp100nM = ModelGRfit{3};
199+
ModelGRDusp100nM.species = {'offGene';'onGene';'cytGR';'nucGR';'rna'};
200+
ModelGRDusp100nM.initialCondition = [2;0;24;1;5];
201+
ModelGRDusp100nM.propensityFunctions = {'kon*offGene*nucGR';'koff*onGene';
202+
'(kcn0 + (t>0)*kcn1*IDex/(MDex+IDex)) * cytGR';'knc*nucGR';'kg1';'gg1*cytGR';'gg2*nucGR';...
203+
'kr*onGene';'gr*rna'};
204+
ModelGRDusp100nM.stoichiometry = [-1,1,0,0,0,0,0,0,0;...
205+
1,-1,0,0,0,0,0,0,0;...
206+
0,0,-1,1,1,-1,0,0,0;...
207+
0,0,1,-1,0,0,-1,0,0;...
208+
0,0,0,0,0,0,0,1,-1];
209+
ModelGRDusp100nM.useHybrid = true;
210+
ModelGRDusp100nM.hybridOptions.upstreamODEs = {'cytGR','nucGR'};
211+
ModelGRDusp100nM.solutionScheme = 'FSP';
212+
ModelGRDusp100nM.customConstraintFuns = [];
213+
ModelGRDusp100nM.fspOptions.bounds = [0;0;0;2;2;400];
214+
ModelGRDusp100nM.fittingOptions.modelVarsToFit = 1:4;
215+
ModelGRDusp100nM = ModelGRDusp100nM.formPropensitiesGeneral('EricModDusp1');
216+
duspLogPrior = @(x)-sum((log10(x(:))'-log10PriorMean(1:4)).^2./(2*log10PriorStd(1:4).^2));
217+
ModelGRDusp100nM.fittingOptions.logPrior = duspLogPrior;
218+
219+
%% STEP 2.A.2. -- Load pre-fit parameters into model.
220+
if loadPrevious
221+
load('EricModelDusp1_MMDex','DUSP1pars')
222+
else
223+
%% Pull the DUSP1 parameters from the Model
224+
% Find the indices of the desired parameter names
225+
knc_idx = find(strcmp(ModelGRDusp100nM.parameters(:,1), 'knc'));
226+
kg1_idx = find(strcmp(ModelGRDusp100nM.parameters(:,1), 'kg1'));
227+
gg1_idx = find(strcmp(ModelGRDusp100nM.parameters(:,1), 'gg1'));
228+
gg2_idx = find(strcmp(ModelGRDusp100nM.parameters(:,1), 'gg2'));
229+
230+
% Extract the values and store them in DUSP1pars
231+
DUSP1pars = [ModelGRDusp100nM.parameters{knc_idx, 2}, ...
232+
ModelGRDusp100nM.parameters{kg1_idx, 2}, ...
233+
ModelGRDusp100nM.parameters{gg1_idx, 2}, ...
234+
ModelGRDusp100nM.parameters{gg2_idx, 2}];
235+
end
236+
ModelGRDusp100nM.parameters(ModelGRDusp100nM.fittingOptions.modelVarsToFit,2) = num2cell(DUSP1pars);
237+
%% STEP 2.A.3. -- Load and Associate with DUSP1 smFISH Data (100nM Dex Only)
238+
% The commented code below would be needed to fit multiple conditions,
239+
% but that is not used in this case. It is left here in case it is
240+
% needed in later stages of the project.
241+
% % % Dusp1FitCases = {'100','100',201,'DUSP1 Fit (100nM Dex)'};
242+
% % % ModelDusp1Fit = cell(size(Dusp1FitCases,1),1);
243+
% % % ModelDusp1parameterMap = cell(1,size(GRfitCases,1));
244+
% % % for i = 1:size(Dusp1FitCases,1)
245+
% % % ModelDusp1Fit{i} = ModelGRDusp100nM.loadData('EricData/pdoCalibrationData_EricIntensity_DexSweeps.csv',...
246+
% % % {'rna','totalNucRNA'},...
247+
% % % {'Dex_Conc','100'});
248+
% % % ModelDusp1Fit{i}.inputExpressions = {'IDex','Dex0*exp(-gDex*t)'};
249+
% % %
250+
% % % ModelDusp1parameterMap{i} = (1:4);
251+
% % % % Set Dex concentration.
252+
% % % ModelDusp1Fit{i}.parameters{13,2} = str2num(Dusp1FitCases{i,1});
253+
% % % ModelDusp1Fit{i} = ModelDusp1Fit{i}.formPropensitiesGeneral(['EricModDusp1_',num2str(i),'_FSP']);
254+
% % % end
255+
% % % DUSP1pars = [ModelDusp1Fit{i}.parameters{ModelGRDusp100nM.fittingOptions.modelVarsToFit,2}];
256+
ModelGRDusp100nM = ModelGRDusp100nM.loadData('EricData/pdoCalibrationData_EricIntensity_DexSweeps.csv',...
257+
{'rna','totalNucRNA'},{'Dex_Conc','100'});
258+
DUSP1pars = [ModelGRDusp100nM.parameters{ModelGRDusp100nM.fittingOptions.modelVarsToFit,2}];
259+
end
260+
%% STEP 2.B. -- Fit DUSP1 model at 100nM Dex.
261+
for i = 1:fitIters
262+
fitOptions.suppressFSPExpansion = true;
263+
DUSP1pars = ModelGRDusp100nM.maximizeLikelihood(...
264+
DUSP1pars, fitOptions);
265+
ModelGRDusp100nM.parameters(1:4,2) = num2cell(DUSP1pars);
266+
save('EricModelDusp1_MMDex','GRpars','DUSP1pars')
267+
end
268+
%% STEP 2.C. -- Plot predictions for other Dex concentrations.
269+
showCases = [1,1,1,1];
270+
makePlotsDUSP1({ModelGRDusp100nM},ModelGRDusp100nM,DUSP1pars,Dusp1FitCases,showCases)
271+
272+
%% STEP 2.D. -- Sample uncertainty for Dusp1 Parameters
273+
% STEP 2.D.1. -- Compute sensitivity of the FSP solution
274+
ModelGRDusp100nM.solutionScheme = 'fspSens';
275+
sensSoln = ModelGRDusp100nM.solve();
276+
ModelGRDusp100nM.solutionScheme = 'FSP';
277+
% STEP 2.D.2. -- Compute FIM
278+
% define which species in model are not observed.
279+
ModelGRDusp100nM.pdoOptions.unobservedSpecies = {'offGene';'onGene'};
280+
% compute the FIM
281+
fimResults = ModelGRDusp100nM.computeFIM(sensSoln.sens,'log');
282+
% In the following, the log-prior is used as a prior co-variance matrix.
283+
% This will be used in the FIM calculation as an FIM without new evidence
284+
% being set equal to the inverse of this covariance matrix. More rigorous
285+
% justification is needed to support this heuristic.
286+
fimTotal = ModelGRDusp100nM.evaluateExperiment(fimResults,ModelGRDusp100nM.dataSet.nCells,...
287+
diag(log10PriorStd(1:13).^2));
288+
FIMfree = fimTotal{1}(1:4,1:4);
289+
if min(eig(FIMfree))<1
290+
disp('Warning -- FIM has one or more small eigenvalues. Reducing proposal width to 10x in those directions. MH Convergence may be slow.')
291+
FIMfree = FIMfree + 1*eye(length(FIMfree));
292+
end
293+
covFree = FIMfree^-1;
294+
covFree = 0.5*(covFree+covFree');
295+
296+
%% STEP 2.D.3. -- Run Metropolis Hastings Search
297+
if loadPrevious
298+
MHDusp1File = 'MHDusp1_Jul22';
299+
load(MHDusp1File)
300+
else
301+
MHFitOptions.proposalDistribution=@(x)mvnrnd(x,covFree);
302+
MHFitOptions.thin=1;
303+
MHFitOptions.numberOfSamples=1000;
304+
MHFitOptions.burnIn=0;
305+
MHFitOptions.progress=true;
306+
MHFitOptions.numChains = 1;
307+
MHFitOptions.saveFile = 'TMPEricMHDusp1.mat';
308+
[DUSP1pars,~,MHResultsDusp1] = ModelGRDusp100nM.maximizeLikelihood(...
309+
[], MHFitOptions, 'MetropolisHastings');
310+
delete('TMPEricMHDusp1.mat')
311+
ModelGRDusp100nM.parameters(1:4,2) = num2cell(DUSP1pars);
312+
end
313+
%% STEP 1.B. -- Compute FIM for GR parameters.
314+
combinedGRModel = combinedGRModel.computeFIMs([],'log');
315+
fimGR_withPrior = combinedGRModel.FIM.totalFIM+... % the FIM in log space.
316+
diag(1./(log10PriorStd(ModelGR.fittingOptions.modelVarsToFit)*log(10)).^2); % Add prior in log space.
317+
%% STEP 1.C. -- Run MH on GR Models.
318+
%GRpars = GRpars';
319+
MHFitOptions.proposalDistribution=@(x)mvnrnd(x,covFree);
320+
MHFitOptions.thin=1;
321+
MHFitOptions.numberOfSamples=1000;
322+
MHFitOptions.burnIn=0;
323+
MHFitOptions.progress=true;
324+
MHFitOptions.numChains = 1;
325+
MHFitOptions.useFIMforMetHast = true;
326+
MHFitOptions.saveFile = 'TMPEricMHGR.mat';
327+
[GRpars,~,MHResultsGR] = combinedGRModel.maximizeLikelihood(...
328+
GRpars, MHFitOptions, 'MetropolisHastings');
329+
delete(MHFitOptions.saveFile)
330+
%%
331+
figNew = figure;
332+
ModelGR.plotMHResults(MHResultsGR,[],'log',[],figNew)
333+
for i = 1:7
334+
for j = i+1:7
335+
subplot(7,7,(i-1)*7+j-1)
336+
CH = get(gca,'Children');
337+
CH(1).Color=[1,0,1]; %
338+
CH(1).LineWidth = 3;
339+
end
340+
end
341+
170342
%% STEP 1.D. -- Make Plots of GR Fit Results
171343
makeGRPlots(combinedGRModel,GRpars)
172344
%%

0 commit comments

Comments
 (0)