Skip to content

Commit 20d98ed

Browse files
alexpopingaalexpopinga
committed
Fit GR-alpha, -beta, and DUSP1 obj functions simul.
1 parent 969c7f2 commit 20d98ed

File tree

1 file changed

+72
-45
lines changed

1 file changed

+72
-45
lines changed

‎WorkSpace/EricModel/GR_2_Electric_Boogaloo.m

+72-45
Original file line numberDiff line numberDiff line change
@@ -411,7 +411,7 @@
411411
%ModelGRDusp100nM.fspOptions.bounds = [0;0;0;2;2;400];
412412
ModelGRDusp100nM.fittingOptions.modelVarsToFit = 10:13;
413413
ModelGRDusp100nM = ModelGRDusp100nM.formPropensitiesGeneral('EricModDusp1_1');
414-
log10PriorMean = [-2 -5 -6 -5 0.5 -3 -3 -2 -1 -1 0 -2 -2 -5];
414+
log10PriorMean = [-2 -5 -6 -5 0.5 -3 -3 -2 2 -1 -1 0 -2 -2 -5];
415415
log10PriorStd = 2*ones(1,15);
416416
duspLogPrior = @(x)-sum((log10(x(:))'-log10PriorMean(9:12)).^2./(2*log10PriorStd(9:12).^2));
417417
ModelGRDusp100nM.fittingOptions.logPrior = duspLogPrior;
@@ -429,7 +429,7 @@
429429
ModelDusp1Fit{i}.inputExpressions = {'IDex','Dex0*exp(-gDex*t)'};
430430
ModelDusp1parameterMap{i} = (1:4);
431431
% Set Dex concentration.
432-
ModelDusp1Fit{i}.parameters{13,2} = str2num(Dusp1FitCases{i,1});
432+
ModelDusp1Fit{i}.parameters{9,2} = str2num(Dusp1FitCases{i,1});
433433
ModelDusp1Fit{i} = ModelDusp1Fit{i}.formPropensitiesGeneral(['EricModDusp1_',num2str(i),'_FSP']);
434434
end
435435
DUSP1pars = [ModelDusp1Fit{i}.parameters{ModelGRDusp100nM.fittingOptions.modelVarsToFit,2}];
@@ -453,26 +453,26 @@
453453
save('EricModelDusp1_MMDex','GRpars_a','GRpars_b','DUSP1pars')
454454
end
455455

456-
%% STEP 2.C. -- Plot predictions for other Dex concentrations.
457-
showCases = [1,1,1,1];
458-
makePlotsDUSP1({ModelGRDusp100nM},ModelGRDusp100nM,DUSP1pars,Dusp1FitCases,showCases)
456+
%% Plot predictions for other Dex concentrations.
457+
%showCases = [1,1,1,1];
458+
%makePlotsDUSP1({ModelGRDusp100nM},ModelGRDusp100nM,DUSP1pars,Dusp1FitCases,showCases)
459459

460-
%% STEP 2.D. -- Sample uncertainty for Dusp1 Parameters
461-
% STEP 2.D.1. -- Compute sensitivity of the FSP solution
460+
%% Sample uncertainty for Dusp1 Parameters
461+
% Compute sensitivity of the FSP solution
462462
ModelGRDusp100nM.solutionScheme = 'fspSens';
463463
sensSoln = ModelGRDusp100nM.solve();
464464
ModelGRDusp100nM.solutionScheme = 'FSP';
465-
% STEP 2.D.2. -- Compute FIM
465+
% Compute FIM
466466
% define which species in model are not observed.
467467
ModelGRDusp100nM.pdoOptions.unobservedSpecies = {'offGene';'onGene'};
468468
% compute the FIM
469469
fimResults = ModelGRDusp100nM.computeFIM(sensSoln.sens,'log');
470470
% In the following, the log-prior is used as a prior co-variance matrix.
471471
% This will be used in the FIM calculation as an FIM without new evidence
472472
% being set equal to the inverse of this covariance matrix. More rigorous
473-
% justification is needed to support this heuristic.
473+
%% justification is needed to support this heuristic.
474474
fimTotal = ModelGRDusp100nM.evaluateExperiment(fimResults,ModelGRDusp100nM.dataSet.nCells,...
475-
diag(log10PriorStd(1:13).^2));
475+
diag(log10PriorStd(1:15).^2));
476476
FIMfree = fimTotal{1}(1:4,1:4);
477477
if min(eig(FIMfree))<1
478478
disp('Warning -- FIM has one or more small eigenvalues. Reducing proposal width to 10x in those directions. MH Convergence may be slow.')
@@ -499,7 +499,7 @@
499499
ModelGRDusp100nM.parameters(1:4,2) = num2cell(DUSP1pars);
500500
end
501501

502-
save('workspace_Feb4_2024.mat','ModelGRDusp100nM','DUSP1pars','fimTotal','sensSoln','combinedGRModel','MHResultsDusp1')
502+
save('workspace_Feb4_2024.mat','ModelGRDusp100nM','DUSP1pars','fimTotal','sensSoln','combinedGRModel_a','MHResultsDusp1')
503503

504504
%% Plot the MH results
505505
figNew = figure;
@@ -514,15 +514,11 @@
514514
end
515515
end
516516

517-
%
518-
% %% Switch solver to ODE and generate model codes
519-
% ModelGRDusp100nM_a.solutionScheme = 'fsp';
520-
% ModelGRDusp100nM_b.solutionScheme = 'fsp';
521-
% ModelGRDusp100nM_a.useHybrid = false;
522-
% ModelGRDusp100nM_b.useHybrid = false;
523-
% ModelGRDusp100nM_a = ModelGRDusp100nM_a.formPropensitiesGeneral('fspDUSP1_a');
524-
% ModelGRDusp100nM_a = ModelGRDusp100nM_a.formPropensitiesGeneral('fspDUSP1_b');
525-
% %% Change parameters manually, solve, and make plots
517+
%% Switch solver to ODE
518+
ModelGRDusp100nM.solutionScheme = 'fsp';
519+
ModelGRDusp100nM.useHybrid = false;
520+
ModelGRDusp100nM = ModelGRDusp100nM.formPropensitiesGeneral('fspDUSP1_a');
521+
%% Change parameters manually, solve, and make plots
526522
% [Dusp_soln_a,ModelGRDusp100nM_a.fspOptions.bounds] = ModelGRDusp100nM_a.solve;
527523
% [Dusp_soln_b,ModelGRDusp100nM_b.fspOptions.bounds] = ModelGRDusp100nM_b.solve;
528524
% %plotODEresults(ModelGRDusp100nM_a,Dusp_soln_a,ModelGRfit_a{3})
@@ -619,27 +615,27 @@
619615
save('EricModelDusp1_MMDex','GRpars_a','GRpars_b','DUSP1pars')
620616
end
621617

622-
%% STEP 2.C. -- Plot predictions for other Dex concentrations.
623-
showCases = [1,1,1,1];
624-
makePlotsDUSP1({ModelGRDusp100nM},ModelGRDusp100nM,DUSP1pars,Dusp1FitCases,showCases)
618+
%% Plot predictions for other Dex concentrations.
619+
%showCases = [1,1,1,1];
620+
%makePlotsDUSP1({ModelGRDusp100nM},ModelGRDusp100nM,DUSP1pars,Dusp1FitCases,showCases)
625621

626-
%% STEP 2.D. -- Sample uncertainty for Dusp1 Parameters
627-
% STEP 2.D.1. -- Compute sensitivity of the FSP solution
622+
%% Sample uncertainty for Dusp1 Parameters
623+
% Compute sensitivity of the FSP solution
628624
ModelGRDusp100nM.solutionScheme = 'fspSens';
629625
sensSoln = ModelGRDusp100nM.solve();
630626
ModelGRDusp100nM.solutionScheme = 'FSP';
631-
% STEP 2.D.2. -- Compute FIM
627+
% Compute FIM
632628
% define which species in model are not observed.
633629
ModelGRDusp100nM.pdoOptions.unobservedSpecies = {'offGene';'onGene'};
634630
% compute the FIM
635631
fimResults = ModelGRDusp100nM.computeFIM(sensSoln.sens,'log');
636632
% In the following, the log-prior is used as a prior co-variance matrix.
637633
% This will be used in the FIM calculation as an FIM without new evidence
638634
% being set equal to the inverse of this covariance matrix. More rigorous
639-
% justification is needed to support this heuristic.
635+
%% justification is needed to support this heuristic.
640636
fimTotal = ModelGRDusp100nM.evaluateExperiment(fimResults,ModelGRDusp100nM.dataSet.nCells,...
641-
diag(log10PriorStd(1:13).^2));
642-
FIMfree = fimTotal{1}(1:4,1:4);
637+
diag(log10PriorStd(1:12).^2));
638+
FIMfree = fimTotal{1}(10:13,10:13);
643639
if min(eig(FIMfree))<1
644640
disp('Warning -- FIM has one or more small eigenvalues. Reducing proposal width to 10x in those directions. MH Convergence may be slow.')
645641
FIMfree = FIMfree + 1*eye(length(FIMfree));
@@ -649,7 +645,7 @@
649645

650646
%% STEP 2.D.3. -- Run Metropolis Hastings Search
651647
if loadPrevious
652-
MHDusp1File = 'MHDusp1_Dec92024';
648+
MHDusp1File = 'MHDusp1_Feb_2025';
653649
load(MHDusp1File)
654650
else
655651
MHFitOptions.proposalDistribution=@(x)mvnrnd(x,covFree);
@@ -662,12 +658,12 @@
662658
[DUSP1pars,~,MHResultsDusp1] = ModelGRDusp100nM.maximizeLikelihood(...
663659
[], MHFitOptions, 'MetropolisHastings');
664660
delete('TMPEricMHDusp1.mat')
665-
ModelGRDusp100nM.parameters(1:4,2) = num2cell(DUSP1pars);
661+
ModelGRDusp100nM.parameters(10:13,2) = num2cell(DUSP1pars);
666662
end
667663

668664
save('workspace_Feb4_2024.mat','ModelGRDusp100nM','DUSP1pars','fimTotal','sensSoln','combinedGRModel','MHResultsDusp1')
669665

670-
%% STEP 2.D.4. -- Plot the MH results
666+
%% Plot the MH results
671667
figNew = figure;
672668
ModelGRDusp100nM.plotMHResults(MHResultsDusp1,[],'log',[],figNew)
673669
for j = 1:3
@@ -681,41 +677,72 @@
681677
end
682678

683679
%% Save results
684-
varNames = unique({'ModelGR'
680+
varNames = unique({'ModelGR_a'
681+
'ModelGR_b'
685682
'GRfitCases'
686683
'log10PriorMean'
687684
'log10PriorStd'
688-
'GRpars'
689-
'MHResultsGR'
685+
'GRpars_a'
686+
'GRpars_b'
690687
'ModelGRparameterMap'
691-
'ModelGRfit'
692-
'boundGuesses'
688+
'combinedGR_a'
693689
'ModelGRDusp100nM'
694-
'GRfitCases'
695-
'log10PriorMean'
696-
'log10PriorStd'
697690
'duspLogPrior'
698691
'DUSP1pars'
699-
'ModelGRfit'
692+
'ModelGRfit_a'
693+
'ModelGRfit_b'
700694
'fimResults'
701695
'MHResultsDusp1'
702-
'sensSoln'
696+
'ModelGRDusp100nM'
703697
});
704698

705-
save('workspaceDec9_2024',varNames{:}) % WARNING: THIS OVERWRITE THE PREVIOUSLY SAVED WORKSPACE - TODO: FIX
699+
save('workspace_Feb_2025',varNames{:}) % WARNING: THIS OVERWRITE THE PREVIOUSLY SAVED WORKSPACE - TODO: FIX
706700
end
707701

708702
save('conv2solnTensor_postData','GRpars_a','combinedGRModel_a', 'ModelGRfit_a','ModelGRfit_b','GRpars_b','GR_b_fspSoln', 'fspSolnsSMM',...
709703
'fimGR_a_withPrior','fimGR_b_withPrior','ModelGR_b_fimResults','fimGR_a_covFree','fimGR_b_covFree','MHResultsGR_a','MHFitOptions','ModelGRfit');
710-
704+
705+
%% Create new objective function combining all of the previous ones.
706+
% Remove all priors from individual models.
707+
for i=1:3
708+
ModelGRfit_a{i}.fittingOptions.logPrior = [];
709+
ModelGRfit_a{i}.tSpan = ModelGRfit_a{i}.dataSet.times;
710+
end
711+
ModelGRfit_b{1}.fittingOptions.logPrior = [];
712+
ModelGRfit_b{1}.tSpan = ModelGRfit_b{1}.dataSet.times;
713+
ModelGRDusp100nM.fittingOptions.logPrior = [];
714+
fullPars = [ModelGRDusp100nM.parameters{1:15,2}];
715+
716+
%% Fit all objective functions at once.
717+
log10PriorMean = [-2 -5 -6 ... % GR pars
718+
-5 0.5 -3 -3 -2 ... % GR-alpha pars
719+
2 ... % Dex
720+
-1 -1 0 -2 ... % Dusp1 pars
721+
-2 -5]; % GR-beta pars
722+
log10PriorStd = 2*ones(1,15);
723+
%duspLogPrior = @(x)-sum((log10(x(:))'-log10PriorMean(9:12)).^2./(2*log10PriorStd(9:12).^2));
724+
logPriorAll = @(x)-sum((log10(x)-log10PriorMean(1:15)).^2./(2*log10PriorStd(1:15).^2));
725+
726+
% extendedMod.fittingOptions.modelVarsToFit = [1:12,14,15];
727+
Organization = {ModelGRfit_a{3},[1:8],[1:8],'computeLikelihood',1;...
728+
ModelGRfit_b{1},[1:5],[1:5],'computeLikelihood',1;...
729+
ModelGRDusp100nM,[1:8,10:15],[1:8,10:15],'computeLikelihood',1};
730+
%extendedMod,[1:12,14:15],[1:14],'computeLikelihoodODE',0.01};
731+
Organization = getTotalFitErr(Organization,fullPars,true);
732+
getTotalFitErr(Organization,fullPars,false)
733+
objAll = @(x)-getTotalFitErr(Organization,exp(x),false)-logPriorAll(exp(x));
734+
for jj = 1:fitIters
735+
fullParsLog = log(fullPars);
736+
fullPars = exp(fminsearch(objAll,fullParsLog,fitOptions));
737+
end
711738

712739
%% compute likelihood
713740
% logLikelihood = sum(log(conv2solnTensor{f}(PAs>0)).*PAs((PAs>0)),"all");
714741
%logLikelihood = sum(log(conv2solnTensor_postData{t}(fspsoln_sptensor_a_postData{t}>0)).*fspsoln_sptensor_a_postData{t}((fspsoln_sptensor_a_postData{t}>0)),"all");
715742

716743
case 2
717744

718-
%% STEP 3.A.1 -- Extend model to include nuclear and cytoplasmic RNA
745+
%% Extend model to include nuclear and cytoplasmic RNA
719746
extendedMod = ModelGRDusp100nM;
720747
extendedMod = extendedMod.addSpecies({'rCyt'},0);
721748
% Adjust the final reaction to be nuc -> cyt transport.

0 commit comments

Comments
 (0)