diff --git a/N803_collector_2.m b/N803_collector_2.m new file mode 100644 index 0000000..629a510 --- /dev/null +++ b/N803_collector_2.m @@ -0,0 +1,231 @@ +% collector.m - collects outputs from model calibration and UQ +addpath Data + +% This script will load the MCMC parameter samples, thin them by selecting +% every 100th sample, solve the ODE model, and store the results. + +%% ------------------------------------------------------------------------ +% INPUTS +% ------------------------------------------------------------------------ + +name = 'N803_MCMC' ;% file name for storing output + +nWork = 4 ;% number of parallel workers +nIter = 1e6 ;% total number of MCMC iteration +nThin = 100 ;% thin sample by every 'nThin' samples + +% scenario descriptions +Scenario{1} = 'Fitting to both (C1)' ;% Simultaneous fitting - Cohort 1 +Scenario{2} = 'Fitting to both (C2)' ;% Simultaneous fitting - Cohort 2 +Scenario{3} = 'Fitting to Cohort 1 only' ;% Seperate fitting - Cohort 1 +Scenario{4} = 'Fitting to Cohort 2 only' ;% Seperate fitting - Cohort 2 +Scenario{5} = 'Fitting no N to both (C1)' ;% Model w/o non-SIV-spec CD8 T +Scenario{6} = 'Fitting no N to both (C2)' ;% Model w/o non-SIV-spec CD8 T +Scenario{7} = 'Webb18 Validation (naive)' ;% Validation of shared fitting +Scenario{8} = 'Webb18 Validation (SIV)' ;% Validation of shared fitting + +% x-values (days) corresponding to outputs +X{1} = [ (- 4:0)*7 (.1:.1:6.9) (7:.5:44*7) ] ; +X{2} = [ (-22:0)*7 (.1:.1:6.9) (7:.5:10*7) ] ; +X{3} = [ (- 4:0)*7 (1:.5:44*7) ] ; +X{4} = [ (-22:0)*7 (1:.5:10*7) ] ; +X{5} = [ (- 4:0)*7 (1:.5:44*7) ] ; +X{6} = [ (-22:0)*7 (1:.5:10*7) ] ; +X{7} = -7:.1: 4*7 ; +X{8} = -7:.2:18*7 ; + +% x-values (days) for doses given +Dose{1} = 7*[0:3 5:8 37:40] ; +Dose{2} = 7*[0 2 4] ; +Dose{3} = 7*[0:3 5:8 37:40] ; +Dose{4} = 7*[0 2 4] ; +Dose{5} = 7*[0:3 5:8 37:40] ; +Dose{6} = 7*[0 2 4] ; +Dose{7} = 0 ; +Dose{8} = 7*[0 1 2 7] ; + +% Descriptions of additional model outputs (see 'N803_model_2.m') +Output{01} = 'virions [log fold change]' ; +Output{02} = 'CD8+ T cells [fold change]' ; +Output{03} = 'V [log wrt C1 pre]' ; +Output{04} = 'S_0 [log wrt C1 pre]' ; +Output{05} = 'S_1 [log wrt C1 pre]' ; +Output{06} = 'S_2 [log wrt C1 pre]' ; +Output{07} = 'S_3 [log wrt C1 pre]' ; +Output{08} = 'S_4 [log wrt C1 pre]' ; +Output{09} = 'S_5 [log wrt C1 pre]' ; +Output{10} = 'S_6 [log wrt C1 pre]' ; +Output{11} = 'S_7 [log wrt C1 pre]' ; +Output{12} = 'S_8 [log wrt C1 pre]' ; +Output{13} = 'N_0 [log wrt C1 pre]' ; +Output{14} = 'N_1 [log wrt C1 pre]' ; +Output{15} = 'N_2 [log wrt C1 pre]' ; +Output{16} = 'X [pmol/kg]' ; +Output{17} = 'C [pM]' ; +Output{18} = 'R_1 [log wrt C1 pre]' ; +Output{19} = 'R_2 [log wrt C1 pre]' ; +Output{20} = 'C/(C+C_50)' ; +Output{21} = 'CD8+ T cells [#/無]' ; +Output{22} = 'Sa = sum(S1:S8) [#/無]' ; +Output{23} = 'Na = N1+N2 [#/無]' ; +Output{24} = 'S = sum(S0:S8) [#/無]' ; +Output{25} = 'N = N0+N1+N2 [#/無]' ; +Output{26} = 'N/(S+N) frequency' ; +Output{27} = 'Sa/S frequency' ; +Output{28} = 'Na/N frequency' ; +Output{29} = 'Sa [log wrt C1 pre]' ; +Output{30} = 'Na [log wrt C1 pre]' ; +Output{31} = 'S0 proliferation [log wrt C1 pre]' ; +Output{32} = 'N0 proliferation [log wrt C1 pre]' ; +Output{33} = 'S0 death [log wrt C1 pre]' ; +Output{34} = 'N0 death [log wrt C1 pre]' ; +Output{35} = 'S0,N0 prolif. rate (p*) [log wrt C1 pre]' ; +Output{36} = 'C multiplier for p* [log wrt C1 pre]' ; +Output{37} = 'R multiplier for p* [log wrt C1 pre]' ; +Output{38} = 'density dep for p* [log wrt C1 pre]' ; +Output{39} = 'S0 activation [log wrt C1 pre]' ; +Output{40} = 'N0 activation [log wrt C1 pre]' ; +Output{41} = 'Sa proliferation [log wrt C1 pre]' ; +Output{42} = 'Na proliferation [log wrt C1 pre]' ; +Output{43} = 'Sa death [log wrt C1 pre]' ; +Output{44} = 'Na death [log wrt C1 pre]' ; +Output{45} = 'Sa reversion [log wrt C1 pre]' ; +Output{46} = 'Na reversion [log wrt C1 pre]' ; +Output{47} = 'S0 activ. rate (aS*) [log wrt C1 pre]' ; +Output{48} = 'N0 activ. rate (aN*) [log wrt C1 pre]' ; +Output{49} = 'C multiplier for aS* [log wrt C1 pre]' ; +Output{50} = 'C multiplier for aN* [log wrt C1 pre]' ; +Output{51} = 'R multiplier for aS* [log wrt C1 pre]' ; +Output{52} = 'R multiplier for aN* [log wrt C1 pre]' ; +Output{53} = 'V multiplier for aS* [log wrt C1 pre]' ; +Output{54} = 'V multiplier for aN* [log wrt C1 pre]' ; +Output{55} = 'V growth [log wrt C1 pre]' ; +Output{56} = 'S killing [log wrt C1 pre]' ; +Output{57} = 'N killing [log wrt C1 pre]' ; +Output{58} = 'S killing rate (gS*) [log wrt C1 pre]' ; +Output{59} = 'N killing rate (gN*) [log wrt C1 pre]' ; +Output{60} = 'R multiplier for gS* [log wrt C1 pre]' ; +Output{61} = 'R multiplier for gN* [log wrt C1 pre]' ; +Output{62} = 'N killing fraction' ; +Output{63} = 'R generation by aS [log wrt C1 pre]' ; +Output{64} = 'R generation by aN [log wrt C1 pre]' ; +Output{65} = 'R gen by aN fraction' ; + +%% ------------------------------------------------------------------------ +% BODY +% ------------------------------------------------------------------------ + +% Percolate results structure +N = length(Dose) ; +C = struct('Scenario',Scenario) ; +[C(:).X] = deal(X{:}) ;% x-values (days) corresponding to outputs +[C(:).Dose] = deal(Dose{:}) ;% x-values (days) for doses given +[C(:).YSummary] = deal([]) ;% output summary statistics +[C(:).ParSample] = deal([]) ;% parameter samples + +% load(name,'C') + +disp('Collecting...') ; + +% For each model... +for n = 1:8 + + disp(Scenario{n}) ; + % Declare model function handle and results file name + if (n==1) || (n==2) + saveFile = 'N803_shared_Results' ; + Doses = { Dose{1} , Dose{2} } ; + modelFun = @(P) N803_shared(X{n}, Doses, P, {[],[]}, n,... + 'fullOut',true) ; + elseif (n==3) || (n==4) + saveFile = ['N803_cohort_' num2str(n-2) '_Results'] ; + modelFun = @(P) N803_single(X{n}, Dose{n}, P, []) ; + elseif (n==5) || (n==6) + saveFile = 'N803_shared_noN_Results' ; + Doses = { Dose{5} , Dose{6} } ; + modelFun = @(P) N803_shared_noN(X{n}, Doses, P, {[],[]}, n-4) ; + elseif (n==7) + saveFile = 'N803_shared_Results' ; + modelFun = @(P) N803_shared(X{n}, {Dose{n},[]}, P, {[],[]}, 1,... + 'fullOut',true,'ivDosing',880*[.06 0]/.03,'isNaive',true) ; + elseif (n==8) + saveFile = 'N803_shared_Results' ; + modelFun = @(P) N803_shared(X{n}, {Dose{n},[]}, P, {[],[]}, 1,... + 'fullOut',true,'ivDosing',880*[.06 .06 .06 1 0]/.03) ; + end + + % Load model parameters, etc., from MSLS and MCMC data files + warning off + load(saveFile,'Results') + load([saveFile '_MCMC'],'res') + warning on + SamID = nThin:nThin:nIter ;% thin sample index (within MCMC) + SamID(end) = nIter-1 ;% using 2nd-to-last because of an MCMC glitch !! + nSam = length(SamID) ;% thinned number of MCMC samples + PARS = repmat(Results(1).ParameterFit,nSam,1) ;% sample of parameters + ParID = Results(1).setup.SampleIndex ;% index within par array + PARS(:,ParID) = 10.^res.par(:,SamID)' ;% sample of parameters + + % Evaluate model and collect terms + Ysam = cell(nSam,1) ;% sample of model outputs + Pcoh = cell(nSam,1) ;% sample of cohort parameters + passID = true(nSam,1) ;% index of models that did not fail + parfor (sam = 1:nSam, nWork) + try + [Ysam{sam},Pcoh{sam}] = modelFun(PARS(sam,:)) ; + if any( Pcoh{sam} < 0 ) + error('Negative parameters or initial values.') + end + catch + passID(sam) = false ; + end + end + nSam = sum(passID) ; + Ysam = Ysam(passID) ; + Pcoh = Pcoh(passID) ; + PCOH = cat(1,Pcoh{:}) ;% sample of calculated parameters + + C(n).NLL = - res.logPost( SamID(passID) ) ; + + % Store confidence regions (OR samples) + OUTPUTS = zeros(length(X{n}),nSam) ;% storing output sample (temp) + nOuts = size(Ysam{1},2) ;% number of outputs + Ystats = struct('Output',Output) ; + [Ystats(:).min] = deal([]) ;% minimum values across samples + [Ystats(:).lo95] = deal([]) ;% 2.5th percentile + [Ystats(:).med] = deal([]) ;% median + [Ystats(:).hi95] = deal([]) ;% 97.5th percentile + [Ystats(:).max] = deal([]) ;% maximum + + if n == 1 + NORMFACTS = ones(nOuts,nSam) ;% matrix of normalization factors + end + for out = 1:nOuts + for sam = 1:nSam + OUTPUTS(:,sam) = Ysam{sam}(:,out) ; + end + if any( n == 1:2 ) && any( out == [ 3:15 , 18:19 , 29:61 , 63:64 ] ) + if n == 1 + NORMFACTS(out,:) = OUTPUTS(1,:) ; + end + OUTPUTS = OUTPUTS - NORMFACTS(out,:) ; + end + Ystats(out).min = prctile(OUTPUTS', 0 )' ; + Ystats(out).lo95 = prctile(OUTPUTS', 2.5)' ; + Ystats(out).med = prctile(OUTPUTS', 50 )' ; + Ystats(out).hi95 = prctile(OUTPUTS', 97.5)' ; + Ystats(out).max = prctile(OUTPUTS',100 )' ; + end + + C(n).YSummary = Ystats ;% output summary statistics + C(n).ParSample = PCOH ;% sample of calculated parameters +end + +%% ------------------------------------------------------------------------ +% OUTPUTS +% ------------------------------------------------------------------------ + +save([name '.mat'],'C') ; + +disp('Done.') ; +disp(' ') ; diff --git a/main_shared_noN_calibration_2022_11_02.m b/main_shared_noN_calibration_2022_11_02.m index c84a889..9453dcc 100644 --- a/main_shared_noN_calibration_2022_11_02.m +++ b/main_shared_noN_calibration_2022_11_02.m @@ -87,7 +87,7 @@ % costFunction(GuessBOUNDS(1,:)) ;% uncomment to debug your costFunction % filePrefix = string to prefix to file names ----------------------------- -filePrefix = 'N803_shared_2S' ; +filePrefix = 'N803_shared_noN' ; % timeLimit = number of minutes allowed for a single optimization run ----- timeLimit = 30 ;