Skip to content

Commit 6231139

Browse files
alexpopingaalexpopinga
committed
Plot 3 Dex conc data on GR-beta model for comparison with GR-alpha
1 parent f2e7d9e commit 6231139

File tree

2 files changed

+86
-49
lines changed

2 files changed

+86
-49
lines changed

‎WorkSpace/EricModel/GR_2_Electric_Boogaloo.m

+81-44
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@
1616
loadPrevious = false;
1717
savedWorkspace = 'workspace_Boogaloo';
1818

19-
fitOptions = optimset('Display','iter','MaxIter',300);
19+
fitOptions = optimset('Display','iter','MaxIter',200);
2020

2121
%%
2222
GR = input('(1) Convolve GR-alpha + GR-beta with FSP;\n(2) ODES for GR-alpha + GR-beta;\n(3) SSAs\nChoose your destiny: ');
@@ -92,7 +92,7 @@
9292
data = readtable('../EricModel/EricData/GR_ALL_gated_with_CytoArea_and_normGR_Feb2825_03.csv');
9393

9494
% Filter the data if desired for particular Time_index and/or dex_conc
95-
filteredData = data(data.time == 180 & data.dex_conc == 100, :);
95+
filteredData = data(data.time == 10 & data.dex_conc == 100, :);
9696

9797
% Randomly sample 1000 points from the filtered data
9898
numSamples = min(1000, height(filteredData)); % Ensure we don’t sample more than available data
@@ -129,7 +129,6 @@
129129
GRpars_a = cell2mat(ModelGR_a.parameters(1:8,2))';
130130
GRpars_b = cell2mat(ModelGR_b.parameters(1:5,2))';
131131

132-
133132
% The log prior will be applied to the fit to multiple models as an additional constraint.
134133
log10PriorMean_a = [-2 -5 -6 -5 0.5 -3 -3 -2];
135134
log10PriorMean_b = [-2 -5 -6 -2 -5];
@@ -146,23 +145,43 @@
146145
'100','100',103,'GR Fit (100nM Dex)'};
147146

148147
ModelGRparameterMap_a = cell(1,size(GRfitCases,1));
148+
ModelGRparameterMap_b = cell(1,size(GRfitCases,1));
149149
ModelGRfit_a = cell(1,size(GRfitCases,1));
150-
ModelGRfit_b = cell(1,1);
151-
152-
% Load data for GR-beta
153-
ModelGRfit_b{1} = ModelGR_b.loadData("../EricModel/EricData/GR_ALL_gated_with_CytoArea_and_normGR_Feb2825_03.csv",...
154-
{'nucGR_b','normGRnuc';'cytGR_b','normGRcyt'},{'dex_conc','100'});
150+
ModelGRfit_b = cell(1,size(GRfitCases,1));
151+
152+
%% Load data for GR-alpha for 3 experimental conditions (Dex concentration 1,10,100nM)
153+
% NEW GR DATA
154+
ModelGRfit_a{1} = ModelGR_a.loadData("EricData/GR_ALL_gated_with_CytoArea_and_normGR_Feb2825_03_dex_conc_1.csv",...
155+
{'nucGR_a','normGRnuc';'cytGR_a','normGRcyt'});
156+
ModelGRfit_a{2} = ModelGR_a.loadData("EricData/GR_ALL_gated_with_CytoArea_and_normGR_Feb2825_03_dex_conc_10.csv",...
157+
{'nucGR_a','normGRnuc';'cytGR_a','normGRcyt'});
158+
ModelGRfit_a{3} = ModelGR_a.loadData("EricData/GR_ALL_gated_with_CytoArea_and_normGR_Feb2825_03_dex_conc_100.csv",...
159+
{'nucGR_a','normGRnuc';'cytGR_a','normGRcyt'});
160+
161+
%% Load data for GR-beta (Dex conc. does not affect GR-beta, but it's easier to compare with GR-alpha
162+
%% to compare data from each of the 3 conditions)
163+
% NEW GR DATA
164+
ModelGRfit_b{1} = ModelGR_b.loadData("EricData/GR_ALL_gated_with_CytoArea_and_normGR_Feb2825_03_dex_conc_1.csv",...
165+
{'nucGR_b','normGRnuc';'cytGR_b','normGRcyt'});
166+
ModelGRfit_b{2} = ModelGR_b.loadData("EricData/GR_ALL_gated_with_CytoArea_and_normGR_Feb2825_03_dex_conc_10.csv",...
167+
{'nucGR_b','normGRnuc';'cytGR_b','normGRcyt'});
168+
ModelGRfit_b{3} = ModelGR_b.loadData("EricData/GR_ALL_gated_with_CytoArea_and_normGR_Feb2825_03_dex_conc_100.csv",...
169+
{'nucGR_b','normGRnuc';'cytGR_b','normGRcyt'});
155170

156171
% ModelGRODEfit = cell(1,size(GRfitCases,1));
157172

158173
% Load data, fit 3 Dex concentrations to GR-alpha's 'Dex0' parameter,
159174
% and set solution scheme to FSP.
160175
for i=1:3
161-
ModelGRfit_a{i} = ModelGR_a.loadData("../EricModel/EricData/GR_ALL_gated_with_CytoArea_and_normGR_Feb2825_03.csv",...
162-
{'nucGR_a','normGRnuc';'cytGR_a','normGRcyt'},...
163-
{'dex_conc',GRfitCases{i,2}});
176+
%ModelGRfit_a{i} = ModelGR_a.loadData("../EricModel/EricData/GR_ALL_gated_with_CytoArea_and_normGR_Feb2825_03.csv",...
177+
% {'nucGR_a','normGRnuc';'cytGR_a','normGRcyt'},...
178+
% {'dex_conc',GRfitCases{i,2}});
179+
%ModelGRfit_b{i} = ModelGR_b.loadData("../EricModel/EricData/GR_ALL_gated_with_CytoArea_and_normGR_Feb2825_03.csv",...
180+
% {'nucGR_b','normGRnuc';'cytGR_b','normGRcyt'},...
181+
% {'dex_conc',GRfitCases{i,2}});
164182
ModelGRfit_a{i}.parameters(9,:) = {'Dex0', str2num(GRfitCases{i,1})};
165-
ModelGRparameterMap_a(i) = {(1:8)};
183+
ModelGRparameterMap_a(i) = {(1:8)};
184+
ModelGRparameterMap_b(i) = {(1:5)};
166185
end
167186

168187
% Make Guesses for the FSP bounds
@@ -188,16 +207,17 @@
188207
% TODO: Automate with statistics.
189208

190209
fitIters = 3; % 3 concentrations of Dex (1,10,100nM)
191-
fitMHiters = 2; % Adjust as needed
210+
fitMHiters = 1; % Adjust as needed
192211

193212
for GR = 1:fitMHiters
194213
% Specify dataset time points for GR-alpha.
195214
for i = 1:3
196215
ModelGRfit_a{i}.tSpan = ModelGRfit_a{i}.dataSet.times;
216+
ModelGRfit_b{i}.tSpan = ModelGRfit_b{i}.dataSet.times;
197217
end
198218

199219
% Specify dataset time points for GR-beta.
200-
ModelGRfit_b{1}.tSpan = ModelGRfit_b{1}.dataSet.times;
220+
%ModelGRfit_b{1}.tSpan = ModelGRfit_b{1}.dataSet.times;
201221

202222

203223
%%
@@ -267,17 +287,23 @@
267287
combinedGRModel_a = combinedGRModel_a.updateModels(GRpars_a,false);
268288
GRpars_a = combinedGRModel_a.maximizeLikelihood(GRpars_a, fitOptions);
269289
save('combinedGRModel_a','GRpars_a')
290+
%% Solve for GR-beta model.
291+
combinedGRModel_b = SSITMultiModel(ModelGRfit_b,ModelGRparameterMap_b,logPriorGR_b);
292+
combinedGRModel_b = combinedGRModel_b.initializeStateSpaces(boundGuesses);
293+
combinedGRModel_b = combinedGRModel_b.updateModels(GRpars_b,false);
294+
GRpars_b = combinedGRModel_b.maximizeLikelihood(GRpars_b, fitOptions);
295+
save('combinedGRModel_b','GRpars_b')
270296
end
271297

272298
%% Solve for GR-beta model.
273-
ModelGRfit_b{1}.fspOptions.fspTol = 1e-4;
274-
ModelGRfit_b{1}.fspOptions.bounds = boundGuesses{1};
275-
[GR_b_fspSoln,ModelGRfit_b{1}.fspOptions.bounds] = ModelGRfit_b{1}.solve;
276-
[GR_b_fspSoln,ModelGRfit_b{1}.fspOptions.bounds] = ModelGRfit_b{1}.solve(GR_b_fspSoln.stateSpace);
277-
GRpars_b = ModelGRfit_b{1}.maximizeLikelihood(GRpars_b, fitOptions);
278-
save('ModelGRfit_b','GRpars_b','GR_b_fspSoln')
279-
280-
save('GRpars_a','combinedGRModel_a', 'ModelGRfit_a','ModelGRfit_b','GRpars_b','GR_b_fspSoln')
299+
% ModelGRfit_b{1}.fspOptions.fspTol = 1e-4;
300+
% ModelGRfit_b{1}.fspOptions.bounds = boundGuesses{1};
301+
% [GR_b_fspSoln,ModelGRfit_b{1}.fspOptions.bounds] = ModelGRfit_b{1}.solve;
302+
% [GR_b_fspSoln,ModelGRfit_b{1}.fspOptions.bounds] = ModelGRfit_b{1}.solve(GR_b_fspSoln.stateSpace);
303+
% GRpars_b = ModelGRfit_b{1}.maximizeLikelihood(GRpars_b, fitOptions);
304+
% save('ModelGRfit_b','GRpars_b','GR_b_fspSoln')
305+
%
306+
% save('GRpars_a','combinedGRModel_a', 'ModelGRfit_a','ModelGRfit_b','GRpars_b','GR_b_fspSoln')
281307

282308

283309
%% Compute FIM
@@ -288,6 +314,10 @@
288314
diag(1./(log10PriorStd_a(ModelGR_a.fittingOptions.modelVarsToFit)*log(10)).^2); % Add prior in log space.
289315

290316
% for ModelGR_b parameters.
317+
combinedGRModel_b = combinedGRModel_b.computeFIMs([],'log');
318+
319+
fimGR_b_withPrior = combinedGRModel_b.FIM.totalFIM+... % the FIM in log space.
320+
diag(1./(log10PriorStd_b(ModelGR_b.fittingOptions.modelVarsToFit)*log(10)).^2); % Add prior in log space.
291321
%ModelGR_b_fimResults = ModelGRfit_b{1}.computeFIM([],'log'); % Compute individual FIMs
292322

293323
%fimGR_b_withPrior = ModelGRfit_b{1}.totalFim(ModelGR_b_fimResults,ModelGRfit_b{1}.dataSet.nCells,diag(1./(log10PriorStd_b(ModelGRfit_b{1}.fittingOptions.modelVarsToFit)*log(10)).^2)); % Add prior in log space.
@@ -300,13 +330,12 @@
300330
fimGR_a_covFree = fimGR_a_withPrior^-1;
301331
fimGR_a_covFree = 0.5*(fimGR_a_covFree+fimGR_a_covFree');
302332

303-
% if min(eig(fimGR_b_withPrior{1}))<1
304-
% disp('Warning -- FIM has one or more small eigenvalues. Reducing proposal width to 10x in those directions. MH Convergence may be slow.')
305-
% fimGR_b_withPrior{1} = fimGR_b_withPrior{1} + 1*eye(length(fimGR_b_withPrior{1}));
306-
% end
307-
% fimGR_b_covFree = fimGR_b_withPrior{1}^-1;
308-
% fimGR_b_covFree = 0.5*(fimGR_b_covFree+fimGR_b_covFree');
309-
333+
if min(eig(fimGR_b_withPrior))<1
334+
disp('Warning -- FIM has one or more small eigenvalues. Reducing proposal width to 10x in those directions. MH Convergence may be slow.')
335+
fimGR_b_withPrior = fimGR_b_withPrior + 1*eye(length(fimGR_b_withPrior));
336+
end
337+
fimGR_b_covFree = fimGR_b_withPrior^-1;
338+
fimGR_b_covFree = 0.5*(fimGR_b_covFree+fimGR_b_covFree');
310339

311340
%% Run MH on GR Models.
312341
%GRpars = GRpars';
@@ -339,20 +368,20 @@
339368
end
340369
end
341370
%% GR_beta
342-
% MHFitOptions_b.thin=1;
343-
% MHFitOptions_b.numberOfSamples=1000;
344-
% MHFitOptions_b.burnIn=100;
345-
% MHFitOptions_b.progress=true;
346-
% MHFitOptions_b.numChains = 1;
371+
% MHfitOptions_b.thin=1;
372+
% MHfitOptions_b.numberOfSamples=1000;
373+
% MHfitOptions_b.burnIn=100;
374+
% MHfitOptions_b.progress=true;
375+
% MHfitOptions_b.numChains = 1;
347376
%
348377
% % Use FIM computed above rather than making SSIT call 'useFIMforMetHast'
349378
% % which forces SSIT.m to compute it within (no prior, etc.)
350-
% MHFitOptions_b.useFIMforMetHast = false;
351-
% MHFitOptions_b.proposalDistribution=@(y)mvnrnd(y,fimGR_b_covFree);
379+
% MHfitOptions_b.useFIMforMetHast = false;
380+
% MHfitOptions_b.proposalDistribution=@(x)mvnrnd(x,fimGR_b_covFree);
352381
%
353-
% MHFitOptions_b.saveFile = 'TMPEricMHGR_b.mat';
354-
% [~,~,MHResultsGR_b] = ModelGRfit_b{1}.maximizeLikelihood(...
355-
% GRpars_b, MHFitOptions_b, 'MetropolisHastings');
382+
% MHfitOptions_b.saveFile = 'TMPEricMHGR_b.mat';
383+
% [~,~,MHResultsGR_b] = combinedGRModel_b.maximizeLikelihood(...
384+
% GRpars_b, MHfitOptions_b, 'MetropolisHastings');
356385
% %delete(MHFitOptions.saveFile)
357386
% %MHResultsGR
358387
% %
@@ -376,6 +405,7 @@
376405

377406
for model=1:size(GRfitCases,1)
378407
dataDex = data(data.dex_conc == str2double(GRfitCases{model}), :);
408+
%for t=1:max(numel(fspSolnsSMM(model).fsp),numel(GR_b_fspSoln.fsp))
379409
for t=1:max(numel(fspSolnsSMM(model).fsp),numel(GR_b_fspSoln.fsp))
380410
% Figure index
381411
figIndex = (model - 1) * 7 + t;
@@ -392,7 +422,7 @@
392422
randomIndices = randperm(height(dataDexTime), numSamples);
393423
normGRnuc = dataDexTime.normGRnuc(randomIndices);
394424
normGRcyt = dataDexTime.normGRcyt(randomIndices);
395-
425+
396426
% t is time point, so solution tensors are FSP probabilities across states for each time point
397427
conv2solnTensor_postData{t} = conv2(double(fspSolnsSMM(model).fsp{t}.p.data),double(GR_b_fspSoln.fsp{t}.p.data));
398428
figure(figIndex)
@@ -408,12 +438,12 @@
408438
contourf(log10(fspsoln_sptensor_a_postData{t}))
409439
hold on
410440
scatter(normGRcyt, normGRnuc, 'r', 'filled') % Red dots
411-
441+
412442
subplot(1,3,2)
413443
contourf(log10(fspsoln_sptensor_b_postData{t}))
414444
hold on
415445
scatter(normGRcyt, normGRnuc, 'r', 'filled') % Red dots
416-
446+
417447
subplot(1,3,3)
418448
contourf(log10(conv2solnTensor_postData{t}))
419449
hold on
@@ -426,7 +456,7 @@
426456
data = readtable('../EricModel/EricData/GR_ALL_gated_with_CytoArea_and_normGR_Feb2825_03.csv');
427457

428458
% Filter the data if desired for particular Time_index and/or dex_conc
429-
filteredData = data(data.time == 150 & data.dex_conc == 10, :);
459+
filteredData = data(data.time == 180 & data.dex_conc == 10, :);
430460

431461
% Randomly sample 1000 points from the filtered data
432462
numSamples = min(1000, height(filteredData)); % Ensure we don’t sample more than available data
@@ -440,7 +470,14 @@
440470
for logL=1:max(size(conv2solnTensor_postData))
441471
logLikelihood = sum(log(conv2solnTensor_postData{logL}(normgr>0)).*normgr((normgr>0)),"all")
442472
end
443-
%%
473+
%% STEP 1.F. -- Make Plots of GR Fit Results
474+
makeGRPlots(combinedGRModel_a,GRpars_a)
475+
476+
% ModelGRfit_b{1}.makeFitPlot([],1,[],true,'STD')
477+
478+
makeGRPlots(combinedGRModel_b,GRpars_b)
479+
480+
444481
save('conv2solnTensor_postData','GRpars_a','combinedGRModel_a', 'ModelGRfit_a','ModelGRfit_b','GRpars_b','GR_b_fspSoln', 'fspSolnsSMM')
445482

446483
dusp1 = input('(0) GR-beta turns off the DUSP1 gene;\n(1) GR-beta has no effect on DUSP1;\nChoose your destiny: ');

‎WorkSpace/EricModel/a0_Fit_GR_and_DUSP1_models.m

+5-5
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@
1010
clear
1111
addpath(genpath('../../src'));
1212
loadPrevious = false;
13-
savedWorkspace = 'workspaceDec9_2024';
13+
savedWorkspace = 'workspaceMar9_2025';
1414
addpath('tmpPropensityFunctions');
1515

1616
%% GR model
@@ -406,7 +406,7 @@
406406

407407
%% STEP 2.D.3. -- Run Metropolis Hastings Search
408408
if loadPrevious
409-
MHDusp1File = 'MHDusp1_Dec92024';
409+
MHDusp1File = 'MHDusp1_Dec92025';
410410
load(MHDusp1File)
411411
else
412412
MHFitOptions.proposalDistribution=@(x)mvnrnd(x,covFree);
@@ -422,7 +422,7 @@
422422
ModelGRDusp100nM.parameters(1:4,2) = num2cell(DUSP1pars);
423423
end
424424

425-
save('workspaceDec9_2024.mat','ModelGRDusp100nM','DUSP1pars','fimTotal','sensSoln','combinedGRModel','MHResultsDusp1')
425+
save('workspaceMar9_2025.mat','ModelGRDusp100nM','DUSP1pars','fimTotal','sensSoln','combinedGRModel','MHResultsDusp1')
426426

427427
%% STEP 2.D.4. -- Plot the MH results
428428
figNew = figure;
@@ -459,7 +459,7 @@
459459
'sensSoln'
460460
});
461461

462-
save('workspaceDec9_2024',varNames{:}) % WARNING: THIS OVERWRITE THE PREVIOUSLY SAVED WORKSPACE - TODO: FIX
462+
save('workspaceMar9_2025',varNames{:}) % WARNING: THIS OVERWRITE THE PREVIOUSLY SAVED WORKSPACE - TODO: FIX
463463

464464
%% STEP 3. -- Model Extensions using ODE Analyses
465465
if loadPrevious
@@ -723,7 +723,7 @@
723723
'sensSoln'
724724
});
725725

726-
save('workspaceDec9_2024',varNames{:})
726+
save('workspaceMar9_2025',varNames{:})
727727

728728
%% Extra Function
729729

0 commit comments

Comments
 (0)