Add files via upload

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 [#/µL]' ;
+Output{22} = 'Sa = sum(S1:S8) [#/µL]' ;
+Output{23} = 'Na = N1+N2 [#/µL]' ;
+Output{24} = 'S = sum(S0:S8) [#/µL]' ;
+Output{25} = 'N = N0+N1+N2 [#/µL]' ;
+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(' ') ;

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 ;