diff --git a/N803_collector_2.m b/N803_collector_2.m index 69def17..279dec4 100644 --- a/N803_collector_2.m +++ b/N803_collector_2.m @@ -1,245 +1,286 @@ % 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. +% This script will load the MCMC parameter samples or MSLS parameter sets, +% solve the ODE model, and store the results. %% ======================================================================== % INPUTS % ======================================================================== -name = 'N803_MCMC' ;% file name for storing output +name = 'N803_ForPlots' ;% 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 +nWork = 4 ;% number of parallel workers +nMSLS = 10 ;% if storing MSLS sets, store top 'nMSLS' +nThin = 10 ;% if storing entire MCMC sample, thin to every 'nThin' % 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} = 'Webb18 Validation (naive)' ;% Validation of shared fitting -Scenario{6} = 'Webb18 Validation (SIV)' ;% Validation of shared fitting +Scenario{ 1} = 'Shared constants (Cohort1) MCMC' ; +Scenario{ 2} = 'Shared constants (Cohort2) MCMC' ; +Scenario{ 3} = 'Webb18 Comparison (naive)' ; +Scenario{ 4} = 'Webb18 Comparison (SIV)' ; +Scenario{ 5} = 'Shared constants (Cohort1) MSLS' ; +Scenario{ 6} = 'Shared constants (Cohort2) MSLS' ; +Scenario{ 7} = 'Seperate constants (Cohort1) MSLS' ; +Scenario{ 8} = 'Seperate constants (Cohort2) MSLS' ; +Scenario{ 9} = 'Model A (Cohort1) MSLS' ; +Scenario{10} = 'Model A (Cohort2) MSLS' ; +Scenario{11} = 'Model B (Cohort1) MSLS' ; +Scenario{12} = 'Model B (Cohort2) MSLS' ; +Scenario{13} = 'Model C (Cohort1) MSLS' ; +Scenario{14} = 'Model C (Cohort2) MSLS' ; +Scenario{15} = 'Model D (Cohort1) MSLS' ; +Scenario{16} = 'Model D (Cohort2) MSLS' ; +Scenario{17} = 'Long term behavior (Cohort1)' ; +Scenario{18} = 'Long term behavior (Cohort2)' ; +Scenario{19} = 'Three 2-week doses (Cohort2)' ; +Scenario{20} = 'Three 3-week doses (Cohort2)' ; +Scenario{21} = 'Three 4-week doses (Cohort2)' ; +Scenario{22} = '2-week dosing for 48w (Cohort2)' ; +Scenario{23} = '3-week dosing for 48w (Cohort2)' ; +Scenario{24} = '4-week dosing for 48w (Cohort2)' ; +Scenario{25} = 'Three 2-week doses (fold change C2)' ; +Scenario{26} = 'Three 3-week doses (fold change C2)' ; % 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} = -7:.1: 4*7 ; -X{6} = -7:.2:18*7 ; +X{3} = -7:.1: 4*7 ; +X{4} = -7:.2:18*7 ; +X(5:2:15) = { [ (- 4:0)*7 (1:.5:44*7) ] } ; +X(6:2:16) = { [ (-22:0)*7 (1:.5:10*7) ] } ; +X{17} = [ (-.5:.5:44) 44+(1:250) ] *7 ; +X{18} = [ (-.5:.5:10) 10+(1:250) ] *7 ; +X(19:21) = { -7:1:7*24 } ; +X(22:24) = { -7:1:7*52 } ; +X(25:26) = { -7:1:7*24 } ; % 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} = 0 ; -Dose{6} = 7*[0 1 2 7] ; - -% Save the parameter sample or initial condition sample -SavePars{1} = true ; SaveYics{1} = true ; -SavePars{2} = false ; SaveYics{2} = true ; -SavePars{3} = false ; SaveYics{3} = false ; -SavePars{4} = false ; SaveYics{4} = false ; -SavePars{5} = false ; SaveYics{5} = false ; -SavePars{6} = false ; SaveYics{6} = false ; - -% Save the following outputs (see list below and in 'N803_model_2.m') -SaveOutput{1} = [1 2 4 27 29 39 49 51 53 58 60] ; -SaveOutput{2} = [1 2 4 27 29 39 49 51 53 58 60] ; -SaveOutput{3} = [1 2] ; -SaveOutput{4} = [1 2] ; -SaveOutput{5} = 21 ; -SaveOutput{6} = [3 21] ; - -% 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' ; - -% Add figure labels to the following outputs (for Scenarios 1-6) -Figure = cell(size(Output)) ; -Figure{01} = {'2A','2D','S3A','S3D', [], []} ; -Figure{02} = {'2B','2E','S3B','S3E', [], []} ; -Figure{03} = { [], [], [], [], [],'S4C'} ; -Figure{04} = {'4D','4D', [], [], [], []} ; -Figure{21} = { [], [], [], [],'S4A','S4D'} ; -Figure{27} = {'3A','3C', [], [], [], []} ; -Figure{29} = {'4B','4B', [], [], [], []} ; -Figure{39} = {'4F','4F', [], [], [], []} ; -Figure{49} = {'4I','4I', [], [], [], []} ; -Figure{51} = {'4H','4H', [], [], [], []} ; -Figure{53} = {'4G','4G', [], [], [], []} ; -Figure{58} = {'4A','4A', [], [], [], []} ; -Figure{60} = {'4C','4C', [], [], [], []} ; +Dose{3} = 0 ; +Dose{4} = 7*[0 1 2 7] ; +Dose(5:2:17) = { 7*[0:3 5:8 37:40] } ; +Dose(6:2:18) = { 7*[0 2 4] } ; +Dose{19} = (0:2:4 )*7 ; +Dose{20} = (0:3:6 )*7 ; +Dose{21} = (0:4:8 )*7 ; +Dose{22} = (0:2:48)*7 ; +Dose{23} = (0:3:48)*7 ; +Dose{24} = (0:4:48)*7 ; +Dose{25} = (0:2:4 )*7 ; +Dose{26} = (0:3:6 )*7 ; + +% Save the parameter sample +SavePars{1} = true ; +SavePars{2} = true ; +SavePars(3:26) = { false } ; + +% Use MCMC sample (or MSLS sets) Save entire sample (or only 95% CI) +UseMCMC{1} = true ; SaveSample{1} = false ; +UseMCMC{2} = true ; SaveSample{2} = false ; +UseMCMC{3} = true ; SaveSample{3} = false ; +UseMCMC{4} = true ; SaveSample{4} = false ; +UseMCMC(5:16) = { false } ; SaveSample(5:26) = { true } ; +UseMCMC(17:26) = { true } ; + +% Figure label, output ID (for model function), and output description +% LFC1 = log fold change wrt cohort 1 pre-treatment value +F(1).Fig = '3A' ; F(1).Out = 1 ; F(1).Desc = 'Virions [log CEQ/mL]' ; +F(2).Fig = '3B' ; F(2).Out = 2 ; F(2).Desc = 'CD8+ T cells [#/uL]' ; +F(3).Fig = '4A' ; F(3).Out = 26 ; F(3).Desc = 'Active T8 freq. in S' ; +F(4).Fig = '5A' ; F(4).Out = 55 ; F(4).Desc = 'Viral decay via S [LFC1]' ; +F(5).Fig = '5B' ; F(5).Out = 23 ; F(5).Desc = 'Active S [LFC1]' ; +F(6).Fig = '5C' ; F(6).Out = 56 ; F(6).Desc = 'Regulation of killing [LFC1]' ; +F(7).Fig = '5D' ; F(7).Out = 4 ; F(7).Desc = 'Resting S [LFC1]' ; +F(8).Fig = '5F' ; F(8).Out = 38 ; F(8).Desc = 'S activation [LFC1]' ; +F(9).Fig = '5G' ; F(9).Out = 51 ; F(9).Desc = 'Viral term (S act.) [LFC1]' ; +F(10).Fig= '5H' ; F(10).Out= 49 ; F(10).Desc= 'Regulation of S act. [LFC1]' ; +F(11).Fig= '5I' ; F(11).Out= 47 ; F(11).Desc= 'N803 term (S act.) [LFC1]' ; +FigData{1} = F ; + +F(1).Fig = '3D' ; +F(2).Fig = '3E' ; +F(3).Fig = '4C' ; +FigData{2} = F ; clear F + +F(1).Fig = 'S9A' ; F(1).Out = 1 ; F(1).Desc = 'CD8+ T cells [#/uL]' ; +FigData{3} = F ; +F(1).Fig = 'S9C' ; F(1).Out = 1 ; F(1).Desc = 'Virions [log CEQ/mL]' ; +F(2).Fig = 'S9D' ; F(2).Out = 2 ; F(2).Desc = 'CD8+ T cells [#/uL]' ; +FigData{4} = F ; + +F(1).Fig = 'S#A' ; F(1).Out = 1 ; F(1).Desc = 'Virions [log CEQ/mL]' ; +F(2).Fig = 'S#B' ; F(2).Out = 2 ; F(2).Desc = 'CD8+ T cells [#/uL]' ; +FigData{5} = F ; +F(1).Fig = 'S#D' ; +F(2).Fig = 'S#E' ; +FigData{6} = F ; + +F(1).Fig = 'S8A' ; F(2).Fig = 'S8B' ; FigData{07} = F ; +F(1).Fig = 'S8D' ; F(2).Fig = 'S8E' ; FigData{08} = F ; +F(1).Fig = 'S2A' ; F(2).Fig = 'S2B' ; FigData{09} = F ; +F(1).Fig = 'S2D' ; F(2).Fig = 'S2E' ; FigData{10} = F ; +F(1).Fig = 'S3A' ; F(2).Fig = 'S3B' ; FigData{11} = F ; +F(1).Fig = 'S3D' ; F(2).Fig = 'S3E' ; FigData{12} = F ; +F(1).Fig = 'S4A' ; F(2).Fig = 'S4B' ; FigData{13} = F ; +F(1).Fig = 'S4D' ; F(2).Fig = 'S4E' ; FigData{14} = F ; +F(1).Fig = 'S5A' ; F(2).Fig = 'S5B' ; FigData{15} = F ; +F(1).Fig = 'S5D' ; F(2).Fig = 'S5E' ; FigData{16} = F ; clear F + +F(1).Fig = 'S11A' ; F(1).Out = 1 ; F(1).Desc = 'Virions [log CEQ/mL]' ; +FigData{17} = F ; +F(1).Fig = 'S11B' ; +FigData{18} = F ; + +F(1).Fig = 'S6A' ; FigData{19} = F ; +F(1).Fig = 'S6B' ; FigData{20} = F ; +F(1).Fig = 'S6C' ; FigData{21} = F ; +F(1).Fig = 'S6D' ; FigData{22} = F ; +F(1).Fig = 'S6E' ; FigData{23} = F ; +F(1).Fig = 'S6F' ; FigData{24} = F ; clear F + +% LFC2 = log fold change wrt cohort 2 pre-treatment value +F(1).Fig = 'S10A' ; F(1).Out = 3 ; F(1).Desc = 'Virions [LFC2]' ; +F(2).Fig = 'S10B' ; F(2).Out = 23 ; F(2).Desc = 'Active S [LFC2]' ; +F(3).Fig = 'S10C' ; F(3).Out = 56 ; F(3).Desc = 'Regulation of killing [LFC2]' ; +FigData{25} = F ; + +F(1).Fig = 'S10E' ; +F(2).Fig = 'S10F' ; +F(3).Fig = 'S10G' ; +FigData{26} = F ; %% ======================================================================== % 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(:).ParSample] = deal([]) ;% parameter samples from Bayesian MCMC -[C(:).YicSample] = deal([]) ;% initial condition samples -[C(:).YSummary] = deal([]) ;% output summary statistics (95% CI) - -% load(name,'C') +[C(:).FigData] = deal(FigData{:}) ;% figure data structure +[C(:).nFails] = deal([]) ;% number of samples that failed (error thrown) disp('Collecting...') ; % For each model... -for n = 1:6 +for n = 1:26 + + disp(C(n).Scenario) ; - 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) - 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==6) - saveFile = 'N803_shared_Results' ; - modelFun = @(P) N803_shared(X{n}, {Dose{n},[]}, P, {[],[]}, 1,... - 'fullOut',true,'ivDosing',880*[.06 .06 .06 1 0]/.03) ; + c = 1+rem(n+1,2) ;% cohort # + switch n + case {1,2,5,6,17,18} + saveFile = 'N803_shared_2_Results' ; + modelFun = @(P) N803_shared_2(X{n}, P, ... + 'DoseTimes', Dose{n}, 'cohort', c, 'fullOut', true) ; + case 3 + saveFile = 'N803_shared_2_Results' ; + modelFun = @(P) N803_shared_2(X{n}, P, ... + 'DoseTimes', Dose{n}, 'cohort', 1, 'fullOut', true, ... + 'IVdoses', 880*[.06 0]/.03,'isNaive',true) ; + case 4 + saveFile = 'N803_shared_2_Results' ; + modelFun = @(P) N803_shared_2(X{n}, P, ... + 'DoseTimes', Dose{n}, 'cohort', 1, 'fullOut', true, ... + 'IVdoses', 880*[.06 .06 .06 1 0]/.03) ; + case {7,8} + saveFile = ['N803_cohort_' num2str(c) '_Results'] ; + modelFun = @(P) N803_single_2(X{n}, P, c, ... + 'DoseTimes', Dose{n}, 'fullOut', true) ; + case {9,10} + saveFile = 'N803_shared_2A_Results' ; + modelFun = @(P) N803_shared_2A(X{n}, P, ... + 'DoseTimes', Dose{n}, 'cohort', c, 'fullOut', true) ; + case {11,12} + saveFile = 'N803_shared_2B_Results' ; + modelFun = @(P) N803_shared_2B(X{n}, P, ... + 'DoseTimes', Dose{n}, 'cohort', c, 'fullOut', true) ; + case {13,14} + saveFile = 'N803_shared_2C_Results' ; + modelFun = @(P) N803_shared_2C(X{n}, P, ... + 'DoseTimes', Dose{n}, 'cohort', c, 'fullOut', true) ; + case {15,16} + saveFile = 'N803_shared_2D_Results' ; + modelFun = @(P) N803_shared_2D(X{n}, P, ... + 'DoseTimes', Dose{n}, 'cohort', c, 'fullOut', true) ; + case {19,20,21,22,23,24,25,26} + saveFile = 'N803_shared_2_Results' ; + modelFun = @(P) N803_shared_2(X{n}, P, ... + 'DoseTimes', Dose{n}, 'cohort', 2, 'fullOut', true) ; end - 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) ;% number of thinned MCMC samples - - % Load parameters for initial condition and parameter calculation - % functions: 'N803_shared.m', 'N803_single.m' + % Load parameters for initial condition and constant calculation + % functions (listed above) warning off load(saveFile,'Results') - load([saveFile '_MCMC'],'res') + if UseMCMC{n} + load([saveFile '_MCMC'],'res') + nSam = length(res.logPost(1:end-1)) ;% number of MCMC samples + PARS = repmat(Results(1).ParameterFit,nSam,1) ;% parameter samples + ParID = Results(1).setup.SampleIndex ; + PARS(:,ParID) = 10.^res.par(:,1:end-1)' ; + + if SaveSample{n} % only save every 'nThin'th sample + PARS = PARS(nThin:nThin:nSam,:) ; + nSam = nSam/nThin ; + end + else + nSam = nMSLS ;% use top 'nMSLS' from MSLS + PARS = vertcat(Results(1:nSam).ParameterFit) ; + end warning on - PARS = repmat(Results(1).ParameterFit,nSam,1) ;% sample of parameters - ParID = Results(1).setup.SampleIndex ; - PARS(:,ParID) = 10.^res.par(:,SamID)' ; - % Evaluate model and collect outputs, along with parameter and initial - % conditions (ICs) for 'N803_model_2.m' (called by 'shared' & 'single') - Ysam = cell(nSam,1) ;% sample of model outputs - Pcoh = cell(nSam,1) ;% sample of cohort [ parameters , ICs ] + % Evaluate model and collect outputs, along with constants + Ysam = cell(nSam,1) ;% sample of model outputs + Pcoh = cell(nSam,1) ;% sample of cohort [ constants , ICs ] 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 + C(n).nFails = sum(~passID) ; nSam = sum(passID) ; Ysam = Ysam(passID) ; Pcoh = Pcoh(passID) ; PCOH = cat(1,Pcoh{:}) ; - if SavePars{n} ; C(n).ParSample = PCOH(:, 1:31 ) ; end - if SaveYics{n} ; C(n).YicSample = PCOH(:,32:end) ; end + if SavePars{n} ; C(n).ParSample = PCOH ; end % 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(:).Figure] = deal([]) ;% figure label - [Ystats(:).lo95] = deal([]) ;% 2.5th percentile - [Ystats(:).hi95] = deal([]) ;% 97.5th percentile + OUTPUTS = zeros(length(X{n}),nSam) ;% temp storage for output sample + F = C(n).FigData ;% structure for storing figure data + Outs = vertcat(F(:).Out) ;% output index (in model function) - if n == 1 - NORMFACTS = ones(nOuts,nSam) ;% matrix of normalization factors - end - for out = SaveOutput{n} - for sam = 1:nSam - OUTPUTS(:,sam) = Ysam{sam}(:,out) ; + for o = 1:length(Outs) + for sam = 1:nSam % collect output (from each sample) into 1 matrix + OUTPUTS(:,sam) = Ysam{sam}(:,Outs(o)) ; end - if any( n == 1:2 ) && any( out == [ 3:15, 18:19, 29:61, 63:64 ] ) - if n == 1 - NORMFACTS(out,:) = OUTPUTS(1,:) ; + + if strcmp(F(o).Fig(1),'5') % fold change wrt cohort 1 pre-treat + if c == 1 + NORMFACTS(o,:) = OUTPUTS(1,:) ; %#ok end - OUTPUTS = OUTPUTS - NORMFACTS(out,:) ; + OUTPUTS = log10( OUTPUTS ./ NORMFACTS(o,:) ) ; + end + if any( n == [25,26] ) % fold change wrt cohort 2 pre-treat + OUTPUTS = log10( OUTPUTS ./OUTPUTS(1,:) ) ; + end + + if ~SaveSample{n} + F(o).lo95 = prctile(OUTPUTS', 2.5)' ; + F(o).hi95 = prctile(OUTPUTS', 97.5)' ; + else + F(o).SETS = OUTPUTS ; end - Ystats(out).Figure = Figure{out}{n} ; - Ystats(out).lo95 = prctile(OUTPUTS', 2.5)' ; - Ystats(out).hi95 = prctile(OUTPUTS', 97.5)' ; end - C(n).YSummary = Ystats ;% output summary statistics + C(n).FigData = F ;% storing figure data end %% ======================================================================== diff --git a/N803_model_2.m b/N803_model_2.m index 63a29fc..0810284 100644 --- a/N803_model_2.m +++ b/N803_model_2.m @@ -1,22 +1,19 @@ %% N803_model_2.m - solves ODE model for N-803 treatment of SIV % % /--------------------------------------------------------------\ -% | Date: 01/08/2022 | +% | Date: 05/06/2023 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | % | Pienaar Computational Systems Pharmacology Lab | % \--------------------------------------------------------------/ % -% Nomenclature: V = SIV virions [#/uL] -% T8 = total CD8+ T cells [#/uL] -% S0 = resting SIV-specific CD8+ T cells [#/uL] -% Sa = active SIV-specific CD8+ T cells [#/uL] -% N0 = resting non-SIV-specific CD8+ T cells [#/uL] -% Na = active non-SIV-specific CD8+ T cells [#/uL] -% X = N803 at absorption site [pmol/kg] -% C = N803 plasma concentration [pM] -% R = regulation [] (dimensionless quantity) +% Nomenclature: V = SIV virions [#/mL] +% S = SIV-specific CD8+ T cells [#/uL] +% N = non-SIV-specific CD8+ T cells [#/uL] +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = immune regulation [] (dimensionless quantity) % %% ======================================================================== % INPUTS @@ -25,114 +22,100 @@ % SoluTimes = ascending vector of days at which to evaluate solution % % DoseTimes = ascending vector of days at which to administer doses -% (elements of 'DoseTimes' must also be in 'SoluTimes') +% ( members of 'DoseTimes' must be members of 'SoluTimes' ) % -% Pars = vector of parameters (see list on line 96) +% Pars = vector of parameters (see list on line 90) % -% Yic = vector of initial conditions (see list on line 273) +% Yic = vector of initial conditions (see list on line 184) % %% ======================================================================== % OPTIONS % ======================================================================== % -% Opts = cell vector used to add optional inputs (see list on line 53) +% Opts = cell vector used to add optional inputs (see list on line 51) % EX: Opts = {'AbsTol',1e-2} will set ode solver absolute tolerance to 1e-2 % %% ======================================================================== % OUTPUTS % ======================================================================== % -% Y_OUT(:,1) = V at points in 'SoluTimes' [log fold change] -% Y_OUT(:,2) = T8 at points in 'SoluTimes' [fold change] +% Y_OUT(:,1) = SIV virions [log(#/mL)] at points in 'SoluTimes' +% Y_OUT(:,2) = total CD8+ T cells [#/uL] at points in 'SoluTimes' % %% ======================================================================== % FUNCTION % ======================================================================== function Y_OUT = N803_model_2(SoluTimes,DoseTimes,Pars,Yic,Opts) -% Declare optional input default values -EqualPars = [] ;% vector of indices in 'Pars' for equating parameters - % EX: EqualPars(28) = 27 will set P(28) = P(27) -fullOut = false ;% if true, more outputs are added (see line 183) -ivDosing = [] ;% vector of intravenous dose strengths -isNaive = false ;% if true, all IC are zero (except for N0) -timeOut = 60 ;% scalar time limit [seconds] for ODE solver +% declare optional input default values +fullOut = false ;% if true, more outputs are added to Y_OUT (line 346) +IVdoses = [] ;% vector of IV dose strengths (length of 'DoseTimes') +timeOut = 30 ;% scalar time limit [seconds] for ODE solver % if exceeded, an error is thrown -warnOff = false ;% if true, warning messages are suppressed +warnOff = true ;% if true, warning messages are suppressed AbsTol = 1e-6 ;% absolute error tolerance (see 'odeset') RelTol = 1e-3 ;% relative error tolerance (see 'odeset') odeSolver = '15s' ;% ode solver function: 'ode15s' or 'ode23s' -% Overwrite optional inputs, if specified +% overwrite optional inputs, if specified for n = 1:length(Opts) -if strcmp('EqualPars', Opts(n)) ; EqualPars = Opts{n+1} ; end -if strcmp('fullOut', Opts(n)) ; fullOut = Opts{n+1} ; end -if strcmp('ivDosing', Opts(n)) ; ivDosing = Opts{n+1} ; end -if strcmp('isNaive', Opts(n)) ; isNaive = Opts{n+1} ; end -if strcmp('timeOut', Opts(n)) ; timeOut = Opts{n+1} ; end -if strcmp('warnOff', Opts(n)) ; warnOff = Opts{n+1} ; end -if strcmp('AbsTol', Opts(n)) ; AbsTol = Opts{n+1} ; end -if strcmp('RelTol', Opts(n)) ; RelTol = Opts{n+1} ; end -if strcmp('odeSolver', Opts(n)) ; odeSolver = Opts{n+1} ; end + if strcmp('fullOut', Opts(n)) ; fullOut = Opts{n+1} ; end + if strcmp('IVdoses', Opts(n)) ; IVdoses = Opts{n+1} ; end + if strcmp('timeOut', Opts(n)) ; timeOut = Opts{n+1} ; end + if strcmp('warnOff', Opts(n)) ; warnOff = Opts{n+1} ; end + if strcmp('AbsTol', Opts(n)) ; AbsTol = Opts{n+1} ; end + if strcmp('RelTol', Opts(n)) ; RelTol = Opts{n+1} ; end + if strcmp('odeSolver', Opts(n)) ; odeSolver = Opts{n+1} ; end end -% check inputs vector orientation +% check input vector orientation if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end if size(DoseTimes ,1)==1 ; DoseTimes = DoseTimes' ; end if size(Yic ,2)==1 ; Yic = Yic' ; end +% check that members of 'DoseTimes' are members of 'SoluTimes' +if length( unique( [ SoluTimes ; DoseTimes ] ) ) > length( SoluTimes ) + error('All members of DoseTimes must be members of SoluTimes!') +end + % declare index in 'SoluTimes' for start, dose, and end times Pieces = [ 1 ; find(ismember(SoluTimes,DoseTimes)) ; length(SoluTimes) ] ; -% equate parameters in 'P' based on indices in 'EquatedPars' -if ~isempty(EqualPars) ; Pars = Pars(EqualPars) ; end - % suppress warnings (if requested) if warnOff ; warning off ; end %% ------------------------------------------------------------------------ % Rename inputed parameters ----------------------------------------------- -q = Pars(01) ;% V growth rate (if S+N were absent) [/d] -gS0 = Pars(02) ;% Sa killing rate constant [uL/#-d] -gN0 = Pars(03) ;% Na killing rate constant [uL/#-d] - -V50S = Pars(04) ;% 50% viral stimulation saturation for S [#/uL] -V50N = Pars(05) ;% 50% viral stimulation saturation for N [#/uL] -aS0 = Pars(06) ;% S0 activation rate constant [/d] -aN0 = Pars(07) ;% N0 activation rate constant [/d] -mS = Pars(08) ;% Sa reversion rate constant [/d] -mN = Pars(09) ;% Na reversion rate constant [/d] - -SN50 = Pars(10) ;% 50% S+N proliferation saturation [#/uL] -p0 = Pars(11) ;% S0/N0 proliferation rate constant [/d] -pS = Pars(12) ;% Sa proliferation rate constant [/d] -pN = Pars(13) ;% Na proliferation rate constant [/d] -d = Pars(14) ;% S0/N0 death rate constant [/d] -dA = Pars(15) ;% Sa/Na death rate constant [/d] - -Xi = Pars(16) ;% X initial condition [pmol/kg] -ka = Pars(17) ;% N803 absorption rate constant [/d] -ke = Pars(18) ;% N803 elimination rate constant [/d] -vd = Pars(19) ;% N803 'volume of distribution'/'bioavailability' [L/kg] -C50 = Pars(20) ;% 50% N803 stimulation concentration [pM] -p1 = Pars(21) ;% S0/N0 proliferation stimulation factor [] -aS1 = Pars(22) ;% S activation stimulation factor [] -aN1 = Pars(23) ;% N activation stimulation factor [] - -sS = Pars(24) ;% R generation due to S0 activation [/d] -sN = Pars(25) ;% R generation due to N0 activation [/d] -dR = Pars(26) ;% R decay rate constant [/d] -gS2 = Pars(27) ;% S killing regulation factor [] -gN2 = Pars(28) ;% N killing regulation factor [] -p2 = Pars(29) ;% S0/N0 proliferation regulation factor [] -aS2 = Pars(30) ;% S activation regulation factor [] -aN2 = Pars(31) ;% N activation regulation factor [] - -if isNaive % If SIV-naive, calculate initial N0 - Yic = Yic*0 ; - Yic(11) = SN50*( p0/d - 1 ) ; -end +q = Pars(01) ;% V growth rate constant [/d] +g = Pars(02) ;% S killing rate constant [uL/#-d] +V50S = Pars(03) ;% 50% V saturation for S activation [#/mL] +V50N = Pars(04) ;% 50% V saturation for N activation [#/mL] +aS = Pars(05) ;% S activation rate constant [/d] +aN = Pars(06) ;% N activation rate constant [/d] +mS = Pars(07) ;% S reversion rate constant [/d] +mN = Pars(08) ;% N reversion rate constant [/d] +S50 = Pars(09) ;% 50% S saturation for proliferation [#/uL] +N50 = Pars(10) ;% 50% N saturation for proliferation [#/uL] +p = Pars(11) ;% S,N resting proliferation rate constant [/d] +pA = Pars(12) ;% S active proliferation rate constant [/d] +d = Pars(13) ;% S,N resting death rate constant [/d] +dA = Pars(14) ;% S,N active death rate constant [/d] +Xi = Pars(15) ;% X initial condition (at dose) [pmol/kg] +ka = Pars(16) ;% N803 absorption rate constant [/d] +ke = Pars(17) ;% N803 elimination rate constant [/d] +vd = Pars(18) ;% N803 'volume of distribution'/'bioavailability' [L/kg] +C50 = Pars(19) ;% 50% N803 saturation (drug effects) [pM] +ro = Pars(20) ;% S,N proliferation stimulation factor [] +alS = Pars(21) ;% S activation stimulation factor [] +alN = Pars(22) ;% N activation stimulation factor [] +sS = Pars(23) ;% R generation rate constant (S activation) [/d] +sN = Pars(24) ;% R generation rate constant (N activation) [/d] +dR = Pars(25) ;% R decay rate constant [/d] +la = Pars(26) ;% S killing regulation factor [] +ph = Pars(27) ;% S,N proliferation regulation factor [] +zeS = Pars(28) ;% S activation regulation factor [] +zeN = Pars(29) ;% N activation regulation factor [] %% ------------------------------------------------------------------------ % Solve for each solution period ------------------------------------------ @@ -159,278 +142,268 @@ end Y(Piece,:) = Y_TEMP ; % add 'Y_TEMP' to full solution 'Y' - - Yic = Y_TEMP(end,:) ;% update ICs for next solution piece - if isempty(ivDosing) - Yic(14) = Yic(14) + Xi ;% add SC drug dose as initial condition + Yic = Y_TEMP(end,:) ;% update ICs for next solution piece + + if isempty(IVdoses) + Yic(13) = Yic(13) + Xi ;% add SC drug dose as initial condition else - Yic(15) = Yic(15) + ivDosing(n)/vd ;% add IV drug dose + Yic(14) = Yic(14) + IVdoses(n)/vd ;% add IV drug dose end end %% ------------------------------------------------------------------------ % Prepare main outputs 'Y_OUT' -------------------------------------------- -Y = abs(Y) ;% removing negatives resulting from numerical error +% Extremely small values can create numerical artifacts where the value +% becomes negative (mathematically impossible). This occurs for 'bad' +% parameter sets during calibration, and it occurs for the drug (where +% extremely small values are irrelevant). Taking the absolute value of Y +% prevents future errors with the 'log10' function. +Y = abs(Y) ; -Y_OUT = [ log10( Y(:,1) / Y(1,1) )... V [log fold change] - , sum(Y(:,2:13),2) / sum(Y(1,2:13)) ];% S+N [fold change] +Y_OUT = [ log10( Y(:,1) ) ... SIV virions [log(#/mL)] + , sum( Y(:,2:12), 2) ];% total CD8+ T cells [#/uL] %% ------------------------------------------------------------------------ % Add columns to 'Y_OUT' (if requested) ----------------------------------- -if ~fullOut ; return ; end - -Y_OUT = [ Y_OUT, Y, zeros(Pieces(end),46) ] ; - -V = Y(:,01) ;% SIV virions [#/uL] -S0 = Y(:,02) ;% resting specific CD8+ T cells [#/uL] -N0 = Y(:,11) ;% resting non-SIV-specific CD8+ T cells [#/uL] -Sa = sum( Y(:,03:10) ,2) ;% active specific CD8+ T cells [#/uL] -Na = sum( Y(:,12:13) ,2) ;% active non-SIV-specific T cells [#/uL] -Saf = Y(:,10) ;% last generation of Sa -Naf = Y(:,13) ;% last generation of Na -Sap = Sa - Saf ;% Sa that is proliferating -Nap = Y(:,12) ;% Na that is proliferating - -F = terms(Y) ;% see 'terms' function -p = F(:,13) ;% S0/N0 proliferation rate [/d] -aS = F(:,14) ;% S0 activation rate [/d] -aN = F(:,15) ;% N0 activation rate [/d] -gS = F(:,16) ;% Sa killing rate [#/uL-d] -gN = F(:,17) ;% Na killing rate [#/uL-d] - -Y_OUT(:,20) = F(:,01) ;% N803 effect [] -Y_OUT(:,21) = (S0+N0+Sa+Na) ;% CD8+ T cells [#/uL] -Y_OUT(:,22) = Sa ;% Sa -Y_OUT(:,23) = Na ;% Na -Y_OUT(:,24) = (S0+Sa) ;% S -Y_OUT(:,25) = (N0+Na) ;% N -Y_OUT(:,26) = (N0+Na)./(S0+N0+Sa+Na) ;% N frequency in (S+N) -Y_OUT(:,27) = (Sa)./(S0+Sa) ;% activated frequency in S -Y_OUT(:,28) = (Na)./(N0+Na) ;% activated frequency in N -Y_OUT(:,29) = log10( Sa ) ;% Sa (log) -Y_OUT(:,30) = log10( Na ) ;% Na (log) - -Y_OUT(:,31) = p.*S0 ;% S0 proliferation term [#/uL/d] -Y_OUT(:,32) = p.*N0 ;% N0 proliferation term [#/uL/d] -Y_OUT(:,33) = d.*S0 ;% S0 death term [#/uL/d] -Y_OUT(:,34) = d.*N0 ;% N0 death term [#/uL/d] -Y_OUT(:,35) = p ;% S0/N0 proliferation rate [/d] -Y_OUT(:,36) = F(:,02) ;% C multiplier for p* (log) -Y_OUT(:,37) = F(:,05) ;% R multiplier for p* (log) -Y_OUT(:,38) = F(:,10) ;% density-dep multiplier for p* (log) - -Y_OUT(:,39) = aS.*S0 ;% S0 activation term [#/uL/d] -Y_OUT(:,40) = aN.*N0 ;% N0 activation term [#/uL/d] -Y_OUT(:,41) = pS*Sap ;% Sa proliferation term [#/uL/d] -Y_OUT(:,42) = pN*Nap ;% Na proliferation term [#/uL/d] -Y_OUT(:,43) = dA*Sa ;% Sa death term [#/uL/d] -Y_OUT(:,44) = dA*Na ;% Na death term [#/uL/d] -Y_OUT(:,45) = mS*Saf ;% Sa reversion term [#/uL/d] -Y_OUT(:,46) = mN*Naf ;% Na reversion term [#/uL/d] -Y_OUT(:,47) = aS ;% S0 activation rate [/d] -Y_OUT(:,48) = aN ;% N0 activation rate [/d] -Y_OUT(:,49) = F(:,03) ;% C multiplier for aS* -Y_OUT(:,50) = ( F(:,04) + F(:,12) )./F(:,12) ;% C multiplier for aN* -Y_OUT(:,51) = F(:,06) ;% R multiplier for aS* -Y_OUT(:,52) = F(:,07) ;% R multiplier for aN* -Y_OUT(:,53) = F(:,11) ;% V multiplier for aS* -Y_OUT(:,54) = F(:,12) ;% V multiplier for aN* - -Y_OUT(:,55) = q*V ;% V growth term [#/uL/d] -Y_OUT(:,56) = gS.*Sa.*V ;% S killing term [#/uL/d] -Y_OUT(:,57) = gN.*Na.*V ;% N killing term [#/uL/d] -Y_OUT(:,58) = gS.*Sa ;% S killing rate [/d] -Y_OUT(:,59) = gN.*Na ;% N killing rate [/d] -Y_OUT(:,60) = F(:,08) ;% R multiplier for gS* -Y_OUT(:,61) = F(:,09) ;% R multiplier for gN* -Y_OUT(:,62) = gN.*Na ./ (gS.*Sa + gN.*Na) ;% N killing fraction - -Y_OUT(:,63) = F(:,18) ;% R gen by aS -Y_OUT(:,64) = F(:,19) ;% R gen by aN -Y_OUT(:,65) = F(:,19) ./ (F(:,18) + F(:,19)) ;% R gen by aN fraction - -% Converting most outputs to log-value -logID = [ 3:15 , 18:19 , 31:61 , 63:64 ] ; -Y_OUT(:,logID) = log10( Y_OUT(:,logID) ) ; - -% Converting virus to #/mL -Y_OUT(:,3) = Y_OUT(:,3) + 3 ; +if fullOut + Y_OUT = outputs(Y_OUT,Y) ; +end %% ======================================================================== -% NESTED FUNCTIONS +% 'system' nested function % ======================================================================== - -function dY = system(~,Y) - - if toc(solvertime) >= timeOut - error('ODE solver stalled.') + function dY = system(~,Y) + + if toc(solvertime) >= timeOut + error('ODE solver stalled.') + end + + V = Y(01) ;% SIV virions [#/mL] + S0 = Y(02) ;% resting SIV-specific CD8+ T cells [#/uL] + S = Y(03:10) ;% active SIV-specific CD8+ T cells [#/uL] (S1-S8) + N0 = Y(11) ;% resting non-SIV-specific CD8+ T cells [#/uL] + N1 = Y(12) ;% active non-SIV-specific CD8+ T cells [#/uL] + X = Y(13) ;% N803 at absorption site [pmol/kg] + C = Y(14) ;% N803 plasma concentration [pM] + R = Y(15:16) ;% immune regulation [] + F = terms(Y') ;% see 'terms' function + + PS = F(13) ;% S resting proliferation rate [/d] + PN = F(14) ;% N resting proliferation rate [/d] + AS = F(15) ;% S activation rate [/d] + AN = F(16) ;% N activation rate [/d] + GS = F(17) ;% S killing rate [/d] + SR = F(18)+F(19) ;% regulation generation rate [/d] + + dY = Y ;% dY/dt + + % calculate rates of change for V,S,N,X,C,R + dY(01) = q*V - GS*V ;% dV/dt + dY(02) = PS*S0 - d *S0 + mS*S(8) - AS*S0 ;% dS0/dt + dY(03) = - dA*S(1) + 2*AS*S0 - pA*S(1) ;% dS1/dt + dY(4:9)= - dA*S(2:7) + 2*pA*S(1:6) - pA*S(2:7) ;% dS2-7/dt + dY(10) = - dA*S(8) + 2*pA*S(7) - mS*S(8) ;% dS8/dt + dY(11) = PN*N0 - d *N0 + mN*N1 - AN*N0 ;% dN0/dt + dY(12) = - dA*N1 + 2*AN*N0 - mN*N1 ;% dN1/dt + dY(13) = - ka*X ;% dX/dt + dY(14) = ka*X/vd - ke*C ;% dC/dt + dY(15) = SR - dR*R(1) ;% dR1/dt + dY(16) = dR*R(1) - dR*R(2) ;% dR2/dt + end +%% ======================================================================== +% 'jacob' nested function +% ======================================================================== + function J = jacob(~,Y) + + V = Y(01) ;% SIV virions [#/mL] + S0 = Y(02) ;% resting SIV-specific CD8+ T cells [#/uL] + N0 = Y(11) ;% resting non-SIV-specific CD8+ T cells [#/uL] + C = Y(14) ;% N803 plasma concentration [pM] + F = terms(Y') ;% see 'terms' function + + PS = F(13) ;% S resting proliferation rate [/d] + PN = F(14) ;% N resting proliferation rate [/d] + AS = F(15) ;% S activation rate [/d] + AN = F(16) ;% N activation rate [/d] + GS = F(17) ;% S killing rate [/d] + + J = zeros(16) ; + + % Partial derivatives of PS,PN,AS,AN,GS w.r.t. V,S,N,C,R + dPSdR = - PS * ph * F(05) ; + dPNdR = - PN * ph * F(05) ; + dASdR = - AS * zeS * F(06) ; + dANdR = - AN * zeN * F(07) ; + dGSdR = - GS * la * F(08) ; + dPSdX = - PS * F(09) / S50 ; + dPNdX = - PN * F(10) / N50 ; + dASdV = AS * V50S / V / (V50S+V) ; + dANdV = aN * V50N / (V50N+V)^2 * F(07) ; + dOdC = C50 / (C50+C)^2 ; + dPSdC = p * F(09) * ro * dOdC * F(05) ; + dPNdC = p * F(10) * ro * dOdC * F(05) ; + dASdC = aS * F(11) * alS * dOdC * F(06) ; + dANdC = aN * alN * dOdC * F(07) ; + + % Partial derivatives of dV/dt + J(01,01 ) = q - GS ;% wrt V + J(01,03:10) = - g * F(08) * V ;% wrt SA + J(01,16 ) = - dGSdR * V ;% wrt R + + % Partial derivatives of dS0/dt + J(02,01 ) = - dASdV*S0 ;% wrt V + J(02,02 ) = dPSdX*S0 + PS - d - AS ;% wrt S0 + J(02,03:09) = dPSdX*S0 ;% wrt S1-S7 + J(02,10 ) = dPSdX*S0 + mS ;% wrt S8 + J(02,14 ) = dPSdC*S0 - dASdC*S0 ;% wrt C + J(02,16 ) = dPSdR*S0 - dASdR*S0 ;% wrt R + + % Partial derivatives of dS1-8/dt + J(03,01 ) = 2*dASdV*S0 ;% dS1/dt wrt V + J(03,02 ) = 2*AS ;% dS1/dt wrt S0 + J(03,14 ) = 2*dASdC*S0 ;% dS1/dt wrt C + J(03,16 ) = 2*dASdR*S0 ;% dS1/dt wrt R + J(35:17:137) = - dA - pA ;% dSi/dt wrt Si for i=1-7 + J(36:17:138) = 2*pA ;% dSi/dt wrt Si-1 for i=2-8 + J(10,10 ) = - dA - mS ;% dS8/dt wrt S8 + + % Partial derivatives of dN0/dt + J(11,01) = - dANdV*N0 ;% wrt V + J(11,11) = dPNdX*N0 + PN - d - AN ;% wrt N0 + J(11,12) = dPNdX*N0 + mN ;% wrt N1 + J(11,14) = dPNdC*N0 - dANdC*N0 ;% wrt C + J(11,16) = dPNdR*N0 - dANdR*N0 ;% wrt R + + % Partial derivatives of dNA/dt + J(12,01) = 2*dANdV*N0 ;% dN1/dt wrt V + J(12,11) = 2*AN ;% dN1/dt wrt N0 + J(12,12) = - dA - mN ;% dN1/dt wrt N1 + J(12,14) = 2*dANdC*N0 ;% dN1/dt wrt C + J(12,16) = 2*dANdR*N0 ;% dN1/dt wrt R + + % Partial derivatives of dX/dt + J(13,13) = -ka ;% wrt X + + % Partial derivatives of dC/dt + J(14,13) = ka/vd ;% wrt X + J(14,14) = -ke ;% wrt C + + % Partial derivatives of dR/dt + J(15,01) = sS*dASdV/aS/F(06) + sN*dANdV/aN/F(07) ;% dR1/dt wrt V + J(15,14) = sS*dASdC/aS/F(06) + sN*dANdC/aN/F(07) ;% dR1/dt wrt C + J(15,15) = -dR ;% dR1/dt wrt R1 + J(16,15) = dR ;% dR2/dt wrt R1 + J(16,16) = -dR ;% dR2/dt wrt R2 + end +%% ======================================================================== +% 'terms' nested function +% ======================================================================== + function F = terms(Y) + + V = Y(:,01) ;% SIV virions [#/mL] + S = sum(Y(:,02:10),2) ;% SIV-specific CD8+ T cells [#/uL] + N = sum(Y(:,11:12),2) ;% non-SIV-specific CD8+ T cells [#/uL] + SA = sum(Y(:,03:10),2) ;% active S [#/uL] + C = Y(:,14) ;% N803 plasma concentration [pM] + R = Y(:,16) ;% regulatory inhibition [] + + % percolate vector for storing terms + F = zeros(size(Y,1),19) ; + + F(:,01) = C./ (C50+C) ;% N803 effect [] + + F(:,02) = 1 + ro*F(:,1) ;% S,N proliferation stimulation [] + F(:,03) = 1 + alS*F(:,1) ;% S activation stimulation [] + F(:,04) = alN*F(:,1) ;% N803-induced N activation [] + + F(:,05) = 1./ (1+ ph*R) ;% S,N proliferation regulation [] + F(:,06) = 1./ (1+zeS*R) ;% S activation regulation [] + F(:,07) = 1./ (1+zeN*R) ;% N activation regulation [] + F(:,08) = 1./ (1+la*R) ;% S killing regulation [] + + F(:,09) = S50./(S50+S) ;% S proliferation density dependence [] + F(:,10) = N50./(N50+N) ;% N proliferation density dependence [] + F(:,11) = V./ (V50S+V) ;% S activation viral dependence [] + F(:,12) = V./ (V50N+V) ;% N activation viral dependence [] + + F(:,13) = p * F(:,09).* F(:,02).* F(:,05) ;% S proliferation [/d] + F(:,14) = p * F(:,10).* F(:,02).* F(:,05) ;% N proliferation [/d] + F(:,15) = aS* F(:,11).* F(:,03).* F(:,06) ;% S activation [/d] + F(:,16) = aN* (F(:,12)+F(:,04)).* F(:,07) ;% N activation [/d] + F(:,17) = g * SA.* F(:,08) ;% S killing rate [#/uL-d] + + F(:,18) = sS* F(:,11).* F(:,03) ;% R generation by AS [/d] + F(:,19) = sN* (F(:,12) + F(:,04)) ;% R generation by AN [/d] + end +%% ======================================================================== +% 'outputs' nested function +% ======================================================================== + function Y_OUT = outputs(Y_OUT,Y) + + Y_OUT = [ Y_OUT, Y, zeros( size(Y,1) , 41) ] ; + + V = Y(:,01) ;% SIV virions [#/mL] + S0 = Y(:,02) ;% resting specific CD8+ T cells [#/uL] + N0 = Y(:,11) ;% resting non-SIV-specific CD8+ T cells [#/uL] + SA = sum( Y(:,03:10) ,2) ;% active specific CD8+ T cells [#/uL] + NA = Y(:,12) ;% active non-SIV-specific T cells [#/uL] + Saf = Y(:,10) ;% last generation of SA + Sap = SA - Saf ;% SA that is proliferating + + F = terms(Y) ;% see 'terms' function + PS = F(:,13) ;% S resting proliferation rate [/d] + PN = F(:,14) ;% N resting proliferation rate [/d] + AS = F(:,15) ;% S activation rate [/d] + AN = F(:,16) ;% N activation rate [/d] + GS = F(:,17) ;% S killing rate [/d] + + Y_OUT(:,19) = F(:,01) ;% N803 effect [] + Y_OUT(:,20) = (S0+N0+SA+NA) ;% total CD8+ T cells [#/uL] + Y_OUT(:,21) = (S0+SA) ;% total S + Y_OUT(:,22) = (N0+NA) ;% total N + Y_OUT(:,23) = SA ;% active S + Y_OUT(:,24) = NA ;% active N + Y_OUT(:,25) = (S0+SA)./(S0+N0+SA+NA) ;% S frequency in (S+N) + Y_OUT(:,26) = SA./(S0+SA) ;% activated frequency in S + Y_OUT(:,27) = NA./(N0+NA) ;% activated frequency in N + + Y_OUT(:,28) = PS.*S0 ;% S0 proliferation term [#/uL/d] + Y_OUT(:,29) = PN.*N0 ;% N0 proliferation term [#/uL/d] + Y_OUT(:,30) = d.*S0 ;% S0 death term [#/uL/d] + Y_OUT(:,31) = d.*N0 ;% N0 death term [#/uL/d] + Y_OUT(:,32) = PS ;% S0 proliferation rate [/d] + Y_OUT(:,33) = PN ;% N0 proliferation rate [/d] + Y_OUT(:,34) = F(:,02) ;% C multiplier for PS,PN + Y_OUT(:,35) = F(:,05) ;% R multiplier for PS,PN + Y_OUT(:,36) = F(:,09) ;% density-dep multiplier for PS + Y_OUT(:,37) = F(:,10) ;% density-dep multiplier for PN + + Y_OUT(:,38) = AS.*S0 ;% S0 activation term [#/uL/d] + Y_OUT(:,39) = AN.*N0 ;% N0 activation term [#/uL/d] + Y_OUT(:,40) = pA*Sap ;% SA proliferation term [#/uL/d] + Y_OUT(:,41) = dA*SA ;% SA death term [#/uL/d] + Y_OUT(:,42) = dA*NA ;% NA death term [#/uL/d] + Y_OUT(:,43) = mS*Saf ;% SA reversion term [#/uL/d] + Y_OUT(:,44) = mN*NA ;% NA reversion term [#/uL/d] + Y_OUT(:,45) = AS ;% S0 activation rate [/d] + Y_OUT(:,46) = AN ;% N0 activation rate [/d] + Y_OUT(:,47) = F(:,03) ;% C multiplier for AS + Y_OUT(:,48) = ( F(:,04) + F(:,12) )./F(:,12) ;% C multiplier for AN + Y_OUT(:,49) = F(:,06) ;% R multiplier for AS + Y_OUT(:,50) = F(:,07) ;% R multiplier for AN + Y_OUT(:,51) = F(:,11) ;% V multiplier for AS + Y_OUT(:,52) = F(:,12) ;% V multiplier for AN + + Y_OUT(:,53) = q*V ;% V growth term [#/mL/d] + Y_OUT(:,54) = GS.*V ;% S killing term [#/mL/d] + Y_OUT(:,55) = GS ;% S killing rate [/d] + Y_OUT(:,56) = F(:,08) ;% R multiplier for GS + + Y_OUT(:,57) = F(:,18) ;% R gen by AS + Y_OUT(:,58) = F(:,19) ;% R gen by AN + Y_OUT(:,59) = F(:,19) ./ (F(:,18) + F(:,19));% R gen by AN fraction end - - V = Y(01) ;% SIV virions [#/uL] - S0 = Y(02) ;% resting specific CD8+ T cells [#/uL] - Sa = Y(03:10) ;% active specific CD8+ T cells [#/uL] (S1 to S8) - N0 = Y(11) ;% resting non-SIV-specific CD8+ T cells [#/uL] - Na = Y(12:13) ;% active non-SIV-specific T cells [#/uL] (N1 to N2) - X = Y(14) ;% N803 at absorption site [pmol/kg] - C = Y(15) ;% N803 plasma concentration [pM] - R = Y(16:17) ;% immune regulation [] - F = terms(Y') ;% modifiers for 'gS0','gN0','p0','aS0','aN0' - - p = F(13) ;% S0/N0 proliferation rate [/d] - aS = F(14) ;% S0 activation rate [/d] - aN = F(15) ;% N0 activation rate [/d] - gS = F(16) ;% Sa killing rate [#/uL-d] - gN = F(17) ;% Na killing rate [#/uL-d] - s = F(18)+F(19) ;% regulation generation rate [/d] - - dY = Y ;% dY/dt - - % calculate rates of change for V,S,N,X,C,R - dY(01) = q*V - gS*sum(Sa)*V - gN*sum(Na)*V ;% dV/dt - dY(02) = p*S0 - aS*S0 - d*S0 + mS*Sa(8) ;% dS0/dt - dY(03) = + 2*aS*S0 - dA*Sa(1) - pS*Sa(1) ;% dS1/dt - dY(4:9)= 2*pS*Sa(1:6) - dA*Sa(2:7) - pS*Sa(2:7) ;% dS2-7/dt - dY(10) = 2*pS*Sa(7) - dA*Sa(8) - mS*Sa(8) ;% dS8/dt - dY(11) = p*N0 - aN*N0 - d*N0 + mN*Na(2) ;% dN0/dt - dY(12) = + 2*aN*N0 - dA*Na(1) - pN*Na(1) ;% dN1/dt - dY(13) = 2*pN*Na(1) - dA*Na(2) - mN*Na(2) ;% dN2/dt - dY(14) = - ka*X ;% dX/dt - dY(15) = ka*X/vd - ke*C ;% dC/dt - dY(16) = s - dR*R(1) ;% dR1/dt - dY(17) = dR*R(1) - dR*R(2) ;% dR2/dt - -end - -function J = jacob(~,Y) - - V = Y(01) ;% SIV virions [#/uL] - S0 = Y(02) ;% resting specific CD8+ T cells [#/uL] - Sa = sum(Y(03:10)) ;% active specific CD8+ T cells [#/uL] - N0 = Y(11) ;% resting non-SIV-specific CD8+ T cells [#/uL] - Na = Y(12)+Y(13) ;% active non-SIV-specific T cells [#/uL] - C = Y(15) ;% N803 plasma concentration [pM] - F = terms(Y') ;% modifiers for 'gS0','gN0','p0','aS0','aN0' - - p = F(13) ;% S0/N0 proliferation rate [/d] - aS = F(14) ;% S0 activation rate [/d] - aN = F(15) ;% N0 activation rate [/d] - gS = F(16) ;% Sa killing rate [#/uL-d] - gN = F(17) ;% Na killing rate [#/uL-d] - - J = zeros(17) ; - - % Partial derivatives of gS,gN,p,aS,aN w.r.t. V,S,N,C,R - dpdR = - p * p2 * F(05) ; - daSdR = - aS * aS2 * F(06) ; - daNdR = - aN * aN2 * F(07) ; - dgSdR = - gS * gS2 * F(08) ; - dgNdR = - gN * gN2 * F(09) ; - dpdSN = - p * F(10) / SN50 ; - daSdV = aS0 * V50S / (V50S+V)^2 * F(03) * F(06) ; - daNdV = aN0 * V50N / (V50N+V)^2 * F(07) ; - dOdC = C50 / (C50+C)^2 ; - dpdC = p0 * F(10) * p1 * dOdC * F(05) ; - daSdC = aS0 * F(11) * aS1 * dOdC * F(06) ; - daNdC = aN0 * aN1 * dOdC * F(07) ; - - % Partial derivatives of f1 = dV/dt - J(01,01 ) = q - gS*Sa - gN*Na ; - J(01,02:10) = - gS * V ; - J(01,11:13) = - gN * V ; - J(01,17 ) = - ( dgSdR*Sa + dgNdR*Na ) * V ; - - % Partial derivatives of f2 = dS0/dt - J(02,01 ) = - daSdV * S0 ; - J(02,02 ) = p + dpdSN*S0 - d - aS ; - J(02,03:13) = dpdSN*S0 ; - J(02,10 ) = dpdSN*S0 + mS ; - J(02,15 ) = ( dpdC - daSdC ) * S0 ; - J(02,17 ) = ( dpdR - daSdR ) * S0 ; - - % Partial derivatives of f3 = dS1/dt (save f3/dS1) - J(03,01 ) = 2*daSdV * S0 ; - J(03,02 ) = 2*aS ; - J(03,15 ) = 2*daSdC * S0 ; - J(03,17 ) = 2*daSdR * S0 ; - - % Partial derivatives of f4-10 = dS2-8/dt - J(37:18:145) = - dA - pS ; - J(38:18:146) = 2*pS ; - J(10,10) = - dA - mS ; - - % Partial derivatives of f11 = dN0/dt - J(11,01 ) = - daNdV*N0 ; - J(11,02:12) = dpdSN*N0 ; - J(11,11 ) = p + dpdSN*N0 - d - aN ; - J(11,13 ) = dpdSN*N0 + mN ; - J(11,15 ) = ( dpdC - daNdC ) * N0 ; - J(11,17 ) = ( dpdR - daNdR ) * N0 ; - - % Partial derivatives of f12 = dN1/dt - J(12,01) = 2*daNdV * N0 ; - J(12,11) = 2*aN ; - J(12,12) = - dA - pN ; - J(12,15) = 2*daNdC * N0 ; - J(12,17) = 2*daNdR * N0 ; - - % Partial derivatives of f12-13 = dN1-2/dt - J(13,12) = 2*pN ; - J(13,13) = - dA - mN ; - - % Partial derivatives of f14 = dX/dt - J(14,14) = -ka ; - - % Partial derivatives of f15 = dC/dt - J(15,14) = ka/vd ; - J(15,15) = -ke ; - - % Partial derivatives of f16 = dR1/dt - J(16,01) = sS/aS0*daSdV/F(06) + sN/aN0*daNdV/F(07) ; - J(16,15) = sS/aS0*daSdC/F(06) + sN/aN0*daNdC/F(07) ; - J(16,16) = -dR ; - - % Partial derivatives of f17 = dR2/dt - J(17,16) = dR ; - J(17,17) = -dR ; -end - -function F = terms(Y) - - V = Y(:,01) ;% SIV virions [#/uL] - SN = sum(Y(:,02:13),2) ;% CD8+ T cells [#/uL] - C = Y(:,15) ;% N803 plasma concentration [pM] - R = Y(:,17) ;% immune regulation [] - - % percolate vector for storing terms - F = zeros(size(Y,1),19) ; - - F(:,01) = C./ (C50+C) ;% N803 effect [] - - F(:,02) = 1 + p1*F(:,1) ;% S0/N0 proliferation stimulation [] - F(:,03) = 1 + aS1*F(:,1) ;% S activation stimulation [] - F(:,04) = aN1*F(:,1) ;% N803-induced N activation [] - - F(:,05) = 1./ (1+ p2*R) ;% S0/N0 proliferation regulation [] - F(:,06) = 1./ (1+aS2*R) ;% Sa activation regulation [] - F(:,07) = 1./ (1+aN2*R) ;% Na activation regulation [] - F(:,08) = 1./ (1+gS2*R) ;% Sa killing regulation [] - F(:,09) = 1./ (1+gN2*R) ;% Na killing regulation [] - - F(:,10) = SN50./(SN50+SN) ;% S0/N0 proliferation density dependence [] - F(:,11) = V./ (V50S+V) ;% Sa activation antigen dependence [] - F(:,12) = V./ (V50N+V) ;% Na activation antigen dependence [] - - F(:,13) = p0 * F(:,10).* F(:,02).* F(:,05) ;% S0/N0 proliferation [/d] - F(:,14) = aS0* F(:,11).* F(:,03).* F(:,06) ;% S0 activation rate [/d] - F(:,15) = aN0* (F(:,12)+F(:,04)).* F(:,07) ;% N0 activation rate [/d] - F(:,16) = gS0* F(:,08) ;% Sa killing rate [#/uL-d] - F(:,17) = gN0* F(:,09) ;% Na killing rate [#/uL-d] - - F(:,18) = sS* F(:,11).* F(:,03) ;% regulation generation by aS [/d] - F(:,19) = sN* (F(:,12) + F(:,04)) ;% regulation generation by aN [/d] -end - end \ No newline at end of file diff --git a/N803_model_2A.m b/N803_model_2A.m new file mode 100644 index 0000000..9c95f6e --- /dev/null +++ b/N803_model_2A.m @@ -0,0 +1,412 @@ +%% N803_model_2A.m - solves ODE model for N-803 treatment of SIV +% +% /--------------------------------------------------------------\ +% | Date: 05/06/2023 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Model 2A is a version of model 2 where SIV-specific cells divide 6 times +% after activation, instead of 8. +% +% Nomenclature: V = SIV virions [#/mL] +% S = SIV-specific CD8+ T cells [#/uL] +% N = non-SIV-specific CD8+ T cells [#/uL] +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = immune regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% +% SoluTimes = ascending vector of days at which to evaluate solution +% +% DoseTimes = ascending vector of days at which to administer doses +% ( members of 'DoseTimes' must be members of 'SoluTimes' ) +% +% Pars = vector of parameters (see list on line 93) +% +% Yic = vector of initial conditions (see list on line 187) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% +% Opts = cell vector used to add optional inputs (see list on line 54) +% EX: Opts = {'AbsTol',1e-2} will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% +% Y_OUT(:,1) = SIV virions [log(#/mL)] at points in 'SoluTimes' +% Y_OUT(:,2) = total CD8+ T cells [#/uL] at points in 'SoluTimes' +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function Y_OUT = N803_model_2A(SoluTimes,DoseTimes,Pars,Yic,Opts) + +% declare optional input default values +fullOut = false ;% if true, more outputs are added to Y_OUT (line 349) +IVdoses = [] ;% vector of IV dose strengths (length of 'DoseTimes') +timeOut = 30 ;% scalar time limit [seconds] for ODE solver + % if exceeded, an error is thrown +warnOff = true ;% if true, warning messages are suppressed +AbsTol = 1e-6 ;% absolute error tolerance (see 'odeset') +RelTol = 1e-3 ;% relative error tolerance (see 'odeset') +odeSolver = '15s' ;% ode solver function: 'ode15s' or 'ode23s' + +% overwrite optional inputs, if specified +for n = 1:length(Opts) + if strcmp('fullOut', Opts(n)) ; fullOut = Opts{n+1} ; end + if strcmp('IVdoses', Opts(n)) ; IVdoses = Opts{n+1} ; end + if strcmp('timeOut', Opts(n)) ; timeOut = Opts{n+1} ; end + if strcmp('warnOff', Opts(n)) ; warnOff = Opts{n+1} ; end + if strcmp('AbsTol', Opts(n)) ; AbsTol = Opts{n+1} ; end + if strcmp('RelTol', Opts(n)) ; RelTol = Opts{n+1} ; end + if strcmp('odeSolver', Opts(n)) ; odeSolver = Opts{n+1} ; end +end + +% check input vector orientation +if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end +if size(DoseTimes ,1)==1 ; DoseTimes = DoseTimes' ; end +if size(Yic ,2)==1 ; Yic = Yic' ; end + +% check that members of 'DoseTimes' are members of 'SoluTimes' +if length( unique( [ SoluTimes ; DoseTimes ] ) ) > length( SoluTimes ) + error('All members of DoseTimes must be members of SoluTimes!') +end + +% declare index in 'SoluTimes' for start, dose, and end times +Pieces = [ 1 ; find(ismember(SoluTimes,DoseTimes)) ; length(SoluTimes) ] ; + +% suppress warnings (if requested) +if warnOff ; warning off ; end + +%% ------------------------------------------------------------------------ +% Rename inputed parameters ----------------------------------------------- + +q = Pars(01) ;% V growth rate constant [/d] +g = Pars(02) ;% S killing rate constant [uL/#-d] +V50S = Pars(03) ;% 50% V saturation for S activation [#/mL] +V50N = Pars(04) ;% 50% V saturation for N activation [#/mL] +aS = Pars(05) ;% S activation rate constant [/d] +aN = Pars(06) ;% N activation rate constant [/d] +mS = Pars(07) ;% S reversion rate constant [/d] +mN = Pars(08) ;% N reversion rate constant [/d] +S50 = Pars(09) ;% 50% S saturation for proliferation [#/uL] +N50 = Pars(10) ;% 50% N saturation for proliferation [#/uL] +p = Pars(11) ;% S,N resting proliferation rate constant [/d] +pA = Pars(12) ;% S active proliferation rate constant [/d] +d = Pars(13) ;% S,N resting death rate constant [/d] +dA = Pars(14) ;% S,N active death rate constant [/d] +Xi = Pars(15) ;% X initial condition (at dose) [pmol/kg] +ka = Pars(16) ;% N803 absorption rate constant [/d] +ke = Pars(17) ;% N803 elimination rate constant [/d] +vd = Pars(18) ;% N803 'volume of distribution'/'bioavailability' [L/kg] +C50 = Pars(19) ;% 50% N803 saturation (drug effects) [pM] +ro = Pars(20) ;% S,N proliferation stimulation factor [] +alS = Pars(21) ;% S activation stimulation factor [] +alN = Pars(22) ;% N activation stimulation factor [] +sS = Pars(23) ;% R generation rate constant (S activation) [/d] +sN = Pars(24) ;% R generation rate constant (N activation) [/d] +dR = Pars(25) ;% R decay rate constant [/d] +la = Pars(26) ;% S killing regulation factor [] +ph = Pars(27) ;% S,N proliferation regulation factor [] +zeS = Pars(28) ;% S activation regulation factor [] +zeN = Pars(29) ;% N activation regulation factor [] + +%% ------------------------------------------------------------------------ +% Solve for each solution period ------------------------------------------ + +Y = zeros(Pieces(end),length(Yic)) ;% percolate output matrix + +odeOptions = odeset('AbsTol',AbsTol,'RelTol',RelTol,'Jacobian',@jacob) ; + +for n = 1:length(Pieces)-1 + + solvertime = tic ;% timer (for catching a stalled execution) + + Piece = Pieces(n):Pieces(n+1) ;% index for solution piece + + switch odeSolver % solve ODE for time window indexed by Piece + case '15s' + [~,Y_TEMP] = ode15s(@system,SoluTimes(Piece),Yic,odeOptions) ; + case '23s' + [~,Y_TEMP] = ode23s(@system,SoluTimes(Piece),Yic,odeOptions) ; + end + + if length(Piece) == 2 % exclude timepoints added by ODE solver + Y_TEMP = [ Y_TEMP(1,:) ; Y_TEMP(end,:) ] ; + end + + Y(Piece,:) = Y_TEMP ; % add 'Y_TEMP' to full solution 'Y' + + Yic = Y_TEMP(end,:) ;% update ICs for next solution piece + + if isempty(IVdoses) + Yic(11) = Yic(11) + Xi ;% add SC drug dose as initial condition + else + Yic(12) = Yic(12) + IVdoses(n)/vd ;% add IV drug dose + end +end + +%% ------------------------------------------------------------------------ +% Prepare main outputs 'Y_OUT' -------------------------------------------- + +% Extremely small values can create numerical artifacts where the value +% becomes negative (mathematically impossible). This occurs for 'bad' +% parameter sets during calibration, and it occurs for the drug (where +% extremely small values are irrelevant). Taking the absolute value of Y +% prevents future errors with the 'log10' function. +Y = abs(Y) ; + +Y_OUT = [ log10( Y(:,1) ) ... SIV virions [log(#/mL)] + , sum( Y(:,2:10), 2) ];% total CD8+ T cells [#/uL] + +%% ------------------------------------------------------------------------ +% Add columns to 'Y_OUT' (if requested) ----------------------------------- + +if fullOut + Y_OUT = outputs(Y_OUT,Y) ; +end + +%% ======================================================================== +% 'system' nested function +% ======================================================================== + function dY = system(~,Y) + + if toc(solvertime) >= timeOut + error('ODE solver stalled.') + end + + V = Y(01) ;% SIV virions [#/mL] + S0 = Y(02) ;% resting SIV-specific CD8+ T cells [#/uL] + S = Y(03:08) ;% active SIV-specific CD8+ T cells [#/uL] (S1-S6) + N0 = Y(09) ;% resting non-SIV-specific CD8+ T cells [#/uL] + N1 = Y(10) ;% active non-SIV-specific CD8+ T cells [#/uL] + X = Y(11) ;% N803 at absorption site [pmol/kg] + C = Y(12) ;% N803 plasma concentration [pM] + R = Y(13:14) ;% immune regulation [] + F = terms(Y') ;% see 'terms' function + + PS = F(13) ;% S resting proliferation rate [/d] + PN = F(14) ;% N resting proliferation rate [/d] + AS = F(15) ;% S activation rate [/d] + AN = F(16) ;% N activation rate [/d] + GS = F(17) ;% S killing rate [/d] + SR = F(18)+F(19) ;% regulation generation rate [/d] + + dY = Y ;% dY/dt + + % calculate rates of change for V,S,N,X,C,R + dY(01) = q*V - GS*V ;% dV/dt + dY(02) = PS*S0 - d *S0 + mS*S(6) - AS*S0 ;% dS0/dt + dY(03) = - dA*S(1) + 2*AS*S0 - pA*S(1) ;% dS1/dt + dY(4:7)= - dA*S(2:5) + 2*pA*S(1:4) - pA*S(2:5) ;% dS2-5/dt + dY(08) = - dA*S(6) + 2*pA*S(5) - mS*S(6) ;% dS6/dt + dY(09) = PN*N0 - d *N0 + mN*N1 - AN*N0 ;% dN0/dt + dY(10) = - dA*N1 + 2*AN*N0 - mN*N1 ;% dN1/dt + dY(11) = - ka*X ;% dX/dt + dY(12) = ka*X/vd - ke*C ;% dC/dt + dY(13) = SR - dR*R(1) ;% dR1/dt + dY(14) = dR*R(1) - dR*R(2) ;% dR2/dt + end +%% ======================================================================== +% 'jacob' nested function +% ======================================================================== + function J = jacob(~,Y) + + V = Y(01) ;% SIV virions [#/mL] + S0 = Y(02) ;% resting SIV-specific CD8+ T cells [#/uL] + N0 = Y(09) ;% resting non-SIV-specific CD8+ T cells [#/uL] + C = Y(12) ;% N803 plasma concentration [pM] + F = terms(Y') ;% see 'terms' function + + PS = F(13) ;% S resting proliferation rate [/d] + PN = F(14) ;% N resting proliferation rate [/d] + AS = F(15) ;% S activation rate [/d] + AN = F(16) ;% N activation rate [/d] + GS = F(17) ;% S killing rate [/d] + + J = zeros(14) ; + + % Partial derivatives of PS,PN,AS,AN,GS w.r.t. V,S,N,C,R + dPSdR = - PS * ph * F(05) ; + dPNdR = - PN * ph * F(05) ; + dASdR = - AS * zeS * F(06) ; + dANdR = - AN * zeN * F(07) ; + dGSdR = - GS * la * F(08) ; + dPSdX = - PS * F(09) / S50 ; + dPNdX = - PN * F(10) / N50 ; + dASdV = AS * V50S / V / (V50S+V) ; + dANdV = aN * V50N / (V50N+V)^2 * F(07) ; + dOdC = C50 / (C50+C)^2 ; + dPSdC = p * F(09) * ro * dOdC * F(05) ; + dPNdC = p * F(10) * ro * dOdC * F(05) ; + dASdC = aS * F(11) * alS * dOdC * F(06) ; + dANdC = aN * alN * dOdC * F(07) ; + + % Partial derivatives of dV/dt + J(01,01 ) = q - GS ;% wrt V + J(01,03:08) = - g * F(08) * V ;% wrt SA + J(01,14 ) = - dGSdR * V ;% wrt R + + % Partial derivatives of dS0/dt + J(02,01 ) = - dASdV*S0 ;% wrt V + J(02,02 ) = dPSdX*S0 + PS - d - AS ;% wrt S0 + J(02,03:07) = dPSdX*S0 ;% wrt S1-S5 + J(02,08 ) = dPSdX*S0 + mS ;% wrt S6 + J(02,12 ) = dPSdC*S0 - dASdC*S0 ;% wrt C + J(02,14 ) = dPSdR*S0 - dASdR*S0 ;% wrt R + + % Partial derivatives of dS1-6/dt + J(03,01 ) = 2*dASdV*S0 ;% dS1/dt wrt V + J(03,02 ) = 2*AS ;% dS1/dt wrt S0 + J(03,12 ) = 2*dASdC*S0 ;% dS1/dt wrt C + J(03,14 ) = 2*dASdR*S0 ;% dS1/dt wrt R + J(31:15:91) = - dA - pA ;% dSi/dt wrt Si for i=1-5 + J(32:15:92) = 2*pA ;% dSi/dt wrt Si-1 for i=2-6 + J(08,08 ) = - dA - mS ;% dS6/dt wrt S6 + + % Partial derivatives of dN0/dt + J(09,01) = - dANdV*N0 ;% wrt V + J(09,09) = dPNdX*N0 + PN - d - AN ;% wrt N0 + J(09,10) = dPNdX*N0 + mN ;% wrt N1 + J(09,12) = dPNdC*N0 - dANdC*N0 ;% wrt C + J(09,14) = dPNdR*N0 - dANdR*N0 ;% wrt R + + % Partial derivatives of dNA/dt + J(10,01) = 2*dANdV*N0 ;% dN1/dt wrt V + J(10,09) = 2*AN ;% dN1/dt wrt N0 + J(10,10) = - dA - mN ;% dN1/dt wrt N1 + J(10,12) = 2*dANdC*N0 ;% dN1/dt wrt C + J(10,14) = 2*dANdR*N0 ;% dN1/dt wrt R + + % Partial derivatives of dX/dt + J(11,11) = -ka ;% wrt X + + % Partial derivatives of dC/dt + J(12,11) = ka/vd ;% wrt X + J(12,12) = -ke ;% wrt C + + % Partial derivatives of dR/dt + J(13,01) = sS*dASdV/aS/F(06) + sN*dANdV/aN/F(07) ;% dR1/dt wrt V + J(13,12) = sS*dASdC/aS/F(06) + sN*dANdC/aN/F(07) ;% dR1/dt wrt C + J(13,13) = -dR ;% dR1/dt wrt R1 + J(14,13) = dR ;% dR2/dt wrt R1 + J(14,14) = -dR ;% dR2/dt wrt R2 + end +%% ======================================================================== +% 'terms' nested function +% ======================================================================== + function F = terms(Y) + + V = Y(:,01) ;% SIV virions [#/mL] + S = sum(Y(:,02:08),2) ;% SIV-specific CD8+ T cells [#/uL] + N = sum(Y(:,09:10),2) ;% non-SIV-specific CD8+ T cells [#/uL] + SA = sum(Y(:,03:08),2) ;% active S [#/uL] + C = Y(:,12) ;% N803 plasma concentration [pM] + R = Y(:,14) ;% regulatory inhibition [] + + % percolate vector for storing terms + F = zeros(size(Y,1),19) ; + + F(:,01) = C./ (C50+C) ;% N803 effect [] + + F(:,02) = 1 + ro*F(:,1) ;% S,N proliferation stimulation [] + F(:,03) = 1 + alS*F(:,1) ;% S activation stimulation [] + F(:,04) = alN*F(:,1) ;% N803-induced N activation [] + + F(:,05) = 1./ (1+ ph*R) ;% S,N proliferation regulation [] + F(:,06) = 1./ (1+zeS*R) ;% S activation regulation [] + F(:,07) = 1./ (1+zeN*R) ;% N activation regulation [] + F(:,08) = 1./ (1+la*R) ;% S killing regulation [] + + F(:,09) = S50./(S50+S) ;% S proliferation density dependence [] + F(:,10) = N50./(N50+N) ;% N proliferation density dependence [] + F(:,11) = V./ (V50S+V) ;% S activation viral dependence [] + F(:,12) = V./ (V50N+V) ;% N activation viral dependence [] + + F(:,13) = p * F(:,09).* F(:,02).* F(:,05) ;% S proliferation [/d] + F(:,14) = p * F(:,10).* F(:,02).* F(:,05) ;% N proliferation [/d] + F(:,15) = aS* F(:,11).* F(:,03).* F(:,06) ;% S activation [/d] + F(:,16) = aN* (F(:,12)+F(:,04)).* F(:,07) ;% N activation [/d] + F(:,17) = g * SA.* F(:,08) ;% S killing rate [#/uL-d] + + F(:,18) = sS* F(:,11).* F(:,03) ;% R generation by AS [/d] + F(:,19) = sN* (F(:,12) + F(:,04)) ;% R generation by AN [/d] + end +%% ======================================================================== +% 'outputs' nested function +% ======================================================================== + function Y_OUT = outputs(Y_OUT,Y) + + Y_OUT = [ Y_OUT, Y, zeros( size(Y,1) , 48) ] ; + + V = Y(:,01) ;% SIV virions [#/mL] + S0 = Y(:,02) ;% resting specific CD8+ T cells [#/uL] + N0 = Y(:,09) ;% resting non-SIV-specific CD8+ T cells [#/uL] + SA = sum( Y(:,03:08) ,2) ;% active specific CD8+ T cells [#/uL] + NA = Y(:,10) ;% active non-SIV-specific T cells [#/uL] + Saf = Y(:,08) ;% last generation of SA + Sap = SA - Saf ;% SA that is proliferating + + F = terms(Y) ;% see 'terms' function + PS = F(:,13) ;% S resting proliferation rate [/d] + PN = F(:,14) ;% N resting proliferation rate [/d] + AS = F(:,15) ;% S activation rate [/d] + AN = F(:,16) ;% N activation rate [/d] + GS = F(:,17) ;% S killing rate [/d] + + Y_OUT(:,19) = F(:,01) ;% N803 effect [] + Y_OUT(:,20) = (S0+N0+SA+NA) ;% total CD8+ T cells [#/uL] + Y_OUT(:,21) = (S0+SA) ;% total S + Y_OUT(:,22) = (N0+NA) ;% total N + Y_OUT(:,23) = SA ;% active S + Y_OUT(:,24) = NA ;% active N + Y_OUT(:,25) = (S0+SA)./(S0+N0+SA+NA) ;% S frequency in (S+N) + Y_OUT(:,26) = SA./(S0+SA) ;% activated frequency in S + Y_OUT(:,27) = NA./(N0+NA) ;% activated frequency in N + + Y_OUT(:,28) = PS.*S0 ;% S0 proliferation term [#/uL/d] + Y_OUT(:,29) = PN.*N0 ;% N0 proliferation term [#/uL/d] + Y_OUT(:,30) = d.*S0 ;% S0 death term [#/uL/d] + Y_OUT(:,31) = d.*N0 ;% N0 death term [#/uL/d] + Y_OUT(:,32) = PS ;% S0 proliferation rate [/d] + Y_OUT(:,33) = PN ;% N0 proliferation rate [/d] + Y_OUT(:,34) = F(:,02) ;% C multiplier for PS,PN + Y_OUT(:,35) = F(:,05) ;% R multiplier for PS,PN + Y_OUT(:,36) = F(:,09) ;% density-dep multiplier for PS + Y_OUT(:,37) = F(:,10) ;% density-dep multiplier for PN + + Y_OUT(:,38) = AS.*S0 ;% S0 activation term [#/uL/d] + Y_OUT(:,39) = AN.*N0 ;% N0 activation term [#/uL/d] + Y_OUT(:,40) = pA*Sap ;% SA proliferation term [#/uL/d] + Y_OUT(:,41) = dA*SA ;% SA death term [#/uL/d] + Y_OUT(:,42) = dA*NA ;% NA death term [#/uL/d] + Y_OUT(:,43) = mS*Saf ;% SA reversion term [#/uL/d] + Y_OUT(:,44) = mN*NA ;% NA reversion term [#/uL/d] + Y_OUT(:,45) = AS ;% S0 activation rate [/d] + Y_OUT(:,46) = AN ;% N0 activation rate [/d] + Y_OUT(:,47) = F(:,03) ;% C multiplier for AS + Y_OUT(:,48) = ( F(:,04) + F(:,12) )./F(:,12) ;% C multiplier for AN + Y_OUT(:,49) = F(:,06) ;% R multiplier for AS + Y_OUT(:,50) = F(:,07) ;% R multiplier for AN + Y_OUT(:,51) = F(:,11) ;% V multiplier for AS + Y_OUT(:,52) = F(:,12) ;% V multiplier for AN + + Y_OUT(:,53) = q*V ;% V growth term [#/mL/d] + Y_OUT(:,54) = GS.*V ;% S killing term [#/mL/d] + Y_OUT(:,55) = GS ;% S killing rate [/d] + Y_OUT(:,56) = F(:,08) ;% R multiplier for GS + + Y_OUT(:,57) = F(:,18) ;% R gen by AS + Y_OUT(:,58) = F(:,19) ;% R gen by AN + Y_OUT(:,59) = F(:,19) ./ (F(:,18) + F(:,19));% R gen by AN fraction + end +end \ No newline at end of file diff --git a/N803_model_2B.m b/N803_model_2B.m new file mode 100644 index 0000000..1f3900e --- /dev/null +++ b/N803_model_2B.m @@ -0,0 +1,347 @@ +%% N803_model_2B.m - solves ODE model for N-803 treatment of SIV +% +% /--------------------------------------------------------------\ +% | Date: 05/07/2023 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Model 2B is a version of model 2 that models all CD8+ T cells as the +% SIV-specific cells from Model 2, which now only divide 4 times after +% activation. Also, regulation is generated with independent parameters. +% +% Nomenclature: V = SIV virions [#/mL] +% E = CD8+ T cells [#/uL] +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = immune regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% +% SoluTimes = ascending vector of days at which to evaluate solution +% +% DoseTimes = ascending vector of days at which to administer doses +% ( members of 'DoseTimes' must be members of 'SoluTimes' ) +% +% Pars = vector of parameters (see list on line 91) +% +% Yic = vector of initial conditions (see list on line 172) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% +% Opts = cell vector used to add optional inputs (see list on line 54) +% EX: Opts = {'AbsTol',1e-2} will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% +% Y_OUT(:,1) = SIV virions [log(#/mL)] at points in 'SoluTimes' +% Y_OUT(:,2) = total CD8+ T cells [#/uL] at points in 'SoluTimes' +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function Y_OUT = N803_model_2B(SoluTimes,DoseTimes,Pars,Yic,Opts) + +% declare optional input default values +fullOut = false ;% if true, more outputs are added to Y_OUT (line 303) +timeOut = 30 ;% scalar time limit [seconds] for ODE solver + % if exceeded, an error is thrown +warnOff = true ;% if true, warning messages are suppressed +AbsTol = 1e-6 ;% absolute error tolerance (see 'odeset') +RelTol = 1e-3 ;% relative error tolerance (see 'odeset') +odeSolver = '15s' ;% ode solver function: 'ode15s' or 'ode23s' + +% overwrite optional inputs, if specified +for n = 1:length(Opts) + if strcmp('fullOut', Opts(n)) ; fullOut = Opts{n+1} ; end + if strcmp('timeOut', Opts(n)) ; timeOut = Opts{n+1} ; end + if strcmp('warnOff', Opts(n)) ; warnOff = Opts{n+1} ; end + if strcmp('AbsTol', Opts(n)) ; AbsTol = Opts{n+1} ; end + if strcmp('RelTol', Opts(n)) ; RelTol = Opts{n+1} ; end + if strcmp('odeSolver', Opts(n)) ; odeSolver = Opts{n+1} ; end +end + +% check input vector orientation +if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end +if size(DoseTimes ,1)==1 ; DoseTimes = DoseTimes' ; end +if size(Yic ,2)==1 ; Yic = Yic' ; end + +% check that members of 'DoseTimes' are members of 'SoluTimes' +if length( unique( [ SoluTimes ; DoseTimes ] ) ) > length( SoluTimes ) + error('All members of DoseTimes must be members of SoluTimes!') +end + +% declare index in 'SoluTimes' for start, dose, and end times +Pieces = [ 1 ; find(ismember(SoluTimes,DoseTimes)) ; length(SoluTimes) ] ; + +% suppress warnings (if requested) +if warnOff ; warning off ; end + +%% ------------------------------------------------------------------------ +% Rename inputed parameters ----------------------------------------------- + +q = Pars(01) ;% V growth rate constant [/d] +g = Pars(02) ;% E killing rate constant [uL/#-d] +V50E = Pars(03) ;% 50% V saturation for E activation [#/mL] +V50R = Pars(04) ;% 50% V saturation for R generation [#/mL] +a = Pars(05) ;% E activation rate constant [/d] +m = Pars(06) ;% E reversion rate constant [/d] +E50 = Pars(07) ;% 50% E saturation for proliferation [#/uL] +p = Pars(08) ;% E resting proliferation rate constant [/d] +pA = Pars(09) ;% E active proliferation rate constant [/d] +d = Pars(10) ;% E resting death rate constant [/d] +dA = Pars(11) ;% E active death rate constant [/d] +Xi = Pars(12) ;% X initial condition (at dose) [pmol/kg] +ka = Pars(13) ;% N803 absorption rate constant [/d] +ke = Pars(14) ;% N803 elimination rate constant [/d] +vd = Pars(15) ;% N803 'volume of distribution'/'bioavailability' [L/kg] +C50 = Pars(16) ;% 50% N803 saturation (drug effects) [pM] +ro = Pars(17) ;% E proliferation stimulation factor [] +al = Pars(18) ;% E activation stimulation factor [] +si = Pars(19) ;% R generation stimulation factor [] +s = Pars(20) ;% R generation rate constant [/d] +dR = Pars(21) ;% R decay rate constant [/d] +la = Pars(22) ;% E killing regulation factor [] +ph = Pars(23) ;% E proliferation regulation factor [] +ze = Pars(24) ;% E activation regulation factor [] + +%% ------------------------------------------------------------------------ +% Solve for each solution period ------------------------------------------ + +Y = zeros(Pieces(end),length(Yic)) ;% percolate output matrix + +odeOptions = odeset('AbsTol',AbsTol,'RelTol',RelTol,'Jacobian',@jacob) ; + +for n = 1:length(Pieces)-1 + + solvertime = tic ;% timer (for catching a stalled execution) + + Piece = Pieces(n):Pieces(n+1) ;% index for solution piece + + switch odeSolver % solve ODE for time window indexed by Piece + case '15s' + [~,Y_TEMP] = ode15s(@system,SoluTimes(Piece),Yic,odeOptions) ; + case '23s' + [~,Y_TEMP] = ode23s(@system,SoluTimes(Piece),Yic,odeOptions) ; + end + + if length(Piece) == 2 % exclude timepoints added by ODE solver + Y_TEMP = [ Y_TEMP(1,:) ; Y_TEMP(end,:) ] ; + end + + Y(Piece,:) = Y_TEMP ; % add 'Y_TEMP' to full solution 'Y' + + Yic = Y_TEMP(end,:) ;% update ICs for next solution piece + Yic(07) = Yic(07) + Xi ;% add drug dose as initial condition +end + +% Extremely small values can create numerical artifacts where the value +% becomes negative (mathematically impossible). This occurs for 'bad' +% parameter sets during calibration, and it occurs for the drug (where +% extremely small values are irrelevant). Taking the absolute value of Y +% prevents future errors with the 'log10' function. +Y = abs(Y) ; + +Y_OUT = [ log10( Y(:,1) ) ... SIV virions [log(#/mL)] + , sum( Y(:,2:6), 2) ];% total CD8+ T cells [#/uL] + +%% ------------------------------------------------------------------------ +% Add columns to 'Y_OUT' (if requested) ----------------------------------- + +if fullOut + Y_OUT = outputs(Y_OUT,Y) ; +end + +%% ======================================================================== +% 'system' nested function +% ======================================================================== + function dY = system(~,Y) + + if toc(solvertime) >= timeOut + error('ODE solver stalled.') + end + + V = Y(01) ;% SIV virions [#/mL] + E0 = Y(02) ;% resting CD8+ T cells [#/uL] + EA = Y(03:06) ;% active CD8+ T cells [#/uL] (E1 to E4) + X = Y(07) ;% N803 at absorption site [pmol/kg] + C = Y(08) ;% N803 plasma concentration [pM] + R = Y(09:10) ;% immune regulation [] + F = terms(Y') ;% see 'terms' function + + P = F(11) ;% E resting proliferation rate [/d] + A = F(12) ;% E activation rate [/d] + G = F(13) ;% E killing rate [#/uL-d] + S = F(14) ;% R generation rate [/d] + + dY = Y ;% dY/dt + + % calculate rates of change for V,E,X,C,R + dY(01) = q*V - G*sum(EA)*V ;% dV/dt + dY(02) = P*E0 - d *E0 + m *EA(4) - A *E0 ;% dE0/dt + dY(03) = - dA*EA(1) + 2*A *E0 - pA*EA(1) ;% dE1/dt + dY(4:5)= - dA*EA(2:3) + 2*pA*EA(1:2) - pA*EA(2:3) ;% dE2-3/dt + dY(06) = - dA*EA(4) + 2*pA*EA(3) - m *EA(4) ;% dE4/dt + dY(07) = - ka*X ;% dX/dt + dY(08) = ka*X/vd - ke*C ;% dC/dt + dY(09) = S - dR*R(1) ;% dR1/dt + dY(10) = dR*R(1) - dR*R(2) ;% dR2/dt + end +%% ======================================================================== +% 'jacob' nested function +% ======================================================================== + function J = jacob(~,Y) + + V = Y(01) ;% SIV virions [#/mL] + E0 = Y(02) ;% resting CD8+ T cells [#/uL] + EA = sum(Y(03:06)) ;% active CD8+ T cells [#/uL] + C = Y(08) ;% N803 plasma concentration [pM] + F = terms(Y') ;% see 'terms' function + + P = F(11) ;% E proliferation rate [/d] + A = F(12) ;% E activation rate [/d] + G = F(13) ;% E killing rate [#/uL-d] + + J = zeros(10) ; + + % Partial derivatives of P,A,G,S w.r.t. V,E,C,R + dPdR = - P * ph * F(05) ; + dAdR = - A * ze * F(06) ; + dGdR = - G * la * F(07) ; + dPdE = - P * F(08) / E50 ; + dAdV = A * V50E / V / (V50E+V) ; + dSdV = s * V50R / (V50R+V)^2 ; + dOdC = C50 / (C50+C)^2 ; + dPdC = p * F(08) * ro * dOdC * F(05) ; + dAdC = a * F(09) * al * dOdC * F(06) ; + dSdC = s * si * dOdC ; + + % Partial derivatives of dV/dt + J(01,01 ) = q - G*EA ;% wrt V + J(01,03:06) = - G *V ;% wrt E + J(01,10 ) = - dGdR*EA*V ;% wrt R + + % Partial derivatives of dE0/dt + J(02,01 ) = - dAdV*E0 ;% wrt V + J(02,02 ) = dPdE*E0 + P - d - A ;% wrt E0 + J(02,03:05) = dPdE*E0 ;% wrt E1-E3 + J(02,06 ) = dPdE*E0 + m ;% wrt E4 + J(02,08 ) = dPdC*E0 - dAdC*E0 ;% wrt C + J(02,10 ) = dPdR*E0 - dAdR*E0 ;% wrt R + + % Partial derivatives of dEA/dt + J(03,01 ) = 2*dAdV * E0 ;% dE1/dt wrt V + J(03,02 ) = 2*A ;% dE1/dt wrt E1 + J(03,08 ) = 2*dAdC * E0 ;% dE1/dt wrt C + J(03,10 ) = 2*dAdR * E0 ;% dE1/dt wrt R + J(23:11:45) = - dA - pA ;% dEi/dt wrt Ei for i=1-3 + J(24:11:46) = 2*pA ;% dEi/dt wrt Ei-1 for i=2-4 + J(06,06 ) = - dA - m ;% dE8/dt wrt E4 + + % Partial derivatives of dX/dt + J(07,07) = -ka ;% wrt X + + % Partial derivatives of dC/dt + J(08,07) = ka/vd ;% wrt X + J(08,08) = -ke ;% wrt C + + % Partial derivatives of dR/dt + J(09,01) = dSdV ;% dR1/dt wrt V + J(09,08) = dSdC ;% dR1/dt wrt C + J(09,09) = -dR ;% dR1/dt wrt R1 + J(10,09) = dR ;% dR2/dt wrt R1 + J(10,10) = -dR ;% dR2/dt wrt R2 + end +%% ======================================================================== +% 'terms' nested function +% ======================================================================== + function F = terms(Y) + + V = Y(:,01) ;% SIV virions [#/mL] + E = sum(Y(:,02:06),2) ;% CD8+ T cells [#/uL] + C = Y(:,08) ;% N803 plasma concentration [pM] + R = Y(:,10) ;% immune regulation [] + + % percolate vector for storing terms + F = zeros(size(Y,1),14) ; + + F(:,01) = C./ (C50+C) ;% N803 effect [] + + F(:,02) = 1 + ro*F(:,1) ;% E proliferation stimulation [] + F(:,03) = 1 + al*F(:,1) ;% E activation stimulation [] + F(:,04) = si*F(:,1) ;% N803-induced R generation [] + + F(:,05) = 1./ (1+ph*R) ;% E proliferation regulation [] + F(:,06) = 1./ (1+ze*R) ;% E activation regulation [] + F(:,07) = 1./ (1+la*R) ;% E killing regulation [] + + F(:,08) = E50./(E50+E) ;% E proliferation density dependence [] + F(:,09) = V./ (V50E+V) ;% E activation antigen dependence [] + F(:,10) = V./ (V50R+V) ;% R generation antigen dependence [] + + F(:,11) = p* F(:,08).* F(:,02).* F(:,05) ;% E proliferation [/d] + F(:,12) = a* F(:,09).* F(:,03).* F(:,06) ;% E activation rate [/d] + F(:,13) = g* F(:,07) ;% E killing rate [#/uL-d] + F(:,14) = s* (F(:,10)+F(:,04)) ;% R generation rate [/d] + end +%% ======================================================================== +% 'outputs' nested function +% ======================================================================== + function Y_OUT = outputs(Y_OUT,Y) + + % NOTE: There are some zero elements of 'Y_OUT' in order to keep + % the indices matched with 'N803_model_2.m'. + + Y_OUT = [ Y_OUT, Y, zeros( size(Y,1) , 47) ] ; + + V = Y(:,01) ;% SIV virions [#/mL] + E0 = Y(:,02) ;% resting CD8+ T cells [#/uL] + EA = sum( Y(:,3:6) ,2) ;% active CD8+ T cells [#/uL] + Eaf = Y(:,6) ;% last generation of EA + Eap = EA - Eaf ;% EA that is proliferating + + F = terms(Y) ;% see 'terms' function + P = F(:,11) ;% E resting proliferation rate [/d] + A = F(:,12) ;% E activation rate [/d] + G = F(:,13) ;% E killing rate [#/uL-d] + S = F(:,14) ;% R generation rate [/d] + + Y_OUT(:,19) = F(:,01) ;% N803 effect [] + Y_OUT(:,20) = E0+EA ;% total CD8+ T cells [#/uL] + Y_OUT(:,21) = E0+EA ;% total E + Y_OUT(:,23) = EA ;% active E + Y_OUT(:,26) = EA./(E0+EA) ;% activated frequency in E + + Y_OUT(:,28) = P.*E0 ;% E0 proliferation term [#/uL/d] + Y_OUT(:,30) = d.*E0 ;% E0 death term [#/uL/d] + Y_OUT(:,32) = P ;% E0 proliferation rate [/d] + Y_OUT(:,34) = F(:,02) ;% C multiplier for P + Y_OUT(:,35) = F(:,05) ;% R multiplier for P + Y_OUT(:,36) = F(:,08) ;% density-dep multiplier for P + + Y_OUT(:,38) = A.*E0 ;% E0 activation term [#/uL/d] + Y_OUT(:,40) = pA*Eap ;% EA proliferation term [#/uL/d] + Y_OUT(:,41) = dA*EA ;% EA death term [#/uL/d] + Y_OUT(:,43) = m*Eaf ;% EA reversion term [#/uL/d] + Y_OUT(:,45) = A ;% E0 activation rate [/d] + Y_OUT(:,47) = F(:,03) ;% C multiplier for A + Y_OUT(:,49) = F(:,06) ;% R multiplier for A + Y_OUT(:,51) = F(:,09) ;% V multiplier for A + + Y_OUT(:,53) = q*V ;% V growth term [#/mL/d] + Y_OUT(:,54) = G.*EA.*V ;% E killing term [#/mL/d] + Y_OUT(:,55) = G.*EA ;% E killing rate [/d] + Y_OUT(:,56) = F(:,07) ;% R multiplier for G + + Y_OUT(:,57) = S ;% R generation rate [/d] + Y_OUT(:,59) = s*si*F(:,01)./S ;% R gen by C fraction + end +end \ No newline at end of file diff --git a/N803_model_2C.m b/N803_model_2C.m new file mode 100644 index 0000000..524cf98 --- /dev/null +++ b/N803_model_2C.m @@ -0,0 +1,353 @@ +%% N803_model_2C.m - solves ODE model for N-803 treatment of SIV +% +% /--------------------------------------------------------------\ +% | Date: 05/06/2023 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Model 2C is a version of model 2 where CD8+ T cells are represented by +% one resting and one active group. Proliferation of both groups is +% promoted by N803, limited by cell density, and inhibited by regulation. +% Also, regulation is generated with independent parameters. +% +% Nomenclature: V = SIV virions [#/mL] +% E0 = resting CD8+ T cells [#/uL] +% EA = active CD8+ T cells [#/uL] +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = immune regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% +% SoluTimes = ascending vector of days at which to evaluate solution +% +% DoseTimes = ascending vector of days at which to administer doses +% ( members of 'DoseTimes' must be members of 'SoluTimes' ) +% +% Pars = vector of parameters (see list on line 93) +% +% Yic = vector of initial conditions (see list on line 176) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% +% Opts = cell vector used to add optional inputs (see list on line 56) +% EX: Opts = {'AbsTol',1e-2} will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% +% Y_OUT(:,1) = SIV virions [log(#/mL)] at points in 'SoluTimes' +% Y_OUT(:,2) = total CD8+ T cells [#/uL] at points in 'SoluTimes' +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function Y_OUT = N803_model_2C(SoluTimes,DoseTimes,Pars,Yic,Opts) + +% declare optional input default values +fullOut = false ;% if true, more outputs are added to Y_OUT (line 310) +timeOut = 30 ;% scalar time limit [seconds] for ODE solver + % if exceeded, an error is thrown +warnOff = true ;% if true, warning messages are suppressed +AbsTol = 1e-6 ;% absolute error tolerance (see 'odeset') +RelTol = 1e-3 ;% relative error tolerance (see 'odeset') +odeSolver = '15s' ;% ode solver function: 'ode15s' or 'ode23s' + +% overwrite optional inputs, if specified +for n = 1:length(Opts) + if strcmp('fullOut', Opts(n)) ; fullOut = Opts{n+1} ; end + if strcmp('timeOut', Opts(n)) ; timeOut = Opts{n+1} ; end + if strcmp('warnOff', Opts(n)) ; warnOff = Opts{n+1} ; end + if strcmp('AbsTol', Opts(n)) ; AbsTol = Opts{n+1} ; end + if strcmp('RelTol', Opts(n)) ; RelTol = Opts{n+1} ; end + if strcmp('odeSolver', Opts(n)) ; odeSolver = Opts{n+1} ; end +end + +% check input vector orientation +if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end +if size(DoseTimes ,1)==1 ; DoseTimes = DoseTimes' ; end +if size(Yic ,2)==1 ; Yic = Yic' ; end + +% check that members of 'DoseTimes' are members of 'SoluTimes' +if length( unique( [ SoluTimes ; DoseTimes ] ) ) > length( SoluTimes ) + error('All members of DoseTimes must be members of SoluTimes!') +end + +% declare index in 'SoluTimes' for start, dose, and end times +Pieces = [ 1 ; find(ismember(SoluTimes,DoseTimes)) ; length(SoluTimes) ] ; + +% suppress warnings (if requested) +if warnOff ; warning off ; end + +%% ------------------------------------------------------------------------ +% Rename inputed parameters ----------------------------------------------- + +q = Pars(01) ;% V growth rate constant [/d] +g = Pars(02) ;% E killing rate constant [uL/#-d] +V50E = Pars(03) ;% 50% V saturation for E activation [#/mL] +V50R = Pars(04) ;% 50% V saturation for R generation [#/mL] +a = Pars(05) ;% E activation rate constant [/d] +m = Pars(06) ;% E reversion rate constant [/d] +E50 = Pars(07) ;% 50% E saturation for proliferation [#/uL] +p = Pars(08) ;% E0 proliferation rate constant [/d] +pA = Pars(09) ;% EA proliferation rate constant [/d] +d = Pars(10) ;% E0 death rate constant [/d] +dA = Pars(11) ;% EA death rate constant [/d] +Xi = Pars(12) ;% X initial condition (at dose) [pmol/kg] +ka = Pars(13) ;% N803 absorption rate constant [/d] +ke = Pars(14) ;% N803 elimination rate constant [/d] +vd = Pars(15) ;% N803 'volume of distribution'/'bioavailability' [L/kg] +C50 = Pars(16) ;% 50% N803 saturation (drug effects) [pM] +ro = Pars(17) ;% E0 proliferation stimulation factor [] +roA = Pars(18) ;% EA proliferation stimulation factor [] +al = Pars(19) ;% E activation stimulation factor [] +si = Pars(20) ;% R generation stimulation factor [] +s = Pars(21) ;% R generation rate constant [/d] +dR = Pars(22) ;% R decay rate constant [/d] +la = Pars(23) ;% E killing regulation factor [] +ph = Pars(24) ;% E0 proliferation regulation factor [] +phA = Pars(25) ;% EA proliferation regulation factor [] +ze = Pars(26) ;% E activation regulation factor [] + +%% ------------------------------------------------------------------------ +% Solve for each solution period ------------------------------------------ + +Y = zeros(Pieces(end),length(Yic)) ;% percolate output matrix + +odeOptions = odeset('AbsTol',AbsTol,'RelTol',RelTol,'Jacobian',@jacob) ; + +for n = 1:length(Pieces)-1 + + solvertime = tic ;% timer (for catching a stalled execution) + + Piece = Pieces(n):Pieces(n+1) ;% index for solution piece + + switch odeSolver % solve ODE for time window indexed by Piece + case '15s' + [~,Y_TEMP] = ode15s(@system,SoluTimes(Piece),Yic,odeOptions) ; + case '23s' + [~,Y_TEMP] = ode23s(@system,SoluTimes(Piece),Yic,odeOptions) ; + end + + if length(Piece) == 2 % exclude timepoints added by ODE solver + Y_TEMP = [ Y_TEMP(1,:) ; Y_TEMP(end,:) ] ; + end + + Y(Piece,:) = Y_TEMP ; % add 'Y_TEMP' to full solution 'Y' + + Yic = Y_TEMP(end,:) ;% update ICs for next solution piece + Yic(04) = Yic(04) + Xi ;% add drug dose as initial condition +end + +% Extremely small values can create numerical artifacts where the value +% becomes negative (mathematically impossible). This occurs for 'bad' +% parameter sets during calibration, and it occurs for the drug (where +% extremely small values are irrelevant). Taking the absolute value of Y +% prevents future errors with the 'log10' function. +Y = abs(Y) ; + +Y_OUT = [ log10( Y(:,1) ) ... SIV virions [log(#/mL)] + , Y(:,2)+Y(:,3) ];% total CD8+ T cells [#/uL] + +%% ------------------------------------------------------------------------ +% Add columns to 'Y_OUT' (if requested) ----------------------------------- + +if fullOut + Y_OUT = outputs(Y_OUT,Y) ; +end + +%% ======================================================================== +% 'system' nested function +% ======================================================================== + function dY = system(~,Y) + + if toc(solvertime) >= timeOut + error('ODE solver stalled.') + end + + V = Y(01) ;% SIV virions [#/mL] + E0 = Y(02) ;% resting CD8+ T cells [#/uL] + EA = Y(03) ;% active CD8+ T cells [#/uL] + X = Y(04) ;% N803 at absorption site [pmol/kg] + C = Y(05) ;% N803 plasma concentration [pM] + R = Y(06:07) ;% immune regulation [] + F = terms(Y') ;% see 'terms' function + + P0 = F(13) ;% E0 proliferation rate [/d] + PA = F(14) ;% EA proliferation rate [/d] + A = F(15) ;% E activation rate [/d] + G = F(16) ;% E killing rate [#/uL-d] + S = F(17) ;% R generation rate [/d] + + dY = Y ;% dY/dt + + % calculate rates of change for V,E,X,C,R + dY(01) = q*V - G*EA*V ;% dV/dt + dY(02) = P0*E0 - d *E0 - A*E0 + m*EA ;% dE0/dt + dY(03) = PA*EA - dA*EA + 2*A*E0 - m*EA ;% dE2/dt + dY(04) = - ka*X ;% dX/dt + dY(05) = ka*X/vd - ke*C ;% dC/dt + dY(06) = S - dR*R(1) ;% dR1/dt + dY(07) = dR*R(1) - dR*R(2) ;% dR2/dt + end +%% ======================================================================== +% 'jacob' nested function +% ======================================================================== + function J = jacob(~,Y) + + V = Y(01) ;% SIV virions [#/mL] + E0 = Y(02) ;% resting CD8+ T cells [#/uL] + EA = Y(03) ;% active CD8+ T cells [#/uL] + C = Y(05) ;% N803 plasma concentration [pM] + F = terms(Y') ;% modifiers for 'p0','a0','s0','g0' + + P0 = F(13) ;% E0 proliferation rate [/d] + PA = F(14) ;% EA proliferation rate [/d] + A = F(15) ;% E activation rate [/d] + G = F(16) ;% E killing rate [#/uL-d] + + J = zeros(07) ; + + % Partial derivatives of P,A,G,S w.r.t. V,E,C,R + dP0dR = - P0 * ph * F(06) ; + dPAdR = - PA * phA * F(07) ; + dAdR = - A * ze * F(08) ; + dGdR = - G * la * F(09) ; + dP0dE = - P0 * F(10) / E50 ; + dPAdE = - PA * F(10) / E50 ; + dAdV = A * V50E / V / (V50E+V) ; + dSdV = s * V50R / (V50R+V)^2 ; + dOdC = C50 / (C50+C)^2 ; + dP0dC = p * ro * dOdC * F(06) * F(10) ; + dPAdC = pA * roA * dOdC * F(07) * F(10) ; + dAdC = a * al * dOdC * F(08) ; + dSdC = s * si * dOdC ; + + % Partial derivatives of dV/dt + J(01,01) = q - G*EA ;% wrt V + J(01,03) = - G *V ;% wrt E + J(01,07) = - dGdR*EA*V ;% wrt R + + % Partial derivatives of dE0/dt + J(02,01) = - dAdV*E0 ;% wrt V + J(02,02) = dP0dE*E0 + P0 - d - A ;% wrt E0 + J(02,03) = dP0dE*E0 + m ;% wrt EA + J(02,05) = dP0dC*E0 - dAdC*E0 ;% wrt C + J(02,07) = dP0dR*E0 - dAdR*E0 ;% wrt R + + % Partial derivatives of dEA/dt + J(03,01) = 2*dAdV*E0 ;% dEA/dt wrt V + J(03,02) = dPAdE*EA + 2*A ;% dEA/dt wrt E0 + J(03,03) = dPAdE*EA + PA - dA - m ;% dEA/dt wrt EA + J(03,05) = dPAdC*EA + 2*dAdC*E0 ;% dEA/dt wrt C + J(03,07) = dPAdR*EA + 2*dAdR*E0 ;% dEA/dt wrt R + + % Partial derivatives of dX/dt + J(04,04) = -ka ;% wrt X + + % Partial derivatives of dC/dt + J(05,04) = ka/vd ;% wrt X + J(05,05) = -ke ;% wrt C + + % Partial derivatives of dR/dt + J(06,01) = dSdV ;% dR1/dt wrt V + J(06,05) = dSdC ;% dR1/dt wrt C + J(06,06) = -dR ;% dR1/dt wrt R1 + J(07,06) = dR ;% dR2/dt wrt R1 + J(07,07) = -dR ;% dR2/dt wrt R2 + end +%% ======================================================================== +% 'terms' nested function +% ======================================================================== + function F = terms(Y) + + V = Y(:,01) ;% SIV virions [#/mL] + E = Y(:,02)+Y(:,03) ;% CD8+ T cells [#/uL] + C = Y(:,05) ;% N803 plasma concentration [pM] + R = Y(:,07) ;% immune regulation [] + + % percolate vector for storing terms + F = zeros(size(Y,1),17) ; + + F(:,01) = C./ (C50+C) ;% N803 effect [] + + F(:,02) = 1 + ro *F(:,1) ;% E0 proliferation stimulation [] + F(:,03) = 1 + roA*F(:,1) ;% EA proliferation stimulation [] + F(:,04) = 1 + al *F(:,1) ;% E activation stimulation [] + F(:,05) = si *F(:,1) ;% N803-induced R generation [] + + F(:,06) = 1./ (1+ph *R) ;% E0 proliferation regulation [] + F(:,07) = 1./ (1+phA*R) ;% EA proliferation regulation [] + F(:,08) = 1./ (1+ze *R) ;% E activation regulation [] + F(:,09) = 1./ (1+la *R) ;% E killing regulation [] + + F(:,10) = E50./(E50+E) ;% E proliferation density dependence [] + F(:,11) = V./ (V50E+V) ;% E activation antigen dependence [] + F(:,12) = V./ (V50R+V) ;% R generation antigen dependence [] + + F(:,13) = p* F(:,10).* F(:,02).* F(:,06) ;% E0 prolif rate [/d] + F(:,14) = pA*F(:,10).* F(:,03).* F(:,07) ;% EA prolif modifier [/d] + F(:,15) = a* F(:,11).* F(:,04).* F(:,08) ;% E activation rate [/d] + F(:,16) = g* F(:,09) ;% E killing rate [#/uL-d] + F(:,17) = s* (F(:,12)+F(:,05)) ;% R generation rate [/d] + end +%% ======================================================================== +% 'outputs' nested function +% ======================================================================== + function Y_OUT = outputs(Y_OUT,Y) + + % NOTE: There are some zero elements of 'Y_OUT' in order to keep + % the indices matched with 'N803_model_2.m'. + + Y_OUT = [ Y_OUT, Y, zeros( size(Y,1) , 50) ] ; + + V = Y(:,01) ;% SIV virions [#/mL] + E0 = Y(:,02) ;% resting CD8+ T cells [#/uL] + EA = Y(:,03) ;% active CD8+ T cells [#/uL] + + F = terms(Y) ;% see 'terms' function + P0 = F(:,13) ;% E0 proliferation rate [/d] + PA = F(:,14) ;% EA proliferation rate [/d] + A = F(:,15) ;% E activation rate [/d] + G = F(:,16) ;% E killing rate [#/uL-d] + S = F(:,17) ;% R generation rate [/d] + + Y_OUT(:,19) = F(:,01) ;% N803 effect [] + Y_OUT(:,20) = E0+EA ;% total CD8+ T cells [#/uL] + Y_OUT(:,21) = E0+EA ;% total E + Y_OUT(:,23) = EA ;% active E + Y_OUT(:,26) = EA./(E0+EA) ;% activated frequency in E + + Y_OUT(:,28) = P0.*E0 ;% E0 proliferation term [#/uL/d] + Y_OUT(:,30) = d.*E0 ;% E0 death term [#/uL/d] + Y_OUT(:,32) = P0 ;% E0 proliferation rate [/d] + Y_OUT(:,34) = F(:,02) ;% C multiplier for P0 + Y_OUT(:,35) = F(:,06) ;% R multiplier for P0 + Y_OUT(:,36) = F(:,10) ;% density-dep multiplier for P0 + + Y_OUT(:,38) = A.*E0 ;% E0 activation term [#/uL/d] + Y_OUT(:,40) = PA.*EA ;% EA proliferation term [#/uL/d] + Y_OUT(:,41) = dA*EA ;% EA death term [#/uL/d] + Y_OUT(:,43) = m*EA ;% EA reversion term [#/uL/d] + Y_OUT(:,45) = A ;% E0 activation rate [/d] + Y_OUT(:,47) = F(:,04) ;% C multiplier for A + Y_OUT(:,49) = F(:,08) ;% R multiplier for A + Y_OUT(:,51) = F(:,11) ;% V multiplier for A + + Y_OUT(:,53) = q*V ;% V growth term [#/mL/d] + Y_OUT(:,54) = G.*EA.*V ;% E killing term [#/mL/d] + Y_OUT(:,55) = G.*EA ;% E killing rate [/d] + Y_OUT(:,56) = F(:,09) ;% R multiplier for G + + Y_OUT(:,57) = S ;% R generation rate [/d] + Y_OUT(:,59) = s*si*F(:,01)./S ;% R gen by C fraction + end +end \ No newline at end of file diff --git a/N803_model_2D.m b/N803_model_2D.m new file mode 100644 index 0000000..0ec4e3c --- /dev/null +++ b/N803_model_2D.m @@ -0,0 +1,419 @@ +%% N803_model_2D.m - solves ODE model for N-803 treatment of SIV +% +% /--------------------------------------------------------------\ +% | Date: 05/07/2023 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Model 2D is a version of model 2 that includes infected cell killing by +% non-SIV-specific CD8+ T cells. +% +% Nomenclature: V = SIV virions [#/mL] +% S = SIV-specific CD8+ T cells [#/uL] +% N = non-SIV-specific CD8+ T cells [#/uL] +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = immune regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% +% SoluTimes = ascending vector of days at which to evaluate solution +% +% DoseTimes = ascending vector of days at which to administer doses +% ( members of 'DoseTimes' must be members of 'SoluTimes' ) +% +% Pars = vector of parameters (see list on line 93) +% +% Yic = vector of initial conditions (see list on line 188) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% +% Opts = cell vector used to add optional inputs (see list on line 54) +% EX: Opts = {'AbsTol',1e-2} will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% +% Y_OUT(:,1) = SIV virions [log(#/mL)] at points in 'SoluTimes' +% Y_OUT(:,2) = total CD8+ T cells [#/uL] at points in 'SoluTimes' +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function Y_OUT = N803_model_2D(SoluTimes,DoseTimes,Pars,Yic,Opts) + +% declare optional input default values +fullOut = false ;% if true, more outputs are added to Y_OUT (line 356) +IVdoses = [] ;% vector of IV dose strengths (length of 'DoseTimes') +timeOut = 30 ;% scalar time limit [seconds] for ODE solver + % if exceeded, an error is thrown +warnOff = true ;% if true, warning messages are suppressed +AbsTol = 1e-6 ;% absolute error tolerance (see 'odeset') +RelTol = 1e-3 ;% relative error tolerance (see 'odeset') +odeSolver = '15s' ;% ode solver function: 'ode15s' or 'ode23s' + +% overwrite optional inputs, if specified +for n = 1:length(Opts) + if strcmp('fullOut', Opts(n)) ; fullOut = Opts{n+1} ; end + if strcmp('IVdoses', Opts(n)) ; IVdoses = Opts{n+1} ; end + if strcmp('timeOut', Opts(n)) ; timeOut = Opts{n+1} ; end + if strcmp('warnOff', Opts(n)) ; warnOff = Opts{n+1} ; end + if strcmp('AbsTol', Opts(n)) ; AbsTol = Opts{n+1} ; end + if strcmp('RelTol', Opts(n)) ; RelTol = Opts{n+1} ; end + if strcmp('odeSolver', Opts(n)) ; odeSolver = Opts{n+1} ; end +end + +% check input vector orientation +if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end +if size(DoseTimes ,1)==1 ; DoseTimes = DoseTimes' ; end +if size(Yic ,2)==1 ; Yic = Yic' ; end + +% check that members of 'DoseTimes' are members of 'SoluTimes' +if length( unique( [ SoluTimes ; DoseTimes ] ) ) > length( SoluTimes ) + error('All members of DoseTimes must be members of SoluTimes!') +end + +% declare index in 'SoluTimes' for start, dose, and end times +Pieces = [ 1 ; find(ismember(SoluTimes,DoseTimes)) ; length(SoluTimes) ] ; + +% suppress warnings (if requested) +if warnOff ; warning off ; end + +%% ------------------------------------------------------------------------ +% Rename inputed parameters ----------------------------------------------- + +q = Pars(01) ;% V growth rate constant [/d] +gS = Pars(02) ;% S killing rate constant [uL/#-d] +gN = Pars(03) ;% N killing rate constant [uL/#-d] +V50S = Pars(04) ;% 50% V saturation for S activation [#/mL] +V50N = Pars(05) ;% 50% V saturation for N activation [#/mL] +aS = Pars(06) ;% S activation rate constant [/d] +aN = Pars(07) ;% N activation rate constant [/d] +mS = Pars(08) ;% S reversion rate constant [/d] +mN = Pars(09) ;% N reversion rate constant [/d] +S50 = Pars(10) ;% 50% S saturation for proliferation [#/uL] +N50 = Pars(11) ;% 50% N saturation for proliferation [#/uL] +p = Pars(12) ;% S,N resting proliferation rate constant [/d] +pA = Pars(13) ;% S active proliferation rate constant [/d] +d = Pars(14) ;% S,N resting death rate constant [/d] +dA = Pars(15) ;% S,N active death rate constant [/d] +Xi = Pars(16) ;% X initial condition (at dose) [pmol/kg] +ka = Pars(17) ;% N803 absorption rate constant [/d] +ke = Pars(18) ;% N803 elimination rate constant [/d] +vd = Pars(19) ;% N803 'volume of distribution'/'bioavailability' [L/kg] +C50 = Pars(20) ;% 50% N803 saturation (drug effects) [pM] +ro = Pars(21) ;% S,N proliferation stimulation factor [] +alS = Pars(22) ;% S activation stimulation factor [] +alN = Pars(23) ;% N activation stimulation factor [] +sS = Pars(24) ;% R generation rate constant (S activation) [/d] +sN = Pars(25) ;% R generation rate constant (N activation) [/d] +dR = Pars(26) ;% R decay rate constant [/d] +la = Pars(27) ;% S,N killing regulation factor [] +ph = Pars(28) ;% S,N proliferation regulation factor [] +zeS = Pars(29) ;% S activation regulation factor [] +zeN = Pars(30) ;% N activation regulation factor [] + +%% ------------------------------------------------------------------------ +% Solve for each solution period ------------------------------------------ + +Y = zeros(Pieces(end),length(Yic)) ;% percolate output matrix + +odeOptions = odeset('AbsTol',AbsTol,'RelTol',RelTol,'Jacobian',@jacob) ; + +for n = 1:length(Pieces)-1 + + solvertime = tic ;% timer (for catching a stalled execution) + + Piece = Pieces(n):Pieces(n+1) ;% index for solution piece + + switch odeSolver % solve ODE for time window indexed by Piece + case '15s' + [~,Y_TEMP] = ode15s(@system,SoluTimes(Piece),Yic,odeOptions) ; + case '23s' + [~,Y_TEMP] = ode23s(@system,SoluTimes(Piece),Yic,odeOptions) ; + end + + if length(Piece) == 2 % exclude timepoints added by ODE solver + Y_TEMP = [ Y_TEMP(1,:) ; Y_TEMP(end,:) ] ; + end + + Y(Piece,:) = Y_TEMP ; % add 'Y_TEMP' to full solution 'Y' + + Yic = Y_TEMP(end,:) ;% update ICs for next solution piece + + if isempty(IVdoses) + Yic(13) = Yic(13) + Xi ;% add SC drug dose as initial condition + else + Yic(14) = Yic(14) + IVdoses(n)/vd ;% add IV drug dose + end +end + +%% ------------------------------------------------------------------------ +% Prepare main outputs 'Y_OUT' -------------------------------------------- + +% Extremely small values can create numerical artifacts where the value +% becomes negative (mathematically impossible). This occurs for 'bad' +% parameter sets during calibration, and it occurs for the drug (where +% extremely small values are irrelevant). Taking the absolute value of Y +% prevents future errors with the 'log10' function. +Y = abs(Y) ; + +Y_OUT = [ log10( Y(:,1) ) ... SIV virions [log(#/mL)] + , sum( Y(:,2:12), 2) ];% total CD8+ T cells [#/uL] + +%% ------------------------------------------------------------------------ +% Add columns to 'Y_OUT' (if requested) ----------------------------------- + +if fullOut + Y_OUT = outputs(Y_OUT,Y) ; +end + +%% ======================================================================== +% 'system' nested function +% ======================================================================== + function dY = system(~,Y) + + if toc(solvertime) >= timeOut + error('ODE solver stalled.') + end + + V = Y(01) ;% SIV virions [#/mL] + S0 = Y(02) ;% resting SIV-specific CD8+ T cells [#/uL] + S = Y(03:10) ;% active SIV-specific CD8+ T cells [#/uL] (S1-S8) + N0 = Y(11) ;% resting non-SIV-specific CD8+ T cells [#/uL] + N1 = Y(12) ;% active non-SIV-specific CD8+ T cells [#/uL] + X = Y(13) ;% N803 at absorption site [pmol/kg] + C = Y(14) ;% N803 plasma concentration [pM] + R = Y(15:16) ;% immune regulation [] + F = terms(Y') ;% see 'terms' function + + PS = F(13) ;% S resting proliferation rate [/d] + PN = F(14) ;% N resting proliferation rate [/d] + AS = F(15) ;% S activation rate [/d] + AN = F(16) ;% N activation rate [/d] + GS = F(17) ;% S killing rate [/d] + GN = F(18) ;% N killing rate [/d] + SR = F(19)+F(20) ;% regulation generation rate [/d] + + dY = Y ;% dY/dt + + % calculate rates of change for V,S,N,X,C,R + dY(01) = q*V - GS*V - GN*V ;% dV/dt + dY(02) = PS*S0 - d *S0 + mS*S(8) - AS*S0 ;% dS0/dt + dY(03) = - dA*S(1) + 2*AS*S0 - pA*S(1) ;% dS1/dt + dY(4:9)= - dA*S(2:7) + 2*pA*S(1:6) - pA*S(2:7) ;% dS2-7/dt + dY(10) = - dA*S(8) + 2*pA*S(7) - mS*S(8) ;% dS8/dt + dY(11) = PN*N0 - d *N0 + mN*N1 - AN*N0 ;% dN0/dt + dY(12) = - dA*N1 + 2*AN*N0 - mN*N1 ;% dN1/dt + dY(13) = - ka*X ;% dX/dt + dY(14) = ka*X/vd - ke*C ;% dC/dt + dY(15) = SR - dR*R(1) ;% dR1/dt + dY(16) = dR*R(1) - dR*R(2) ;% dR2/dt + end +%% ======================================================================== +% 'jacob' nested function +% ======================================================================== + function J = jacob(~,Y) + + V = Y(01) ;% SIV virions [#/mL] + S0 = Y(02) ;% resting SIV-specific CD8+ T cells [#/uL] + N0 = Y(11) ;% resting non-SIV-specific CD8+ T cells [#/uL] + C = Y(14) ;% N803 plasma concentration [pM] + F = terms(Y') ;% see 'terms' function + + PS = F(13) ;% S resting proliferation rate [/d] + PN = F(14) ;% N resting proliferation rate [/d] + AS = F(15) ;% S activation rate [/d] + AN = F(16) ;% N activation rate [/d] + GS = F(17) ;% S killing rate [/d] + GN = F(18) ;% S killing rate [/d] + + J = zeros(16) ; + + % Partial derivatives of PS,PN,AS,AN,GS w.r.t. V,S,N,C,R + dPSdR = - PS * ph * F(05) ; + dPNdR = - PN * ph * F(05) ; + dASdR = - AS * zeS * F(06) ; + dANdR = - AN * zeN * F(07) ; + dGSdR = - GS * la * F(08) ; + dGNdR = - GN * la * F(08) ; + dPSdX = - PS * F(09) / S50 ; + dPNdX = - PN * F(10) / N50 ; + dASdV = AS * V50S / V / (V50S+V) ; + dANdV = aN * V50N / (V50N+V)^2 * F(07) ; + dOdC = C50 / (C50+C)^2 ; + dPSdC = p * F(09) * ro * dOdC * F(05) ; + dPNdC = p * F(10) * ro * dOdC * F(05) ; + dASdC = aS * F(11) * alS * dOdC * F(06) ; + dANdC = aN * alN * dOdC * F(07) ; + + % Partial derivatives of dV/dt + J(01,01 ) = q - GS - GN ;% wrt V + J(01,03:10) = - gS * F(08) * V ;% wrt SA + J(01,12) = - gN * F(08) * V ;% wrt NA + J(01,16 ) = - dGSdR * V - dGNdR * V ;% wrt R + + % Partial derivatives of dS0/dt + J(02,01 ) = - dASdV*S0 ;% wrt V + J(02,02 ) = dPSdX*S0 + PS - d - AS ;% wrt S0 + J(02,03:09) = dPSdX*S0 ;% wrt S1-S7 + J(02,10 ) = dPSdX*S0 + mS ;% wrt S8 + J(02,14 ) = dPSdC*S0 - dASdC*S0 ;% wrt C + J(02,16 ) = dPSdR*S0 - dASdR*S0 ;% wrt R + + % Partial derivatives of dS1-8/dt + J(03,01 ) = 2*dASdV*S0 ;% dS1/dt wrt V + J(03,02 ) = 2*AS ;% dS1/dt wrt S0 + J(03,14 ) = 2*dASdC*S0 ;% dS1/dt wrt C + J(03,16 ) = 2*dASdR*S0 ;% dS1/dt wrt R + J(35:17:137) = - dA - pA ;% dSi/dt wrt Si for i=1-7 + J(36:17:138) = 2*pA ;% dSi/dt wrt Si-1 for i=2-8 + J(10,10 ) = - dA - mS ;% dS8/dt wrt S8 + + % Partial derivatives of dN0/dt + J(11,01) = - dANdV*N0 ;% wrt V + J(11,11) = dPNdX*N0 + PN - d - AN ;% wrt N0 + J(11,12) = dPNdX*N0 + mN ;% wrt N1 + J(11,14) = dPNdC*N0 - dANdC*N0 ;% wrt C + J(11,16) = dPNdR*N0 - dANdR*N0 ;% wrt R + + % Partial derivatives of dNA/dt + J(12,01) = 2*dANdV*N0 ;% dN1/dt wrt V + J(12,11) = 2*AN ;% dN1/dt wrt N0 + J(12,12) = - dA - mN ;% dN1/dt wrt N1 + J(12,14) = 2*dANdC*N0 ;% dN1/dt wrt C + J(12,16) = 2*dANdR*N0 ;% dN1/dt wrt R + + % Partial derivatives of dX/dt + J(13,13) = -ka ;% wrt X + + % Partial derivatives of dC/dt + J(14,13) = ka/vd ;% wrt X + J(14,14) = -ke ;% wrt C + + % Partial derivatives of dR/dt + J(15,01) = sS*dASdV/aS/F(06) + sN*dANdV/aN/F(07) ;% dR1/dt wrt V + J(15,14) = sS*dASdC/aS/F(06) + sN*dANdC/aN/F(07) ;% dR1/dt wrt C + J(15,15) = -dR ;% dR1/dt wrt R1 + J(16,15) = dR ;% dR2/dt wrt R1 + J(16,16) = -dR ;% dR2/dt wrt R2 + end +%% ======================================================================== +% 'terms' nested function +% ======================================================================== + function F = terms(Y) + + V = Y(:,01) ;% SIV virions [#/mL] + S = sum(Y(:,02:10),2) ;% SIV-specific CD8+ T cells [#/uL] + N = sum(Y(:,11:12),2) ;% non-SIV-specific CD8+ T cells [#/uL] + SA = sum(Y(:,03:10),2) ;% active S [#/uL] + NA = Y(:,12) ;% active N [#/uL] + C = Y(:,14) ;% N803 plasma concentration [pM] + R = Y(:,16) ;% regulatory inhibition [] + + % percolate vector for storing terms + F = zeros(size(Y,1),19) ; + + F(:,01) = C./ (C50+C) ;% N803 effect [] + + F(:,02) = 1 + ro*F(:,1) ;% S,N proliferation stimulation [] + F(:,03) = 1 + alS*F(:,1) ;% S activation stimulation [] + F(:,04) = alN*F(:,1) ;% N803-induced N activation [] + + F(:,05) = 1./ (1+ ph*R) ;% S,N proliferation regulation [] + F(:,06) = 1./ (1+zeS*R) ;% S activation regulation [] + F(:,07) = 1./ (1+zeN*R) ;% N activation regulation [] + F(:,08) = 1./ (1+la*R) ;% S killing regulation [] + + F(:,09) = S50./(S50+S) ;% S proliferation density dependence [] + F(:,10) = N50./(N50+N) ;% N proliferation density dependence [] + F(:,11) = V./ (V50S+V) ;% S activation viral dependence [] + F(:,12) = V./ (V50N+V) ;% N activation viral dependence [] + + F(:,13) = p * F(:,09).* F(:,02).* F(:,05) ;% S proliferation [/d] + F(:,14) = p * F(:,10).* F(:,02).* F(:,05) ;% N proliferation [/d] + F(:,15) = aS* F(:,11).* F(:,03).* F(:,06) ;% S activation [/d] + F(:,16) = aN* (F(:,12)+F(:,04)).* F(:,07) ;% N activation [/d] + F(:,17) = gS * SA.* F(:,08) ;% S killing rate [#/uL-d] + F(:,18) = gN * NA.* F(:,08) ;% N killing rate [#/uL-d] + + F(:,19) = sS* F(:,11).* F(:,03) ;% R generation by AS [/d] + F(:,20) = sN* (F(:,12) + F(:,04)) ;% R generation by AN [/d] + end +%% ======================================================================== +% 'outputs' nested function +% ======================================================================== + function Y_OUT = outputs(Y_OUT,Y) + + Y_OUT = [ Y_OUT, Y, zeros( size(Y,1) , 41) ] ; + + V = Y(:,01) ;% SIV virions [#/mL] + S0 = Y(:,02) ;% resting specific CD8+ T cells [#/uL] + N0 = Y(:,11) ;% resting non-SIV-specific CD8+ T cells [#/uL] + SA = sum( Y(:,03:10) ,2) ;% active specific CD8+ T cells [#/uL] + NA = Y(:,12) ;% active non-SIV-specific T cells [#/uL] + Saf = Y(:,10) ;% last generation of SA + Sap = SA - Saf ;% SA that is proliferating + + F = terms(Y) ;% see 'terms' function + PS = F(:,13) ;% S resting proliferation rate [/d] + PN = F(:,14) ;% N resting proliferation rate [/d] + AS = F(:,15) ;% S activation rate [/d] + AN = F(:,16) ;% N activation rate [/d] + GS = F(:,17) ;% S killing rate [/d] + + Y_OUT(:,19) = F(:,01) ;% N803 effect [] + Y_OUT(:,20) = (S0+N0+SA+NA) ;% total CD8+ T cells [#/uL] + Y_OUT(:,21) = (S0+SA) ;% total S + Y_OUT(:,22) = (N0+NA) ;% total N + Y_OUT(:,23) = SA ;% active S + Y_OUT(:,24) = NA ;% active N + Y_OUT(:,25) = (S0+SA)./(S0+N0+SA+NA) ;% S frequency in (S+N) + Y_OUT(:,26) = SA./(S0+SA) ;% activated frequency in S + Y_OUT(:,27) = NA./(N0+NA) ;% activated frequency in N + + Y_OUT(:,28) = PS.*S0 ;% S0 proliferation term [#/uL/d] + Y_OUT(:,29) = PN.*N0 ;% N0 proliferation term [#/uL/d] + Y_OUT(:,30) = d.*S0 ;% S0 death term [#/uL/d] + Y_OUT(:,31) = d.*N0 ;% N0 death term [#/uL/d] + Y_OUT(:,32) = PS ;% S0 proliferation rate [/d] + Y_OUT(:,33) = PN ;% N0 proliferation rate [/d] + Y_OUT(:,34) = F(:,02) ;% C multiplier for PS,PN + Y_OUT(:,35) = F(:,05) ;% R multiplier for PS,PN + Y_OUT(:,36) = F(:,09) ;% density-dep multiplier for PS + Y_OUT(:,37) = F(:,10) ;% density-dep multiplier for PN + + Y_OUT(:,38) = AS.*S0 ;% S0 activation term [#/uL/d] + Y_OUT(:,39) = AN.*N0 ;% N0 activation term [#/uL/d] + Y_OUT(:,40) = pA*Sap ;% SA proliferation term [#/uL/d] + Y_OUT(:,41) = dA*SA ;% SA death term [#/uL/d] + Y_OUT(:,42) = dA*NA ;% NA death term [#/uL/d] + Y_OUT(:,43) = mS*Saf ;% SA reversion term [#/uL/d] + Y_OUT(:,44) = mN*NA ;% NA reversion term [#/uL/d] + Y_OUT(:,45) = AS ;% S0 activation rate [/d] + Y_OUT(:,46) = AN ;% N0 activation rate [/d] + Y_OUT(:,47) = F(:,03) ;% C multiplier for AS + Y_OUT(:,48) = ( F(:,04) + F(:,12) )./F(:,12) ;% C multiplier for AN + Y_OUT(:,49) = F(:,06) ;% R multiplier for AS + Y_OUT(:,50) = F(:,07) ;% R multiplier for AN + Y_OUT(:,51) = F(:,11) ;% V multiplier for AS + Y_OUT(:,52) = F(:,12) ;% V multiplier for AN + + Y_OUT(:,53) = q*V ;% V growth term [#/mL/d] + Y_OUT(:,54) = GS.*V ;% S killing term [#/mL/d] + Y_OUT(:,55) = GS ;% S killing rate [/d] + Y_OUT(:,56) = F(:,08) ;% R multiplier for GS + + Y_OUT(:,57) = F(:,18) ;% R gen by AS + Y_OUT(:,58) = F(:,19) ;% R gen by AN + Y_OUT(:,59) = F(:,19) ./ (F(:,18) + F(:,19));% R gen by AN fraction + end +end \ No newline at end of file diff --git a/N803_model_2_naive.m b/N803_model_2_naive.m new file mode 100644 index 0000000..34f44ac --- /dev/null +++ b/N803_model_2_naive.m @@ -0,0 +1,299 @@ +%% N803_model_2_naive.m - solves ODE model for N-803 given to healthy NHP +% +% /--------------------------------------------------------------\ +% | Date: 05/12/2023 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Nomenclature: N = CD8+ T cells [#/uL] +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = immune regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% +% SoluTimes = ascending vector of days at which to evaluate solution +% +% DoseTimes = ascending vector of days at which to administer doses +% ( members of 'DoseTimes' must be members of 'SoluTimes' ) +% +% Pars = vector of parameters (see list on line 84) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% +% Opts = cell vector used to add optional inputs (see list on line 46) +% EX: Opts = {'AbsTol',1e-2} will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% +% Y_OUT(:,1) = total CD8+ T cells [#/uL] at points in 'SoluTimes' +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function Y_OUT = N803_model_2_naive(SoluTimes,DoseTimes,Pars,Opts) + +% declare optional input default values +fullOut = false ;% if true, more outputs are added to Y_OUT (line 270) +IVdoses = [] ;% vector of IV dose strengths (length of 'DoseTimes') +timeOut = 30 ;% scalar time limit [seconds] for ODE solver + % if exceeded, an error is thrown +warnOff = true ;% if true, warning messages are suppressed +AbsTol = 1e-6 ;% absolute error tolerance (see 'odeset') +RelTol = 1e-3 ;% relative error tolerance (see 'odeset') +odeSolver = '15s' ;% ode solver function: 'ode15s' or 'ode23s' + +% overwrite optional inputs, if specified +for n = 1:length(Opts) + if strcmp('fullOut', Opts(n)) ; fullOut = Opts{n+1} ; end + if strcmp('IVdoses', Opts(n)) ; IVdoses = Opts{n+1} ; end + if strcmp('timeOut', Opts(n)) ; timeOut = Opts{n+1} ; end + if strcmp('warnOff', Opts(n)) ; warnOff = Opts{n+1} ; end + if strcmp('AbsTol', Opts(n)) ; AbsTol = Opts{n+1} ; end + if strcmp('RelTol', Opts(n)) ; RelTol = Opts{n+1} ; end + if strcmp('odeSolver', Opts(n)) ; odeSolver = Opts{n+1} ; end +end + +% check input vector orientation +if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end +if size(DoseTimes ,1)==1 ; DoseTimes = DoseTimes' ; end + +% check that members of 'DoseTimes' are members of 'SoluTimes' +if length( unique( [ SoluTimes ; DoseTimes ] ) ) > length( SoluTimes ) + error('All members of DoseTimes must be members of SoluTimes!') +end + +% declare index in 'SoluTimes' for start, dose, and end times +Pieces = [ 1 ; find(ismember(SoluTimes,DoseTimes)) ; length(SoluTimes) ] ; + +% suppress warnings (if requested) +if warnOff ; warning off ; end + +%% ------------------------------------------------------------------------ +% Rename inputed parameters ----------------------------------------------- + +aN = Pars(06) ;% N activation rate constant [/d] +mN = Pars(08) ;% N reversion rate constant [/d] +N50 = Pars(10) ;% 50% N saturation for proliferation [#/uL] +p = Pars(11) ;% N resting proliferation rate constant [/d] +d = Pars(13) ;% N resting death rate constant [/d] +dA = Pars(14) ;% N active death rate constant [/d] +Xi = Pars(15) ;% X initial condition [pmol/kg] +ka = Pars(16) ;% N803 absorption rate constant [/d] +ke = Pars(17) ;% N803 elimination rate constant [/d] +vd = Pars(18) ;% N803 'volume of distribution'/'bioavailability' [L/kg] +C50 = Pars(19) ;% 50% N803 saturation (drug effects) [pM] +ro = Pars(20) ;% N proliferation stimulation factor [] +alN = Pars(22) ;% N activation stimulation factor [] +sN = Pars(24) ;% R generation rate (N activation) [/d] +dR = Pars(25) ;% R decay rate constant [/d] +ph = Pars(27) ;% N proliferation regulation factor [] +zeN = Pars(29) ;% N activation regulation factor [] + +%% ------------------------------------------------------------------------ +% Adjust initial conditions for SIV-naive case ---------------------------- + +Yic = [ N50*( p/d - 1 ) 0 0 0 0 0 ] ; + +%% ------------------------------------------------------------------------ +% Solve for each solution period ------------------------------------------ + +Y = zeros(Pieces(end),length(Yic)) ;% percolate output matrix + +odeOptions = odeset('AbsTol',AbsTol,'RelTol',RelTol,'Jacobian',@jacob) ; + +for n = 1:length(Pieces)-1 + + solvertime = tic ;% timer (for catching a stalled execution) + + Piece = Pieces(n):Pieces(n+1) ;% index for solution piece + + switch odeSolver % solve ODE for time window indexed by Piece + case '15s' + [~,Y_TEMP] = ode15s(@system,SoluTimes(Piece),Yic,odeOptions) ; + case '23s' + [~,Y_TEMP] = ode23s(@system,SoluTimes(Piece),Yic,odeOptions) ; + end + + if length(Piece) == 2 % exclude timepoints added by ODE solver + Y_TEMP = [ Y_TEMP(1,:) ; Y_TEMP(end,:) ] ; + end + + Y(Piece,:) = Y_TEMP ; % add 'Y_TEMP' to full solution 'Y' + + Yic = Y_TEMP(end,:) ;% update ICs for next solution piece + + if isempty(IVdoses) + Yic(3) = Yic(3) + Xi ;% add SC drug dose as initial condition + else + Yic(4) = Yic(4) + IVdoses(n)/vd ;% add IV drug dose + end +end + +%% ------------------------------------------------------------------------ +% Prepare main outputs 'Y_OUT' -------------------------------------------- + +% Extremely small values can create numerical artifacts where the value +% becomes negative (mathematically impossible). This occurs for 'bad' +% parameter sets during calibration, and it occurs for the drug (where +% extremely small values are irrelevant). Taking the absolute value of Y +% prevents future errors with the 'log10' function. +Y = abs(Y) ; + +Y_OUT = sum( Y(:,1:2), 2) ;% total CD8+ T cells [#/無] + +%% ------------------------------------------------------------------------ +% Add columns to 'Y_OUT' (if requested) ----------------------------------- + +if fullOut + Y_OUT = outputs(Y_OUT,Y) ; +end + +%% ======================================================================== +% 'system' nested function +% ======================================================================== + function dY = system(~,Y) + + if toc(solvertime) >= timeOut + error('ODE solver stalled.') + end + + N0 = Y(1) ;% resting CD8+ T cells [#/uL] + NA = Y(2) ;% active CD8+ T cells [#/uL] + X = Y(3) ;% N803 at absorption site [pmol/kg] + C = Y(4) ;% N803 plasma concentration [pM] + R = Y(5:6) ;% immune regulation [] + F = terms(Y') ;% see 'terms' function + + PN = F(7) ;% N resting proliferation rate [/d] + AN = F(8) ;% N activation rate [/d] + SR = F(9) ;% regulation generation rate [/d] + + dY = Y ;% dY/dt + + % calculate rates of change for N,X,C,R + dY(1) = PN*N0 - d *N0 + mN*NA - AN*N0 ;% dN0/dt + dY(2) = - dA*NA + 2*AN*N0 - mN*NA ;% dNA/dt + dY(3) = - ka*X ;% dX/dt + dY(4) = ka*X/vd - ke*C ;% dC/dt + dY(5) = SR - dR*R(1) ;% dR1/dt + dY(6) = dR*R(1) - dR*R(2) ;% dR2/dt + end +%% ======================================================================== +% 'jacob' nested function +% ======================================================================== + function J = jacob(~,Y) + + N0 = Y(1) ;% resting CD8+ T cells [#/uL] + C = Y(4) ;% N803 plasma concentration [pM] + F = terms(Y') ;% see 'terms' function + + PN = F(7) ;% N resting proliferation rate [/d] + AN = F(8) ;% N activation rate [/d] + + J = zeros(6) ; + + % Partial derivatives of PS,PN,AS,AN,GS,GN w.r.t. N,C,R + dPNdR = - PN * ph * F(04) ; + dANdR = - AN * zeN * F(05) ; + dPNdX = - PN * F(06) / N50 ; + dOdC = C50 / (C50+C)^2 ; + dPNdC = p * F(06) * ro * dOdC * F(04) ; + dANdC = aN * alN * dOdC * F(05) ; + + % Partial derivatives of dN0/dt + J(1,1) = dPNdX*N0 + PN - d - AN ;% wrt N0 + J(1,2) = dPNdX*N0 + mN ;% wrt NA + J(1,4) = dPNdC*N0 - dANdC*N0 ;% wrt C + J(1,6) = dPNdR*N0 - dANdR*N0 ;% wrt R + + % Partial derivatives of dNA/dt + J(2,1) = 2*AN ;% dNA/dt wrt N0 + J(2,2) = - dA - mN ;% dNA/dt wrt NA + J(2,4) = 2*dANdC*N0 ;% dNA/dt wrt C + J(2,6) = 2*dANdR*N0 ;% dNA/dt wrt R + + % Partial derivatives of dX/dt + J(3,3) = -ka ;% wrt X + + % Partial derivatives of dC/dt + J(4,3) = ka/vd ;% wrt X + J(4,4) = -ke ;% wrt C + + % Partial derivatives of dR/dt + J(5,4) = sN*alN*dOdC ;% dR1/dt wrt C + J(5,5) = -dR ;% dR1/dt wrt R1 + J(6,5) = dR ;% dR2/dt wrt R1 + J(6,6) = -dR ;% dR2/dt wrt R2 + end +%% ======================================================================== +% 'terms' nested function +% ======================================================================== + function F = terms(Y) + + N = sum(Y(:,1:2),2) ;% CD8+ T cells [#/uL] + C = Y(:,4) ;% N803 plasma concentration [pM] + R = Y(:,6) ;% regulatory inhibition [] + + % percolate vector for storing terms + F = zeros(size(Y,1),9) ; + + F(:,01) = C./ (C50+C) ;% N803 effect [] + + F(:,02) = 1 + ro*F(:,1) ;% N proliferation stimulation [] + F(:,03) = alN*F(:,1) ;% N803-induced N activation [] + + F(:,04) = 1./ (1+ ph*R) ;% N proliferation regulation [] + F(:,05) = 1./ (1+zeN*R) ;% N activation regulation [] + + F(:,06) = N50./(N50+N) ;% N proliferation density dependence [] + + F(:,07) = p * F(:,06).* F(:,02).* F(:,04) ;% N proliferation [/d] + F(:,08) = aN* F(:,03).* F(:,05) ;% N activation [/d] + + F(:,09) = sN* F(:,03) ;% R generation by AN [/d] + end +%% ======================================================================== +% 'outputs' nested function +% ======================================================================== + function Y_OUT = outputs(Y_OUT,Y) + + Y_OUT = [ Y_OUT, Y, zeros( size(Y,1) , 51) ] ; + + N0 = Y(:,1) ;% resting CD8+ T cells [#/uL] + NA = Y(:,2) ;% active CD8+ T cells [#/uL] + + F = terms(Y) ;% see 'terms' function + PN = F(:,7) ;% N resting proliferation rate [/d] + AN = F(:,8) ;% N activation rate [/d] + + Y_OUT(:,19) = F(:,01) ;% N803 effect [] + Y_OUT(:,20) = (N0+NA) ;% CD8+ T cells [#/uL] + Y_OUT(:,24) = NA ;% NA + Y_OUT(:,27) = (NA)./(N0+NA) ;% activated frequency in N + + Y_OUT(:,29) = PN.*N0 ;% N0 proliferation term [#/uL/d] + Y_OUT(:,31) = d.*N0 ;% N0 death term [#/uL/d] + Y_OUT(:,33) = PN ;% N0 proliferation rate [/d] + Y_OUT(:,34) = F(:,02) ;% C multiplier for PS,PN + Y_OUT(:,35) = F(:,04) ;% R multiplier for PS,PN + Y_OUT(:,37) = F(:,06) ;% density-dep multiplier for PN + + Y_OUT(:,39) = AN.*N0 ;% N0 activation term [#/uL/d] + Y_OUT(:,42) = dA*NA ;% NA death term [#/uL/d] + Y_OUT(:,44) = mN*NA ;% NA reversion term [#/uL/d] + Y_OUT(:,46) = AN ;% N0 activation rate [/d] + Y_OUT(:,50) = F(:,05) ;% R multiplier for AN + + Y_OUT(:,58) = F(:,9) ;% R gen by AN + end +end \ No newline at end of file diff --git a/N803_plotter_2.m b/N803_plotter_2.m index ed8b2ea..c6ca6e3 100644 --- a/N803_plotter_2.m +++ b/N803_plotter_2.m @@ -1,7 +1,7 @@ % N803_plotter_2.m - plot outputs from 'N803_model_2.m' addpath Analyzer addpath Data -load('N803_MCMC','C') +load('N803_ForPlots','C') % This script will generate all figure panels using the model outputs saved % by 'collector.m' and using custom plotting function 'gridplot.m'. @@ -14,17 +14,16 @@ close all set(0,'DefaultFigureWindowStyle','normal') -% scenario colors -Color{1} = [ 2 1 0 ]/2.4 ;% 'Fitting to both (C1)' ; -Color{2} = [ 0 1 2 ]/2.4 ;% 'Fitting to both (C2)' ; -Color{3} = [ 1 1 1 ]/5 ;% 'Fitting to Cohort 1 only' ; -Color{4} = [ 1 1 1 ]/5 ;% 'Fitting to Cohort 2 only' ; -Color{5} = [ 1 1 1 ]/5 ;% 'Webb18 Validation (naive)' ; -Color{6} = [ 1 1 1 ]/5 ;% 'Webb18 Validation (SIV)' ; +% Scenario colors +Color{1} = [ 2 1 0 ]/2.4 ;% Shared constants (Cohort1) MCMC +Color{2} = [ 0 1 2 ]/2.4 ;% Shared constants (Cohort2) MCMC +Color([5,17 ]) = { [ 2 1 0 ]/2.4 } ; +Color([6,18:26 ]) = { [ 0 1 2 ]/2.4 } ; +Color([3,4,7:16]) = { [ 1 1 1 ]/5 } ; % Loop through 'Fig' and create each figure % NOTE: 'for' indent is removed -for Fig = 2:.5:14.5 +for Fig = 3:.5:22 clear P L @@ -41,33 +40,50 @@ % Initialize some variables used for setting up 'gridplot.m' inputs DatID = zeros(1,10) ;% used for mapping experimental data to axes OutID = zeros(1,10) ;% used for mapping model outputs to axes -plotMed = false(1) ;% plot median line -plotCI = true(1) ;% plot 95% CI region ( OR plot whole sample ) xSpace = 7 ;% converts days to weeks XTickSpace = .2 ;% inches between X Tick Marks -DataColor = [1 1 1]*.55 ;% color for data points +shade2 = [] ;% shade detection limit for axis (log virus of 2) % Update figure settings for current figure % NOTE: 'switch' indent is removed switch Fig -case {2,2.5} ; c = 1+2*rem(Fig,1) ;% Fitting (fold) ----------------------- +case {3,3.5, 12,12.5, 13,13.5, 14,14.5, 15,15.5, 17,17.5 18,18.5} % ------- + c = 1+2*rem(Fig,1) ;% cohort + cs = [' c' num2str(c)] ;% cohort string - S = c ;% scenario - figName = [ 'Fig2 Fitting c' num2str(c) ] ; - + switch Fig + case { 3, 3.5} ; S = c ; figName = [ 'Fig3 Fitting' cs ] ; + case {12,12.5} ; S = c+[08 4] ; figName = [ 'FigS2 Mod A' cs ] ; + case {13,13.5} ; S = c+[10 4] ; figName = [ 'FigS3 Mod B' cs ] ; + case {14,14.5} ; S = c+[12 4] ; figName = [ 'FigS4 Mod C' cs ] ; + case {15,15.5} ; S = c+[14 4] ; figName = [ 'FigS5 Mod D' cs ] ; + case {17,17.5} ; S = c ; figName = [ 'FigS7 Mean' cs ] ; + case {18,18.5} ; S = c+[06 4] ; figName = [ 'FigS8 Single' cs ] ; + end YLengths = [.1 1.5 1.5] ; - P{1,3}.YLabel.String = {'Virus','fold change'} ; - P{1,2}.YLabel.String = {'{CD8+ T cells}','fold change',' '} ; - [P{1,3}.YTick,P{1,3}.YTickLabel] = tickmaker( -3:1 , [], 1 ) ; - [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( 0:6 ) ; + P{1,3}.YLabel.String = {'Virus','[CEQ/mL]'} ; + P{1,2}.YLabel.String = {'CD8+ T cells','[cells/\muL]'} ; + [P{1,3}.YTick,P{1,3}.YTickLabel] = tickmaker( 1:8 , [], 1 ) ; + [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( 0:500:3500 ) ; OutID(3) = 01 ; DatID(3) = 01 ; OutID(2) = 02 ; DatID(2) = 02 ; + shade2 = 3 ; + load(['N803_DataSet_' num2str(c)]) X = {DataSet(1).X_data, DataSet(3).X_data} ; Y = {DataSet(1).Y_data, DataSet(3).Y_data} ; - [X,Y,Sid,Cid] = prepper(X,Y,{2,[]},{2,[]},{0,0},[1 0]) ; + [X,Y,Sid] = prepper(X,Y,[],[],[],[1 0]) ; + + if any( Fig == [17,17.5] )% load mean data + load('mean_data') + X = { X{2*(c-1)+1} , X{2*(c-1)+2} } ; + Y = { Y{2*(c-1)+1} , Y{2*(c-1)+2} } ; + for o = 1:2 + Sid{o} = ones ( size( X{o} ) ) ; + end + end if c == 1 [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -4:2:16 ) ; @@ -77,10 +93,10 @@ [P{2,1}.XTick,P{2,1}.XTickLabel] = tickmaker( 0:2:8 ) ; end -case {3,3.5} ; c = 1+2*rem(Fig,1) ;% GZB Validation ----------------------- +case {4,4.5} ; c = 1+2*rem(Fig,1) ;% GZB Validation ----------------------- S = c ;% scenario - figName = [ 'Fig3 Validation c' num2str(c) ] ; + figName = [ 'Fig4 Validation c' num2str(c) ] ; Margins([1 3]) = .5 ; AxisProps.XTickLabelRotation = 0 ; % AxisProps.YAxisLocation = 'right' ; @@ -88,7 +104,7 @@ YLengths = [.1 1.5] ; P{1,2}.YTick = 0:.2:1 ; P{1,2}.YTickLabel = {'0%','20%','40%','60%','80%','100%'} ; - OutID(2) = 27 ; DatID(2) = 01 ; + OutID(2) = 3 ; DatID(2) = 01 ; xSpace = 1 ; XTickSpace = .15 ; @@ -102,13 +118,13 @@ X = {DataSet(15).X_data } ; Y = {DataSet(15).Y_data } ; end - [X,Y,Sid,Cid] = prepper(X,Y,{[],[]},{[],[]},{[],[]},[0 0]) ; + [X,Y,Sid] = prepper(X,Y,[],[],[],[0 0]) ; Y{1} = Y{1}/100 ; -case 4 ; c = 1 ;% Killing Weeds (alternate) ------------------------------- +case 5 ; c = 1 ;% Killing Weeds ------------------------------------------- S = [1 2] ;% scenario - figName = 'Fig4 Killing Weeds' ; + figName = 'Fig5 Killing Weeds' ; Margins(1) = .5 ; YLengths = [0.2 .75 .75 1 .75] ; @@ -116,18 +132,18 @@ [P{1,4}.YTick,P{1,4}.YTickLabel] = tickmaker( -1:3 , [], 1 ) ; [P{1,3}.YTick,P{1,3}.YTickLabel] = tickmaker( -3:0 , [], 1 ) ; [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( -2:1 , [], 1 ) ; - OutID(5) = 58 ;% S killing rate (gS*) - OutID(4) = 29 ;% Sa - OutID(3) = 60 ;% R multiplier for gS* - OutID(2) = 4 ;% S0 + OutID(5) = 4 ;% S killing rate (gS*) + OutID(4) = 5 ;% Sa + OutID(3) = 6 ;% R multiplier for gS* + OutID(2) = 7 ;% S0 [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -1:10 ) ; XTickSpace = .15 ; -case 4.5 ; c = 1 ;% Cell Weeds (alternate) -------------------------------- +case 5.5 ; c = 1 ;% Cell Weeds -------------------------------------------- S = [1 2] ;% scenario - figName = 'Fig4 Cell Weeds' ; + figName = 'Fig5 Cell Weeds' ; Margins(1) = .5 ; YLengths = [0.2 [3 3 4 5]*3.25/15 ] ; @@ -135,50 +151,51 @@ [P{1,4}.YTick,P{1,4}.YTickLabel] = tickmaker( -2:2 , [], 1 ) ; [P{1,3}.YTick,P{1,3}.YTickLabel] = tickmaker( -3:0 , [], 1 ) ; [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( 0:3 , [], 1 ) ; - OutID(5) = 39 ;% S0 activation - OutID(4) = 53 ;% V multiplier for aS* - OutID(3) = 51 ;% R multiplier for aS* - OutID(2) = 49 ;% C multiplier for aS* + OutID(5) = 8 ;% S0 activation + OutID(4) = 9 ;% V multiplier for aS* + OutID(3) = 10 ;% R multiplier for aS* + OutID(2) = 11 ;% C multiplier for aS* [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -1:10 ) ; XTickSpace = .15 ; -case {13,13.5} ; c = 1+2*rem(Fig,1) ;% Single Fitting --------------------- - - S = [0 2]+c ;% scenario - figName = [ 'FigS3 Single Fitting c' num2str(c) ] ; - - YLengths = [.1 1.5 1.5] ; - P{1,3}.YLabel.String = {'Virus','fold change'} ; - P{1,2}.YLabel.String = {'{CD8+ T cells}','fold change',' '} ; - [P{1,3}.YTick,P{1,3}.YTickLabel] = tickmaker( -3:1 , [], 1 ) ; - [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( 0:6 ) ; - OutID(3) = 01 ; DatID(3) = 01 ; - OutID(2) = 02 ; DatID(2) = 02 ; +case {6,7,8,9,10,11} ; c = 2 ;% Treatment Exploration --------------------- - load(['N803_DataSet_' num2str(c)]) - X = {DataSet(1).X_data, DataSet(3).X_data} ; - Y = {DataSet(1).Y_data, DataSet(3).Y_data} ; - [X,Y,Sid,Cid] = prepper(X,Y,{2,[]},{2,[]},{0,0},[1 0]) ; - - if c == 1 - [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -4:2:16 ) ; - [P{2,1}.XTick,P{2,1}.XTickLabel] = tickmaker( 32:2:44 ) ; - else - [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -22:6:-4 ) ; - [P{2,1}.XTick,P{2,1}.XTickLabel] = tickmaker( 0:2:8 ) ; + switch Fig + case 6 ; S = 19 ; figName = 'Fig6A short 2-week' ; + case 7 ; S = 20 ; figName = 'Fig6B short 3-week' ; + case 8 ; S = 21 ; figName = 'Fig6C short 4-week' ; + case 9 ; S = 22 ; figName = 'Fig6D long 2-week' ; + case 10 ; S = 23 ; figName = 'Fig6E long 3-week' ; + case 11 ; S = 24 ; figName = 'Fig6F long 4-week' ; end -case {14,14.5} ; c = 1+2*rem(Fig,1) ;% Webb Validation -------------------- + YLengths = [ .1 1.5 ] ; + P{1,2}.YLabel.String = {'Virus','[CEQ/mL]'} ; + [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( 0:7 , [], 1 ) ; + OutID(2) = 1 ; - S = 4+c ;% scenario + Gaps(2) = 0 ; - figName = ['FigS4 Webb ' num2str(c) ] ; + P{1,1}.XTick = tickmaker( 0:2:24 ) ; + if Fig >=9 + P{1,1}.XTick = tickmaker( 0:4:48 ) ; + shade2 = 2 ; + end + if Fig >=10 + [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( -3:2:7 , [], 1 ) ; + end + +case {19,19.5} ; c = 1+2*rem(Fig,1) ;% Webb Validation -------------------- + + S = 2+c ;% scenario + + figName = ['FigS9 Webb ' num2str(c) ] ; YLengths = [.1 1.5 1.25] ; - P{1,2}.YLabel.String = {'CD8+ T cells','per 無 blood'} ; + P{1,2}.YLabel.String = {'CD8+ T cells','[cells/\muL]'} ; [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( 0:500:3000 ) ; - OutID(2) = 21 ; DatID(2) = 3 ; + OutID(2) = 1 ; DatID(2) = 3 ; XTickSpace = .35 ; [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -1:2 ) ; @@ -186,17 +203,51 @@ load('Webb_DataSet') X = {DataSet(1).X_data, DataSet(2).X_data, DataSet(3).X_data} ; Y = {DataSet(1).Y_data, DataSet(2).Y_data, DataSet(3).Y_data} ; - [X,Y,Sid,Cid] = prepper(X,Y,[],[],[],[]) ; + [X,Y,Sid] = prepper(X,Y,[],[],[],[]) ; if c == 2 - P{1,3}.YLabel.String = {'SIV RNA copies','per mL blood'} ; + P{1,3}.YLabel.String = {'Virus','[CEQ/mL]'} ; [P{1,3}.YTick,P{1,3}.YTickLabel] = tickmaker( 0:5 , [], 1 ) ; - OutID(3) = 3 ; DatID(3) = 1 ; - DatID(2) = 2 ; + OutID(3) = 1 ; DatID(3) = 1 ; + OutID(2) = 2 ; DatID(2) = 2 ; XTickSpace = .2 ; [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -1:18 ) ; end +case {20,20.5} ; c = 2 ;% Treatment Exploration Weeds ----------------------- + + switch Fig + case 20 ; S = 25 ; figName = 'FigS10 weeds short 2-week' ; + case 20.5 ; S = 26 ; figName = 'FigS10 weeds short 3-week' ; + end + + YLengths = [0.1 .9 .9 1.2] ; + [P{1,4}.YTick,P{1,4}.YTickLabel] = tickmaker( -3:1 , [], 1 ) ; + [P{1,3}.YTick,P{1,3}.YTickLabel] = tickmaker( -2:2 , [], 1 ) ; + [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( -3:1 , [], 1 ) ; + OutID(4) = 1 ;% Virus + OutID(3) = 2 ;% Sa + OutID(2) = 3 ;% R multiplier for gS* + + P{1,1}.XTick = tickmaker( 0:2:24 ) ; + +case {21,21.5} ; c = 1+2*rem(Fig,1) ;% Long term -------------------------- + + S = 16+c ;% scenario + figName = [ 'FigS11 Longterm c' num2str(c) ] ; + + YLengths = 2 ; + P{1,1}.YLabel.String = {'Virus','[CEQ/mL]'} ; + [P{1,1}.YTick,P{1,1}.YTickLabel] = tickmaker( 2:7 , [], 1 ) ; + OutID(1) = 1 ; + + if c == 1 + P{1,1}.XTick = tickmaker( (0:10:250)+40 ) ; + else + P{1,1}.XTick = tickmaker( (0:10:250)+4 ) ; + end + P{1,1}.XTickLabel = 0:10:250 ; + otherwise continue end @@ -213,7 +264,7 @@ L{1,1}{1}.MarkerSize = 5 ; L{1,1}{1}.MarkerEdgeColor = 'none' ; L{1,1}{1}.MarkerFaceColor = 'k' ; - if any(Fig == [4,4.5]) % for Weeds plots + if any(Fig == [5,5.5]) % for Weeds plots L{1,1}{1}.MarkerFaceColor = Color{S(1)} ; L{1,1}{2} = L{1,1}{1} ; L{1,1}{2}.X = C(S(2)).Dose / xSpace ; @@ -231,20 +282,25 @@ lid = lid + 1 ; L{1,y}{lid}.X = C(s).X / xSpace ; L{1,y}{lid}.Color = Color{s} ; - L{1,y}{lid}.Y = [ C(s).YSummary(o).lo95 C(s).YSummary(o).hi95 ] ; - L{1,y}{lid}.LineWidth = .5 ; - L{1,y}{lid}.FillColor = [ Color{s} .4 ] ; + if isfield(C(s).FigData,'SETS') + L{1,y}{lid}.Y = C(s).FigData(o).SETS ; + L{1,y}{lid}.LineWidth = .75 ; + else + L{1,y}{lid}.Y =[ C(s).FigData(o).lo95 C(s).FigData(o).hi95 ] ; + L{1,y}{lid}.LineWidth = .5 ; + L{1,y}{lid}.FillColor = [ Color{s} .4 ] ; + end end end % Plot data -dataLine.MarkerEdgeColor = DataColor ; -dataLine.MarkerFaceColor = DataColor ; -dataLine.Color = DataColor ; -dataLine.LineWidth = 1 ; +dataLine.LineWidth = .5 ; dataLine.LineStyle = ':' ; -DataMarkers = repmat ( {'o','^','s'} , 1 , 5 ) ; -MarkerSize = repmat ( [ 4 , 4 , 5 ]/2 , 1 , 5 ) ; +dataLine.Marker = 'o' ; +dataLine.MarkerSize = 2 ; +dataLine.MarkerEdgeColor = [1 1 1]/2 ; +dataLine.MarkerFaceColor = [1 1 1]/2 ; +dataLine.Color = [1 1 1]/2 ; for y = 1:length(YLengths) d = DatID(y) ; if d == 0 ; continue ; end @@ -252,11 +308,8 @@ for r = 1:Sid{d}(end) subMods.X = X{d}(Sid{d}==r)/xSpace ; subMods.Y = Y{d}(Sid{d}==r) ; - subMods.CensorID = Cid{d}(Sid{d}==r) ; - subMods.Marker = DataMarkers{r} ; - subMods.MarkerSize = MarkerSize(r) ; L{1,y}{lid+r} = copyfields(subMods,dataLine) ; - if any(Fig == [14,14.5]) + if any(Fig == [19,19.5]) if d == 2 ; L{1,y}{lid+r}.ErrorBars = ... [104,205,0,0,130,157,124,0,170,580] ; elseif d == 3 ; L{1,y}{lid+r}.ErrorBars = ... @@ -270,9 +323,19 @@ end end +% Plot viral detection limit +if ~isempty(shade2) + lid = length( L{1,shade2} ) + 1 ; + L{1,shade2}{lid}.X = [-500;500] ; %#ok<*SAGROW> + L{1,shade2}{lid}.Y = [-10,2;-10,2] ; + L{1,shade2}{lid}.LineWidth = .5 ; + L{1,shade2}{lid}.Color = [ .8 .8 .8 0 ] ; + L{1,shade2}{lid}.FillColor = [ .8 .8 .8 .4 ] ; +end + % Set [x,y] axes lengths bases on number of ticks for x = 1:nXaxes - XLengths(x) = ( length(P{x,1}.XTick) - 1 ) * XTickSpace ; %#ok + XLengths(x) = ( length(P{x,1}.XTick) - 1 ) * XTickSpace ; end % Make and save plots diff --git a/N803_plotter_2_violins.m b/N803_plotter_2_violins.m new file mode 100644 index 0000000..02d3871 --- /dev/null +++ b/N803_plotter_2_violins.m @@ -0,0 +1,160 @@ +% N803_violins.m - makes violin plots of parameter distributions and AIC + +close all +clear all %#ok +addpath Analyzer + +%% AIC ==================================================================== + +nT = 553/4 ;% number of time points (average per variable) +nY = 4 ;% number of variables + +ModLetter = {'','A','B','C','D'} ; + +NLL = zeros(10,5) ; +AIC = NLL ; +nP = zeros(1,5) ; + +for n = 1:5 + load(['N803_shared_2' ModLetter{n} '_Results'],'Results') + nP(n) = length( Results(1).ParameterFit ) ;% number of parameters + + NLL(:,n) = vertcat( Results(1:10).modelCost ) ; + AIC(:,n) = 2*NLL(:,n) + 2*(nP(n)+nY)*nT/(nT-nP(n)-nY-1) ; +end + +%% +Size = [3 3] ;% figure size +figName = 'FigS1 NLL' ; + +f = figure('Name', figName, ... + 'Color', [1 1 1], ... + 'PaperUnits', 'Inches', ... + 'PaperSize', Size, ... + 'Units', 'Inches', ... + 'Position', [0 0 Size]) ; + +parplotter(NLL, ... + 'AxesLims', [0 6 390 450],... + 'AxesPos', [.2 .2 .75 .75],... + 'figHand', f,... + 'Names', {'Eq 1-11','Model A','Model B','Model C','Model D'},... + 'plotVals', 1,... + 'ViolColor', [.5 .5 .5]) ; + +ylabel('NLL') +set(gca,'FontSize',10,'FontName','Arial','LineWidth',1,'XTickLabelRotation',315) +set(gcf,'InvertHardcopy','off') +print(figName,'-dtiff','-r300') + +figName = 'FigS1 AIC' ; + +f = figure('Name', figName, ... + 'Color', [1 1 1], ... + 'PaperUnits', 'Inches', ... + 'PaperSize', Size, ... + 'Units', 'Inches', ... + 'Position', [0 0 Size]) ; + +parplotter(AIC, ... + 'AxesLims', [0 6 860 980],... + 'AxesPos', [.2 .2 .75 .75],... + 'figHand', f,... + 'Names', {'Eq 1-11','Model A','Model B','Model C','Model D'},... + 'plotVals', 1,... + 'ViolColor', [.5 .5 .5]) ; + +ylabel('AIC') +set(gca,'FontSize',10,'FontName','Arial','LineWidth',1,'XTickLabelRotation',315) +set(gcf,'InvertHardcopy','off') +print(figName,'-dtiff','-r300') + +%% Model Constants ======================================================== +load('N803_ForPlots','C') + +% Retrieving boundaries for constants ------------------------------------- +warning off ; load('N803_shared_2_Results','Results') ; warning on +% q vS vN mS mN p pA d dA C50 alS alN dR ph +id1 = [ 7 8 9 10 11 13 14 15 16 21 23 24 25 27 ] ; +id2 = [ 1 3 4 7 8 11 12 13 14 15 17 18 21 23 ] ; +Bounds = NaN*ones(2,25) ; +Bounds(:,id2) = Results(1).setup.GuessBOUNDS(:,id1) ; + +% Plotting shared constants ----------------------------------------------- +PARS = C(1).ParSample(:,[1:14,19:29]) ;% parameter sample +Size = [6.5 3] ;% figure size +figName = 'FigS6 Constants' ; + +f = figure('Name', figName, ... + 'Color', [1 1 1], ... + 'PaperUnits', 'Inches', ... + 'PaperSize', Size, ... + 'Units', 'Inches', ... + 'Position', [0 0 Size]) ; + +Names = {'q','g','V_{50,S}','V_{50,N}','a_S','a_N','m_S','m_N',... + 'S_{50}','N_{50}','p','p_A','d','d_A','C_{50}','\rho',... + '\alpha_S','\alpha_N','s_S','s_N','d_R','\lambda',... + '\phi','\zeta_S','\zeta_N'} ; + +parplotter(PARS, ... + 'AxesLims', [0 26 -3 4],... + 'AxesPos', [.075 .2 .9 .75],... + 'Bounds', Bounds,... + 'Factors', [1 1e3 1e-3 1e-3 ones(1,21)],...g to [nL/#-d] + 'figHand', f,... vS,vN to [#/uL] + 'logScale', 1,... + 'Names', Names,... + 'ViolColor', [.5 .5 .5]) ; + +set(gca,'FontSize',10,'FontName','Arial','LineWidth',1,'XTickLabelRotation',315) +set(gcf,'InvertHardcopy','off') +print(figName,'-dtiff','-r300') + +%% Model Initials ========================================================= + +% Plotting cohort 1 initials ---------------------------------------------- +Bounds = NaN*ones(2,14) ; +Bounds(:,1) = Results(1).setup.GuessBOUNDS(:,1) ; +PARS = C(1).ParSample(:,[30:41,44:45]) ;% parameter sample +Size = [4.5 3] ;% figure size +figName = 'FigS6 Initials' ; + +f = figure('Name', figName, ... + 'Color', [1 1 1], ... + 'PaperUnits', 'Inches', ... + 'PaperSize', Size, ... + 'Units', 'Inches', ... + 'Position', [0 0 Size]) ; + +Names = {'V','S_0','S_1','S_2','S_3','S_4','S_5','S_6','S_7','S_8',... + 'N_0','N_1','R_1','R_2'} ; + +parplotter(PARS, ... + 'AxesLims', [0 15 -3 4],... + 'AxesPos', [.125 .2 .85 .75],... + 'Bounds', Bounds,... + 'Factors', [1e-3 ones(1,13)],...V to [#/uL] + 'figHand', f,... + 'logScale', 1,... + 'Names', Names,... + 'ViolColor', [ 2 1 0 ]/2.4) ; + +set(gca,'FontSize',10,'FontName','Arial','LineWidth',1,'XTickLabelRotation',315) +% Plotting cohort 2 initials ---------------------------------------------- +Bounds(:,1) = Results(1).setup.GuessBOUNDS(:,2) ; +PARS = C(2).ParSample(:,[30:41,44:45]) ;% parameter sample + +parplotter(PARS, ... + 'AxesLims', [0 15 -3 4],... + 'AxesPos', [.125 .2 .85 .75],... + 'Bounds', Bounds,... + 'Factors', [1e-3 ones(1,13)],...V to [#/uL] + 'figHand', f,... + 'logScale', 1,... + 'Names', Names,... + 'ViolColor', [ 0 1 2 ]/2.4) ; + +set(gca,'FontSize',10,'FontName','Arial','LineWidth',1,'XTickLabelRotation',315) +set(gcf,'InvertHardcopy','off') +print(figName,'-dtiff','-r300') diff --git a/N803_shared_2.m b/N803_shared_2.m new file mode 100644 index 0000000..ee38cb3 --- /dev/null +++ b/N803_shared_2.m @@ -0,0 +1,249 @@ +%% N803_shared_2.m - solves model 2 for 2 cohorts with shared parameters +% +% /--------------------------------------------------------------\ +% | Date: 05/06/2023 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Nomenclature: V = SIV virions [#/mL] +% T8 = total CD8+ T cells [#/uL] +% S0 = resting SIV-specific CD8+ T cells [#/uL] +% SA = active SIV-specific CD8+ T cells [#/uL] (S1-S8) +% N0 = resting non-SIV-specific CD8+ T cells [#/uL] +% NA = active non-SIV-specific CD8+ T cells [#/uL] (N1) +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = immune regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% NOTE: The format of 'SoluTimes' input will determine the format of the +% 'Y_OUT' output and determine whether the model will be run twice, from +% two different sets of initials. +% +% SoluTimes = ascending vector of days at which to evaluate solution +% OR +% SoluTimes = 1x4 cell vector, where each element is a vector of days +% corresponding to a data type (see OUTPUTS) +% +% AllPars = vector of parameters (see list in function) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% NOTE: 'varargin' is used to add optional inputs +% EX: N803_shared_2(SoluTimes,AllPars,'Input Name',[Input Value]) +% +% DoseTimes = ascending vector of days at which to administer doses +% DEFAULT: 1x2 cell vector, where each element is a vector of days +% corresponding to doses in Ellis-Connell data +% +% cohort = scalar to run model for just one cohort ('1' or '2') +% NOTE: 'cohort' is ignored if 'SoluTimes' is a cell vector +% +% noNegs = logical scalar +% if true, throw an error if parameters are negative +% +% Additional inputs will also be passed to 'N803_model_2.m' and used to +% define options in the same manner (see function for list) +% EX: N803_shared_2(SoluTimes,DoseTimes,AllPars,'AbsTol',1e-2} +% will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% NOTE: The format of 'SoluTimes' input will determine the format of the +% 'Y_OUT' output and determine whether the model will be run twice, from +% two different sets of initials. +% +% If 'SoluTimes' is a 1x4 cell vector +% Y_OUT{1} = V at points in 'SoluTimes{1}' [log(#/mL)] cohort 1 +% Y_OUT{2} = T8 at points in 'SoluTimes{2}' [#/uL] cohort 1 +% Y_OUT{3} = V at points in 'SoluTimes{3}' [log(#/mL)] cohort 2 +% Y_OUT{4} = T8 at points in 'SoluTimes{4}' [#/uL] cohort 2 +% PARS(1,:) = [ parameters, initials ] for cohort 1 (see code) +% PARS(2,:) = [ parameters, initials ] for cohort 2 (see code) +% +% If 'SoluTimes' is a vector of days +% Y_OUT = model outputs as described in 'N803_model_2.m' +% PARS = [ parameters, initials ] (see code) +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function [Y_OUT,PARS] = N803_shared_2(SoluTimes,AllPars,varargin) + +% Declare optional input default values (then overwrite, if specified) +DoseTimes = { 7*[0:3 5:8 37:40]' ... dose days for cohort 1 + , 7*(0:2:4)' };% dose days for cohort 2 +cohort = 1 ;% run for this cohort (ignored if 'SoluTimes' is a cell) +noNegs = true ;% do not allow negative parameter values +isNaive = false ;% if true, all IC are zero (except for resting N) + +for n = 1:length(varargin) +if strcmp('DoseTimes', varargin(n)) ; DoseTimes = varargin{n+1} ; end +if strcmp('cohort', varargin(n)) ; cohort = varargin{n+1} ; end +if strcmp('noNegs', varargin(n)) ; noNegs = varargin{n+1} ; end +if strcmp('isNaive', varargin(n)) ; isNaive = varargin{n+1} ; end +end + +Y_OUT = SoluTimes ;% initialize response output +PARS = [] ;% initialize parameter output + +% Rename inputed parameters ----------------------------------------------- +Vi(1) = AllPars(01) ;% V initial [#/mL] (cohort 1) +Vi(2) = AllPars(02) ;% V initial [#/mL] (cohort 2) +SNi(1) = AllPars(03) ;% S+N initial [#/uL] (cohort 1) +SNi(2) = AllPars(04) ;% S+N initial [#/uL] (cohort 2) +fS(1) = AllPars(05) ;% initial frequency: S/(S+N) (cohort 1) +fS(2) = AllPars(06) ;% initial frequency: S/(S+N) (cohort 2) +Si = SNi.*fS ;% S initial value [#/uL] [cohort 1,2] +Ni = SNi-Si ;% N initial value [#/uL] [cohort 1,2] + +q = AllPars(07) ;% V growth rate constant [/d] +V50S = AllPars(08) ;% 50% V saturation for S activation [#/mL] +V50N = AllPars(09) ;% 50% V saturation for N activation [#/mL] +mS = AllPars(10) ;% S reversion rate constant [/d] +mN = AllPars(11) ;% N reversion rate constant [/d] + +S50 = AllPars(12)* fS(1) ;% 50% S saturation for proliferation [#/uL] +N50 = AllPars(12)*(1-fS(1)) ;% 50% N saturation for proliferation [#/uL] +p = AllPars(13) ;% S,N resting proliferation rate constant [/d] +pA = AllPars(14) ;% S active proliferation rate constant [/d] +d = AllPars(15) ;% S,N resting death rate constant [/d] +dA = AllPars(16) ;% S,N acting death rate constant [/d] + +Xi = AllPars(17) ;% X initial condition (at dose) [pmol/kg] +ka = AllPars(18) ;% N803 absorption rate constant [/d] +ke = AllPars(19) ;% N803 elimination rate constant [/d] +vd = AllPars(20) ;% N803 'vol. of distribution'/'bioavailability' [L/kg] +C50 = AllPars(21) ;% 50% N803 saturation (drug effects) [pM] +pm = AllPars(22) ;% S,N induced resting expansion rate [] +alS = AllPars(23) ;% S activation stimulation factor [] +alN = AllPars(24) ;% N activation stimulation factor [] + +dR = AllPars(25) ;% R decay rate constant [/d] +sig = AllPars(26) ;% ratio of regulation generation rates sN/sS +ph = AllPars(27) ;% S,N proliferation regulation factor [] + +%% ------------------------------------------------------------------------ +% Calculate some initial conditions & parameters -------------------------- + +YS = Vi ./ (Vi + V50S) ;% initial V/(V50S+V) [cohort 1,2] +YN = Vi ./ (Vi + V50N) ;% initial V/(V50N+V) [cohort 1,2] + +K = YS + sig*YN ;% collection [cohort 1,2] +Ri = [ 1 K(2)/K(1) ] ;% initial regulation [cohort 1,2] +sS = dR/K(1) ;% R generation rate constant (S activation) [/d] +sN = sig*sS ;% R generation rate constant (N activation) [/d] + +ro = pm/p ;% S,N proliferation stimulation factor [] + +% calculate aS,aN,zeS,zeN +PS = p*S50./(S50+Si)./(1+ph*Ri) ;% initial S0 prolif rates [/d] ([c1,c2]) +PN = p*N50./(N50+Ni)./(1+ph*Ri) ;% initial N0 prolif rates [/d] ([c1,c2]) +US = 2*(2*pA/(pA+dA))^7/(mS+dA) ;% initial S8 / (AS*S0) +UN = 2 /(mN+dA) ;% initial N1 / (AN*N0) +AS = (d-PS) / (mS*US-1) ;% initial S0 activation [/d] ([c1,c2]) +AN = (d-PN) / (mN*UN-1) ;% initial N0 activation [/d] ([c1,c2]) +WS = AS(1)/AS(2)*YS(2)/YS(1) ; +WN = AN(1)/AN(2)*YN(2)/YN(1) ; +zeS = (WS-1) / (Ri(2)-WS) ;% S activation regulation factor [] +zeN = (WN-1) / (Ri(2)-WN) ;% N activation regulation factor [] +aS = AS(1)/YS(1) * (1+zeS) ;% S activation rate constant [/d] +aN = AN(1)/YN(1) * (1+zeN) ;% N activation rate constant [/d] + +% calculate initial S0,N0,SA for each cohort +ZS = US ; +for i = 1:7 + ZS = ZS + 2*(2*pA)^(i-1)/(pA+dA)^i ;% initial SA / (AS*S0) +end +S0 = Si./(1+ZS*AS) ;% initial S0 +N0 = Ni./(1+UN*AN) ;% initial N0 +SA = Si - S0 ;% initial SA (S1+S2...+S8) +NA = Ni - N0 ;% initial NA (N1) + +% calculate g and la +la = (SA(2) - SA(1)) / (SA(1)*Ri(2) - SA(2)) ;% S killing reg factor [] +g = q * (1+la) / SA(1) ;% S killing rate constant [uL/#-d] + +%% Do for each cohort (NOT indenting loop) ================================ +for c = 1:2 + +% solve for initial S1-8 +S = zeros(1,8) ;% initial S1-S8 +S(1) = 2*AS(c)*S0(c) / (dA+pA) ;% S1 +for i = 2:7 + S(i) = 2*pA*S(i-1) / (dA+pA) ;% S2 to S7 +end +S(8) = 2*pA*S(7) / (dA+mS) ;% S8 + +%% ------------------------------------------------------------------------ +% Prepare parameter and initial value vectors and call 'N803_model_2' ----- + +% V S0-8 N0 N1 X C R initial values +Yic = [ Vi(c) S0(c) S N0(c) NA(c) 0 0 Ri(c) Ri(c) ] ; + +Pars(01) = q ;% V growth rate constant [/d] +Pars(02) = g ;% S killing rate constant [uL/#-d] +Pars(03) = V50S ;% 50% V saturation for S activation [#/mL] +Pars(04) = V50N ;% 50% V saturation for N activation [#/mL] +Pars(05) = aS ;% S activation rate constant [/d] +Pars(06) = aN ;% N activation rate constant [/d] +Pars(07) = mS ;% S reversion rate constant [/d] +Pars(08) = mN ;% N reversion rate constant [/d] +Pars(09) = S50 ;% 50% S saturation for proliferation [#/uL] +Pars(10) = N50 ;% 50% N saturation for proliferation [#/uL] +Pars(11) = p ;% S,N resting proliferation rate constant [/d] +Pars(12) = pA ;% S active proliferation rate constant [/d] +Pars(13) = d ;% S,N resting death rate constant [/d] +Pars(14) = dA ;% S,N active death rate constant [/d] +Pars(15) = Xi ;% X initial condition (at dose) [pmol/kg] +Pars(16) = ka ;% N803 absorption rate constant [/d] +Pars(17) = ke ;% N803 elimination rate constant [/d] +Pars(18) = vd ;% N803 'volume of distribution'/'bioavailability' [L/kg] +Pars(19) = C50 ;% 50% N803 saturation (drug effects) [pM] +Pars(20) = ro ;% S,N proliferation stimulation factor [] +Pars(21) = alS ;% S activation stimulation factor [] +Pars(22) = alN ;% N activation stimulation factor [] +Pars(23) = sS ;% R generation rate constant (S activation) [/d] +Pars(24) = sN ;% R generation rate constant (N activation) [/d] +Pars(25) = dR ;% R decay rate constant [/d] +Pars(26) = la ;% S killing regulation factor [] +Pars(27) = ph ;% S,N proliferation regulation factor [] +Pars(28) = zeS ;% S activation regulation factor [] +Pars(29) = zeN ;% N activation regulation factor [] + +if any( [ Pars Yic ] < 0 ) && noNegs + error('Negative parameters or initial values.') +end + +if iscell(SoluTimes) % if calibrating to Ellis-Connell data + + VirusTimes = SoluTimes{2*(c-1)+1} ;% times for virus data points + CD8TcTimes = SoluTimes{2*(c-1)+2} ;% times for CD8 T cell data points + + X = sort( unique( [ VirusTimes ; CD8TcTimes ; DoseTimes{c} ] ) ) ; + + Y = N803_model_2( X, DoseTimes{c}, Pars, Yic, varargin ) ; + + Y_OUT{2*(c-1)+1} = interp1( X, Y(:,1), VirusTimes ) ; + Y_OUT{2*(c-1)+2} = interp1( X, Y(:,2), CD8TcTimes ) ; + + PARS = [ PARS ; [ Pars Yic ] ] ; %#ok + +elseif cohort == c % if evaluating for one cohort + + if isNaive + Y_OUT = N803_model_2_naive(SoluTimes, DoseTimes, Pars, varargin ) ; + else + Y_OUT = N803_model_2( SoluTimes, DoseTimes, Pars, Yic, varargin ) ; + PARS = [ Pars Yic ] ; + end +end +end +end \ No newline at end of file diff --git a/N803_shared_2A.m b/N803_shared_2A.m new file mode 100644 index 0000000..479ad21 --- /dev/null +++ b/N803_shared_2A.m @@ -0,0 +1,246 @@ +%% N803_shared_2A.m - solves model 2A for 2 cohorts with shared parameters +% +% /--------------------------------------------------------------\ +% | Date: 05/06/2023 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Model 2A is a version of model 2 where SIV-specific cells divide 6 times +% after activation, instead of 8. +% +% Nomenclature: V = SIV virions [#/mL] +% T8 = total CD8+ T cells [#/uL] +% S0 = resting SIV-specific CD8+ T cells [#/uL] +% SA = active SIV-specific CD8+ T cells [#/uL] (S1-S8) +% N0 = resting non-SIV-specific CD8+ T cells [#/uL] +% NA = active non-SIV-specific CD8+ T cells [#/uL] (N1) +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = immune regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% NOTE: The format of 'SoluTimes' input will determine the format of the +% 'Y_OUT' output and determine whether the model will be run twice, from +% two different sets of initials. +% +% SoluTimes = ascending vector of days at which to evaluate solution +% OR +% SoluTimes = 1x4 cell vector, where each element is a vector of days +% corresponding to a data type (see OUTPUTS) +% +% AllPars = vector of parameters (see list in function) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% NOTE: 'varargin' is used to add optional inputs +% EX: N803_shared_2(SoluTimes,AllPars,'Input Name',[Input Value]) +% +% DoseTimes = ascending vector of days at which to administer doses +% DEFAULT: 1x2 cell vector, where each element is a vector of days +% corresponding to doses in Ellis-Connell data +% +% cohort = scalar to run model for just one cohort ('1' or '2') +% NOTE: 'cohort' is ignored if 'SoluTimes' is a cell vector +% +% noNegs = logical scalar +% if true, throw an error if parameters are negative +% +% Additional inputs will also be passed to 'N803_model_2.m' and used to +% define options in the same manner (see function for list) +% EX: N803_shared_2(SoluTimes,DoseTimes,AllPars,'AbsTol',1e-2} +% will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% NOTE: The format of 'SoluTimes' input will determine the format of the +% 'Y_OUT' output and determine whether the model will be run twice, from +% two different sets of initials. +% +% If 'SoluTimes' is a 1x4 cell vector +% Y_OUT{1} = V at points in 'SoluTimes{1}' [log(#/mL)] cohort 1 +% Y_OUT{2} = T8 at points in 'SoluTimes{2}' [#/uL] cohort 1 +% Y_OUT{3} = V at points in 'SoluTimes{3}' [log(#/mL)] cohort 2 +% Y_OUT{4} = T8 at points in 'SoluTimes{4}' [#/uL] cohort 2 +% PARS(1,:) = [ parameters, initials ] for cohort 1 (see code) +% PARS(2,:) = [ parameters, initials ] for cohort 2 (see code) +% +% If 'SoluTimes' is a vector of days +% Y_OUT = model outputs as described in 'N803_model_2.m' +% PARS = [ parameters, initials ] (see code) +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function [Y_OUT,PARS] = N803_shared_2A(SoluTimes,AllPars,varargin) + +% Declare optional input default values (then overwrite, if specified) +DoseTimes = { 7*[0:3 5:8 37:40]' ... dose days for cohort 1 + , 7*(0:2:4)' };% dose days for cohort 2 +cohort = 1 ;% run for this cohort (ignored if 'SoluTimes' is a cell) +noNegs = true ;% do not allow negative parameter values + +for n = 1:length(varargin) +if strcmp('DoseTimes', varargin(n)) ; DoseTimes = varargin{n+1} ; end +if strcmp('cohort', varargin(n)) ; cohort = varargin{n+1} ; end +if strcmp('noNegs', varargin(n)) ; noNegs = varargin{n+1} ; end +end + +Y_OUT = SoluTimes ;% initialize response output +PARS = [] ;% initialize parameter output + +% Rename inputed parameters ----------------------------------------------- +Vi(1) = AllPars(01) ;% V initial [#/mL] (cohort 1) +Vi(2) = AllPars(02) ;% V initial [#/mL] (cohort 2) +SNi(1) = AllPars(03) ;% S+N initial [#/uL] (cohort 1) +SNi(2) = AllPars(04) ;% S+N initial [#/uL] (cohort 2) +fS(1) = AllPars(05) ;% initial frequency: S/(S+N) (cohort 1) +fS(2) = AllPars(06) ;% initial frequency: S/(S+N) (cohort 2) +Si = SNi.*fS ;% S initial value [#/uL] [cohort 1,2] +Ni = SNi-Si ;% N initial value [#/uL] [cohort 1,2] + +q = AllPars(07) ;% V growth rate constant [/d] +V50S = AllPars(08) ;% 50% V saturation for S activation [#/mL] +V50N = AllPars(09) ;% 50% V saturation for N activation [#/mL] +mS = AllPars(10) ;% S reversion rate constant [/d] +mN = AllPars(11) ;% N reversion rate constant [/d] + +S50 = AllPars(12)* fS(1) ;% 50% S saturation for proliferation [#/uL] +N50 = AllPars(12)*(1-fS(1)) ;% 50% N saturation for proliferation [#/uL] +p = AllPars(13) ;% S,N resting proliferation rate constant [/d] +pA = AllPars(14) ;% S active proliferation rate constant [/d] +d = AllPars(15) ;% S,N resting death rate constant [/d] +dA = AllPars(16) ;% S,N acting death rate constant [/d] + +Xi = AllPars(17) ;% X initial condition (at dose) [pmol/kg] +ka = AllPars(18) ;% N803 absorption rate constant [/d] +ke = AllPars(19) ;% N803 elimination rate constant [/d] +vd = AllPars(20) ;% N803 'vol. of distribution'/'bioavailability' [L/kg] +C50 = AllPars(21) ;% 50% N803 saturation (drug effects) [pM] +pm = AllPars(22) ;% S,N induced resting expansion rate [] +alS = AllPars(23) ;% S activation stimulation factor [] +alN = AllPars(24) ;% N activation stimulation factor [] + +dR = AllPars(25) ;% R decay rate constant [/d] +sig = AllPars(26) ;% ratio of regulation generation rates sN/sS +ph = AllPars(27) ;% S,N proliferation regulation factor [] + +%% ------------------------------------------------------------------------ +% Calculate some initial conditions & parameters -------------------------- + +YS = Vi ./ (Vi + V50S) ;% initial V/(V50S+V) [cohort 1,2] +YN = Vi ./ (Vi + V50N) ;% initial V/(V50N+V) [cohort 1,2] + +K = YS + sig*YN ;% collection [cohort 1,2] +Ri = [ 1 K(2)/K(1) ] ;% initial regulation [cohort 1,2] +sS = dR/K(1) ;% R generation rate constant (S activation) [/d] +sN = sig*sS ;% R generation rate constant (N activation) [/d] + +ro = pm/p ;% S,N proliferation stimulation factor [] + +% calculate aS,aN,zeS,zeN +PS = p*S50./(S50+Si)./(1+ph*Ri) ;% initial S0 prolif rates [/d] ([c1,c2]) +PN = p*N50./(N50+Ni)./(1+ph*Ri) ;% initial N0 prolif rates [/d] ([c1,c2]) +US = 2*(2*pA/(pA+dA))^5/(mS+dA) ;% initial S8 / (AS*S0) +UN = 2 /(mN+dA) ;% initial N1 / (AN*N0) +AS = (d-PS) / (mS*US-1) ;% initial S0 activation [/d] ([c1,c2]) +AN = (d-PN) / (mN*UN-1) ;% initial N0 activation [/d] ([c1,c2]) +WS = AS(1)/AS(2)*YS(2)/YS(1) ; +WN = AN(1)/AN(2)*YN(2)/YN(1) ; +zeS = (WS-1) / (Ri(2)-WS) ;% S activation regulation factor [] +zeN = (WN-1) / (Ri(2)-WN) ;% N activation regulation factor [] +aS = AS(1)/YS(1) * (1+zeS) ;% S activation rate constant [/d] +aN = AN(1)/YN(1) * (1+zeN) ;% N activation rate constant [/d] + +% calculate initial S0,N0,SA for each cohort +ZS = US ; +for i = 1:5 + ZS = ZS + 2*(2*pA)^(i-1)/(pA+dA)^i ;% initial SA / (AS*S0) +end +S0 = Si./(1+ZS*AS) ;% initial S0 +N0 = Ni./(1+UN*AN) ;% initial N0 +SA = Si - S0 ;% initial SA (S1+S2...+S8) +NA = Ni - N0 ;% initial NA (N1) + +% calculate g and la +la = (SA(2) - SA(1)) / (SA(1)*Ri(2) - SA(2)) ;% S killing reg factor [] +g = q * (1+la) / SA(1) ;% S killing rate constant [uL/#-d] + +%% Do for each cohort (NOT indenting loop) ================================ +for c = 1:2 + +% solve for initial S1-6 +S = zeros(1,6) ;% initial S1-S8 +S(1) = 2*AS(c)*S0(c) / (dA+pA) ;% S1 +for i = 2:5 + S(i) = 2*pA*S(i-1) / (dA+pA) ;% S2 to S5 +end +S(6) = 2*pA*S(5) / (dA+mS) ;% S6 + +%% ------------------------------------------------------------------------ +% Prepare parameter and initial value vectors and call 'N803_model_2' ----- + +% V S0-8 N0 N1 X C R initial values +Yic = [ Vi(c) S0(c) S N0(c) NA(c) 0 0 Ri(c) Ri(c) ] ; + +Pars(01) = q ;% V growth rate constant [/d] +Pars(02) = g ;% S killing rate constant [uL/#-d] +Pars(03) = V50S ;% 50% V saturation for S activation [#/mL] +Pars(04) = V50N ;% 50% V saturation for N activation [#/mL] +Pars(05) = aS ;% S activation rate constant [/d] +Pars(06) = aN ;% N activation rate constant [/d] +Pars(07) = mS ;% S reversion rate constant [/d] +Pars(08) = mN ;% N reversion rate constant [/d] +Pars(09) = S50 ;% 50% S saturation for proliferation [#/uL] +Pars(10) = N50 ;% 50% N saturation for proliferation [#/uL] +Pars(11) = p ;% S,N resting proliferation rate constant [/d] +Pars(12) = pA ;% S active proliferation rate constant [/d] +Pars(13) = d ;% S,N resting death rate constant [/d] +Pars(14) = dA ;% S,N active death rate constant [/d] +Pars(15) = Xi ;% X initial condition (at dose) [pmol/kg] +Pars(16) = ka ;% N803 absorption rate constant [/d] +Pars(17) = ke ;% N803 elimination rate constant [/d] +Pars(18) = vd ;% N803 'volume of distribution'/'bioavailability' [L/kg] +Pars(19) = C50 ;% 50% N803 saturation (drug effects) [pM] +Pars(20) = ro ;% S,N proliferation stimulation factor [] +Pars(21) = alS ;% S activation stimulation factor [] +Pars(22) = alN ;% N activation stimulation factor [] +Pars(23) = sS ;% R generation rate constant (S activation) [/d] +Pars(24) = sN ;% R generation rate constant (N activation) [/d] +Pars(25) = dR ;% R decay rate constant [/d] +Pars(26) = la ;% S killing regulation factor [] +Pars(27) = ph ;% S,N proliferation regulation factor [] +Pars(28) = zeS ;% S activation regulation factor [] +Pars(29) = zeN ;% N activation regulation factor [] + +if any( [ Pars Yic ] < 0 ) && noNegs + error('Negative parameters or initial values.') +end + +if iscell(SoluTimes) % if calibrating to Ellis-Connell data + + VirusTimes = SoluTimes{2*(c-1)+1} ;% times for virus data points + CD8TcTimes = SoluTimes{2*(c-1)+2} ;% times for CD8 T cell data points + + X = sort( unique( [ VirusTimes ; CD8TcTimes ; DoseTimes{c} ] ) ) ; + + Y = N803_model_2A( X, DoseTimes{c}, Pars, Yic, varargin ) ; + + Y_OUT{2*(c-1)+1} = interp1( X, Y(:,1), VirusTimes ) ; + Y_OUT{2*(c-1)+2} = interp1( X, Y(:,2), CD8TcTimes ) ; + + PARS = [ PARS ; [ Pars Yic ] ] ; %#ok + +elseif cohort == c % if evaluating for one cohort + + Y_OUT = N803_model_2A( SoluTimes, DoseTimes, Pars, Yic, varargin ) ; + PARS = [ Pars Yic ] ; +end +end +end \ No newline at end of file diff --git a/N803_shared_2B.m b/N803_shared_2B.m new file mode 100644 index 0000000..963ec27 --- /dev/null +++ b/N803_shared_2B.m @@ -0,0 +1,224 @@ +%% N803_shared_2B.m - solves model 2B for 2 cohorts with shared parameters +% +% /--------------------------------------------------------------\ +% | Date: 05/06/2023 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Model 2B is a version of model 2 that models all CD8+ T cells as the +% SIV-specific cells from Model 2, which now only divide 4 times after +% activation. Also, regulation is generated with independent parameters. +% +% Nomenclature: V = SIV virions [#/uL] +% T8 = total CD8+ T cells [#/uL] +% E0 = resting CD8+ T cells [#/uL] +% EA = active CD8+ T cells [#/uL] +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = immune regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% NOTE: The format of 'SoluTimes' input will determine the format of the +% 'Y_OUT' output and determine whether the model will be run twice, from +% two different sets of initials. +% +% SoluTimes = ascending vector of days at which to evaluate solution +% OR +% SoluTimes = 1x4 cell vector, where each element is a vector of days +% corresponding to a data type (see OUTPUTS) +% +% AllPars = vector of parameters (see list in function) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% NOTE: 'varargin' is used to add optional inputs +% EX: N803_shared_2B(SoluTimes,AllPars,'Input Name',[Input Value]) +% +% DoseTimes = ascending vector of days at which to administer doses +% DEFAULT: 1x2 cell vector, where each element is a vector of days +% corresponding to doses in Ellis-Connell data +% +% cohort = scalar to run model for just one cohort ('1' or '2') +% NOTE: 'cohort' is ignored if 'SoluTimes' is a cell vector +% +% noNegs = logical scalar +% if true, throw an error if parameters are negative +% +% Additional inputs will also be passed to 'N803_model_2B.m' and used to +% define options in the same manner (see function for list) +% EX: N803_shared_2B(SoluTimes,DoseTimes,AllPars,'AbsTol',1e-2} +% will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% NOTE: The format of 'SoluTimes' input will determine the format of the +% 'Y_OUT' output and determine whether the model will be run twice, from +% two different sets of initials. +% +% If 'SoluTimes' is a 1x4 cell vector +% Y_OUT{1} = V at points in 'SoluTimes{1}' [log(#/mL)] cohort 1 +% Y_OUT{2} = T8 at points in 'SoluTimes{2}' [#/uL] cohort 1 +% Y_OUT{3} = V at points in 'SoluTimes{3}' [log(#/mL)] cohort 2 +% Y_OUT{4} = T8 at points in 'SoluTimes{4}' [#/uL] cohort 2 +% PARS(1,:) = [ parameters, initials ] for cohort 1 (see code) +% PARS(2,:) = [ parameters, initials ] for cohort 2 (see code) +% +% If 'SoluTimes' is a vector of days +% Y_OUT = model outputs as described in 'N803_model_2B.m' +% PARS = [ parameters, initials ] (see code) +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function [Y_OUT,PARS] = N803_shared_2B(SoluTimes,AllPars,varargin) + +% Declare optional input default values (then overwrite, if specified) +DoseTimes = { 7*[0:3 5:8 37:40]' ... dose days for cohort 1 + , 7*(0:2:4)' };% dose days for cohort 2 +cohort = 1 ;% run for this cohort (ignored if 'SoluTimes' is a cell) +noNegs = true ;% do not allow negative parameter values + +for n = 1:length(varargin) +if strcmp('DoseTimes', varargin(n)) ; DoseTimes = varargin{n+1} ; end +if strcmp('cohort', varargin(n)) ; cohort = varargin{n+1} ; end +if strcmp('noNegs', varargin(n)) ; noNegs = varargin{n+1} ; end +end + +Y_OUT = SoluTimes ;% initialize response output +PARS = [] ;% initialize parameter output + +% Rename inputed parameters ----------------------------------------------- +Vi(1) = AllPars(01) ;% V initial [#/mL] (cohort 1) +Vi(2) = AllPars(02) ;% V initial [#/mL] (cohort 2) +Ei(1) = AllPars(03) ;% E initial [#/uL] (cohort 1) +Ei(2) = AllPars(04) ;% E initial [#/uL] (cohort 2) + +q = AllPars(05) ;% V growth rate constant [/d] +V50E = AllPars(06) ;% 50% V saturation for E activation [#/mL] +V50R = AllPars(07) ;% 50% V saturation for R generation [#/mL] +m = AllPars(08) ;% E reversion rate constant [/d] + +E50 = AllPars(09) ;% 50% E saturation for proliferation [#/uL] +p = AllPars(10) ;% E resting proliferation rate constant [/d] +pA = AllPars(11) ;% E active proliferation rate constant [/d] +d = AllPars(12) ;% E resting death rate constant [/d] +dA = AllPars(13) ;% E active death rate constant [/d] + +Xi = AllPars(14) ;% X initial condition (at dose) [pmol/kg] +ka = AllPars(15) ;% N803 absorption rate constant [/d] +ke = AllPars(16) ;% N803 elimination rate constant [/d] +vd = AllPars(17) ;% N803 'vol. of distribution'/'bioavailability' [L/kg] +C50 = AllPars(18) ;% 50% N803 saturation (drug effects) [pM] +pm = AllPars(19) ;% E0 induced resting expansion rate [] +al = AllPars(20) ;% E activation stimulation factor [] +si = AllPars(21) ;% R generation stimulation factor [] + +dR = AllPars(22) ;% R decay rate constant [/d] +ph = AllPars(23) ;% E proliferation regulation factor [] + +%% ------------------------------------------------------------------------ +% Calculate some initial conditions & parameters -------------------------- + +YE = Vi ./ (Vi+ V50E) ;% initial V/(V50E+V) [cohort 1,2] +YR = Vi ./ (Vi+ V50R) ;% initial V/(V50R+V) [cohort 1,2] + +Ri = [ 1 , YR(2)/YR(1) ] ;% initial regulation [cohort 1,2] +s = dR / YR(1) ;% R generation rate constant [/d] + +ro = pm/p ;% E proliferation stimulation factor + +% calculate a,ze +Pi = p*E50./(E50+Ei)./(1+ph*Ri) ;% initial E0 prolif rates [/d] ([c1,c2]) +U = 2/(m+dA) * (2*pA/(pA+dA))^3 ;% initial E4 / (A*E0) +Ai = (d-Pi) / (m*U-1) ;% initial E activation [/d] ([c1,c2]) +W = Ai(1)/Ai(2)*YE(2)/YE(1) ; +ze = (W-1) / (Ri(2)-W) ;% E activation regulation factor [] +a = Ai(1)/YE(1) * (1+ze) ;% E activation rate constant [/d] + +% calculate initial E0,EA for each cohort +Z = U ; +for i = 1:3 + Z = Z + 2*(2*pA)^(i-1)/(pA+dA)^i ;% initial EA / (A*E0) +end +E0 = Ei./(1+Z*Ai) ;% initial E0 +EA = E0.*(Z*Ai) ;% initial EA (E1+E2+E3+E4) + +% calculate g and la +la = (EA(2) - EA(1)) / (EA(1)*Ri(2) - EA(2)) ;% E killing reg factor [] +g = q * (1+la) / EA(1) ;% E killing rate constant [uL/#-d] + +%% Do for each cohort (NOT indenting loop) ================================ +for c = 1:2 + +% solve for initial E1-4 +E = zeros(1,4) ;% initial E1-E4 +E(1) = 2*Ai(c)*E0(c) / (pA+dA) ;% E1 +for i = 2:3 + E(i) = 2*pA*E(i-1) / (pA+dA) ;% E2 to E3 +end +E(4) = 2*pA*E(3) / (m+dA) ;% E4 + +%% ------------------------------------------------------------------------ +% Prepare parameter and initial value vectors and call 'N803_model_2' ----- + +% V E0-4 X C R initial values +Yic = [ Vi(c) E0(c) E 0 0 Ri(c) Ri(c) ] ; + +Pars(01) = q ;% V growth rate constant [/d] +Pars(02) = g ;% E killing rate constant [uL/#-d] +Pars(03) = V50E ;% 50% V saturation for E activation [#/mL] +Pars(04) = V50R ;% 50% V saturation for R generation [#/mL] +Pars(05) = a ;% E activation rate constant [/d] +Pars(06) = m ;% E reversion rate constant [/d] +Pars(07) = E50 ;% 50% E saturation for proliferation [#/uL] +Pars(08) = p ;% E resting proliferation rate constant [/d] +Pars(09) = pA ;% E active proliferation rate constant [/d] +Pars(10) = d ;% E resting death rate constant [/d] +Pars(11) = dA ;% E active death rate constant [/d] +Pars(12) = Xi ;% X initial condition (at dose) [pmol/kg] +Pars(13) = ka ;% N803 absorption rate constant [/d] +Pars(14) = ke ;% N803 elimination rate constant [/d] +Pars(15) = vd ;% N803 'volume of distribution'/'bioavailability' [L/kg] +Pars(16) = C50 ;% 50% N803 saturation (drug effects) [pM] +Pars(17) = ro ;% E proliferation stimulation factor [] +Pars(18) = al ;% E activation stimulation factor [] +Pars(19) = si ;% R generation stimulation factor [] +Pars(20) = s ;% R generation rate constant [/d] +Pars(21) = dR ;% R decay rate constant [/d] +Pars(22) = la ;% E killing regulation factor [] +Pars(23) = ph ;% E proliferation regulation factor [] +Pars(24) = ze ;% E activation regulation factor [] + +if any( [ Pars Yic ] < 0 ) && noNegs + error('Negative parameters or initial values.') +end + +if iscell(SoluTimes) % if calibrating to Ellis-Connell data + + VirusTimes = SoluTimes{2*(c-1)+1} ;% times for virus data points + CD8TcTimes = SoluTimes{2*(c-1)+2} ;% times for CD8 T cell data points + + X = sort( unique( [ VirusTimes ; CD8TcTimes ; DoseTimes{c} ] ) ) ; + + Y = N803_model_2B( X, DoseTimes{c}, Pars, Yic, varargin ) ; + + Y_OUT{2*(c-1)+1} = interp1( X, Y(:,1), VirusTimes ) ; + Y_OUT{2*(c-1)+2} = interp1( X, Y(:,2), CD8TcTimes ) ; + + PARS = [ PARS ; [ Pars Yic ] ] ; %#ok + +elseif cohort == c % if evaluating for one cohort + + Y_OUT = N803_model_2B( SoluTimes, DoseTimes, Pars, Yic, varargin ) ; + PARS = [ Pars Yic ] ; + +end +end +end \ No newline at end of file diff --git a/N803_shared_2C.m b/N803_shared_2C.m new file mode 100644 index 0000000..a0ef6e1 --- /dev/null +++ b/N803_shared_2C.m @@ -0,0 +1,219 @@ +%% N803_shared_2C.m - solves model 2C for 2 cohorts with shared parameters +% +% /--------------------------------------------------------------\ +% | Date: 05/07/2023 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Model 2C is a version of model 2 where CD8+ T cells are represented by +% one resting and one active group. Proliferation of both groups is +% promoted by N803, limited by cell density, and inhibited by regulation. +% Also, regulation is generated with independent parameters. +% +% Nomenclature: V = SIV virions [#/uL] +% T8 = total CD8+ T cells [#/uL] +% E0 = resting CD8+ T cells [#/uL] +% EA = active CD8+ T cells [#/uL] +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = immune regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% NOTE: The format of 'SoluTimes' input will determine the format of the +% 'Y_OUT' output and determine whether the model will be run twice, from +% two different sets of initials. +% +% SoluTimes = ascending vector of days at which to evaluate solution +% OR +% SoluTimes = 1x4 cell vector, where each element is a vector of days +% corresponding to a data type (see OUTPUTS) +% +% AllPars = vector of parameters (see list in function) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% NOTE: 'varargin' is used to add optional inputs +% EX: N803_shared_2C(SoluTimes,AllPars,'Input Name',[Input Value]) +% +% DoseTimes = ascending vector of days at which to administer doses +% DEFAULT: 1x2 cell vector, where each element is a vector of days +% corresponding to doses in Ellis-Connell data +% +% cohort = scalar to run model for just one cohort ('1' or '2') +% NOTE: 'cohort' is ignored if 'SoluTimes' is a cell vector +% +% noNegs = logical scalar +% if true, throw an error if parameters are negative +% +% Additional inputs will also be passed to 'N803_model_2C.m' and used to +% define options in the same manner (see function for list) +% EX: N803_shared_2C(SoluTimes,DoseTimes,AllPars,'AbsTol',1e-2} +% will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% NOTE: The format of 'SoluTimes' input will determine the format of the +% 'Y_OUT' output and determine whether the model will be run twice, from +% two different sets of initials. +% +% If 'SoluTimes' is a 1x4 cell vector +% Y_OUT{1} = V at points in 'SoluTimes{1}' [log(#/mL)] cohort 1 +% Y_OUT{2} = T8 at points in 'SoluTimes{2}' [#/uL] cohort 1 +% Y_OUT{3} = V at points in 'SoluTimes{3}' [log(#/mL)] cohort 2 +% Y_OUT{4} = T8 at points in 'SoluTimes{4}' [#/uL] cohort 2 +% PARS(1,:) = [ parameters, initials ] for cohort 1 (see code) +% PARS(2,:) = [ parameters, initials ] for cohort 2 (see code) +% +% If 'SoluTimes' is a vector of days +% Y_OUT = model outputs as described in 'N803_model_2C.m' +% PARS = [ parameters, initials ] (see code) +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function [Y_OUT,PARS] = N803_shared_2C(SoluTimes,AllPars,varargin) + +% Declare optional input default values (then overwrite, if specified) +DoseTimes = { 7*[0:3 5:8 37:40]' ... dose days for cohort 1 + , 7*(0:2:4)' };% dose days for cohort 2 +cohort = 1 ;% run for this cohort (ignored if 'SoluTimes' is a cell) +noNegs = true ;% do not allow negative parameter values + +for n = 1:length(varargin) +if strcmp('DoseTimes', varargin(n)) ; DoseTimes = varargin{n+1} ; end +if strcmp('cohort', varargin(n)) ; cohort = varargin{n+1} ; end +if strcmp('noNegs', varargin(n)) ; noNegs = varargin{n+1} ; end +end + +Y_OUT = SoluTimes ;% initialize output vector +PARS = [] ;% initialize parameter matrix + +% Rename inputed parameters ----------------------------------------------- +Vi(1) = AllPars(01) ;% V initial [#/mL] (cohort 1) +Vi(2) = AllPars(02) ;% V initial [#/mL] (cohort 2) +Ei(1) = AllPars(03) ;% E initial [#/uL] (cohort 1) +Ei(2) = AllPars(04) ;% E initial [#/uL] (cohort 2) + +q = AllPars(05) ;% V growth rate constant [/d] +V50E = AllPars(06) ;% 50% V saturation for E activation [#/mL] +V50R = AllPars(07) ;% 50% V saturation for R generation [#/mL] +m = AllPars(08) ;% E reversion rate constant [/d] + +E50 = AllPars(09) ;% 50% E saturation for proliferation [#/uL] +p = AllPars(10) ;% E0 proliferation rate constant [/d] +pn = AllPars(11) ;% normalized EA proliferation rate constant [/d] +d = AllPars(12) ;% E0 death rate constant [/d] +dA = AllPars(13) ;% EA death rate constant [/d] + +Xi = AllPars(14) ;% X initial condition (at dose) [pmol/kg] +ka = AllPars(15) ;% N803 absorption rate constant [/d] +ke = AllPars(16) ;% N803 elimination rate constant [/d] +vd = AllPars(17) ;% N803 'vol. of distribution'/'bioavailability' [L/kg] +C50 = AllPars(18) ;% 50% N803 saturation (drug effects) [pM] +pm = AllPars(19) ;% E0 maximum proliferation rate [/d] +pAm = AllPars(20) ;% EA maximum proliferation rate [/d] +al = AllPars(21) ;% E activation stimulation factor [] +si = AllPars(22) ;% R generation stimulation factor [] + +dR = AllPars(23) ;% R decay rate constant [/d] +ph = AllPars(24) ;% E0 proliferation regulation factor [] +phA = AllPars(25) ;% EA proliferation regulation factor [] + +%% ------------------------------------------------------------------------ +% Calculate some initial conditions & parameters -------------------------- + +YE = Vi ./ (Vi+ V50E) ;% initial V/(V50E+V) [cohort 1,2] +YR = Vi ./ (Vi+ V50R) ;% initial V/(V50R+V) [cohort 1,2] + +Ri = [ 1 , YR(2)/YR(1) ] ;% initial regulation [cohort 1,2] +s = dR / YR(1) ;% R generation rate constant [/d] + +pA = p + pn*(pAm-p) ;% E active proliferation rate constant [/d] +ro = pm/p - 1 ;% E0 proliferation stimulation factor +roA = pAm/pA - 1 ;% EA proliferation stimulation factor +D = E50./(E50+Ei) ;% density-dependence multiplier + +P0i = p *D./(1+ph *Ri) ;% initial E0 prolif rates [/d] ([c1,c2]) +PAi = pA*D./(1+phA*Ri) ;% initial EA prolif rates [/d] ([c1,c2]) +U = 2./(m+dA-PAi) ;% initial EA/(A*E0) ([c1,c2]) +Ai = (d-P0i) ./ (m*U-1) ;% initial E activation [/d] ([c1,c2]) + +W = Ai(1)/Ai(2)*YE(2)/YE(1) ; +ze = (W-1) / (Ri(2)-W) ;% E activation regulation factor [] +a = Ai(1)/YE(1) * (1+ze) ;% E activation rate constant [/d] + +E0 = Ei./(1+U.*Ai) ;% initial E0 +EA = E0.*(U.*Ai) ;% initial EA + +la = (EA(1)-EA(2)) / (EA(2)-EA(1)*Ri(2)) ;% E killing regulation factor [] +g = q * (1+la) / EA(1) ;% E killing rate constant [uL/#-d] + +%% Do for each cohort (NOT indenting loop) ================================ +for c = 1:2 + +%% ------------------------------------------------------------------------ +% Prepare parameter and initial value vectors and call 'N803_model_2' ----- + +% V E0 EA X C R initial values +Yic = [ Vi(c) E0(c) EA(c) 0 0 Ri(c) Ri(c) ] ; + +Pars(01) = q ;% V growth rate constant [/d] +Pars(02) = g ;% E killing rate constant [uL/#-d] +Pars(03) = V50E ;% 50% V saturation for E activation [#/mL] +Pars(04) = V50R ;% 50% V saturation for R generation [#/mL] +Pars(05) = a ;% E activation rate constant [/d] +Pars(06) = m ;% E reversion rate constant [/d] +Pars(07) = E50 ;% 50% E saturation for proliferation [#/uL] +Pars(08) = p ;% E0 proliferation rate constant [/d] +Pars(09) = pA ;% EA proliferation rate constant [/d] +Pars(10) = d ;% E0 death rate constant [/d] +Pars(11) = dA ;% EA death rate constant [/d] +Pars(12) = Xi ;% X initial condition (at dose) [pmol/kg] +Pars(13) = ka ;% N803 absorption rate constant [/d] +Pars(14) = ke ;% N803 elimination rate constant [/d] +Pars(15) = vd ;% N803 'volume of distribution'/'bioavailability' [L/kg] +Pars(16) = C50 ;% 50% N803 saturation (drug effects) [pM] +Pars(17) = ro ;% E0 proliferation stimulation factor [] +Pars(18) = roA ;% EA proliferation stimulation factor [] +Pars(19) = al ;% E activation stimulation factor [] +Pars(20) = si ;% R generation stimulation factor [] +Pars(21) = s ;% R generation rate constant [/d] +Pars(22) = dR ;% R decay rate constant [/d] +Pars(23) = la ;% E killing regulation factor [] +Pars(24) = ph ;% E0 proliferation regulation factor [] +Pars(25) = phA ;% EA proliferation regulation factor [] +Pars(26) = ze ;% E activation regulation factor [] + +if any( [ Pars Yic ] < 0 ) && noNegs + error('Negative parameters or initial values.') +end + +if iscell(SoluTimes) % if calibrating to Ellis-Connell data + + VirusTimes = SoluTimes{2*(c-1)+1} ;% times for virus data points + CD8TcTimes = SoluTimes{2*(c-1)+2} ;% times for CD8 T cell data points + + X = sort( unique( [ VirusTimes ; CD8TcTimes ; DoseTimes{c} ] ) ) ; + + Y = N803_model_2C( X, DoseTimes{c}, Pars, Yic, varargin ) ; + + Y_OUT{2*(c-1)+1} = interp1( X, Y(:,1), VirusTimes ) ; + Y_OUT{2*(c-1)+2} = interp1( X, Y(:,2), CD8TcTimes ) ; + + PARS = [ PARS ; [ Pars Yic ] ] ; %#ok + +elseif cohort == c % if evaluating for one cohort + + Y_OUT = N803_model_2C( SoluTimes, DoseTimes, Pars, Yic, varargin ) ; + PARS = [ Pars Yic ] ; + +end +end +end \ No newline at end of file diff --git a/N803_shared_2D.m b/N803_shared_2D.m new file mode 100644 index 0000000..1b4ec47 --- /dev/null +++ b/N803_shared_2D.m @@ -0,0 +1,250 @@ +%% N803_shared_2D.m - solves model 2D for 2 cohorts with shared parameters +% +% /--------------------------------------------------------------\ +% | Date: 05/06/2023 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Model 2D is a version of model 2 that includes infected cell killing by +% non-SIV-specific CD8+ T cells. +% +% Nomenclature: V = SIV virions [#/mL] +% T8 = total CD8+ T cells [#/uL] +% S0 = resting SIV-specific CD8+ T cells [#/uL] +% SA = active SIV-specific CD8+ T cells [#/uL] (S1-S8) +% N0 = resting non-SIV-specific CD8+ T cells [#/uL] +% NA = active non-SIV-specific CD8+ T cells [#/uL] (N1) +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = immune regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% NOTE: The format of 'SoluTimes' input will determine the format of the +% 'Y_OUT' output and determine whether the model will be run twice, from +% two different sets of initials. +% +% SoluTimes = ascending vector of days at which to evaluate solution +% OR +% SoluTimes = 1x4 cell vector, where each element is a vector of days +% corresponding to a data type (see OUTPUTS) +% +% AllPars = vector of parameters (see list in function) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% NOTE: 'varargin' is used to add optional inputs +% EX: N803_shared_2(SoluTimes,AllPars,'Input Name',[Input Value]) +% +% DoseTimes = ascending vector of days at which to administer doses +% DEFAULT: 1x2 cell vector, where each element is a vector of days +% corresponding to doses in Ellis-Connell data +% +% cohort = scalar to run model for just one cohort ('1' or '2') +% NOTE: 'cohort' is ignored if 'SoluTimes' is a cell vector +% +% noNegs = logical scalar +% if true, throw an error if parameters are negative +% +% Additional inputs will also be passed to 'N803_model_2.m' and used to +% define options in the same manner (see function for list) +% EX: N803_shared_2(SoluTimes,DoseTimes,AllPars,'AbsTol',1e-2} +% will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% NOTE: The format of 'SoluTimes' input will determine the format of the +% 'Y_OUT' output and determine whether the model will be run twice, from +% two different sets of initials. +% +% If 'SoluTimes' is a 1x4 cell vector +% Y_OUT{1} = V at points in 'SoluTimes{1}' [log(#/mL)] cohort 1 +% Y_OUT{2} = T8 at points in 'SoluTimes{2}' [#/uL] cohort 1 +% Y_OUT{3} = V at points in 'SoluTimes{3}' [log(#/mL)] cohort 2 +% Y_OUT{4} = T8 at points in 'SoluTimes{4}' [#/uL] cohort 2 +% PARS(1,:) = [ parameters, initials ] for cohort 1 (see code) +% PARS(2,:) = [ parameters, initials ] for cohort 2 (see code) +% +% If 'SoluTimes' is a vector of days +% Y_OUT = model outputs as described in 'N803_model_2.m' +% PARS = [ parameters, initials ] (see code) +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function [Y_OUT,PARS] = N803_shared_2D(SoluTimes,AllPars,varargin) + +% Declare optional input default values (then overwrite, if specified) +DoseTimes = { 7*[0:3 5:8 37:40]' ... dose days for cohort 1 + , 7*(0:2:4)' };% dose days for cohort 2 +cohort = 1 ;% run for this cohort (ignored if 'SoluTimes' is a cell) +noNegs = true ;% do not allow negative parameter values + +for n = 1:length(varargin) +if strcmp('DoseTimes', varargin(n)) ; DoseTimes = varargin{n+1} ; end +if strcmp('cohort', varargin(n)) ; cohort = varargin{n+1} ; end +if strcmp('noNegs', varargin(n)) ; noNegs = varargin{n+1} ; end +end + +Y_OUT = SoluTimes ;% initialize response output +PARS = [] ;% initialize parameter output + +% Rename inputed parameters ----------------------------------------------- +Vi(1) = AllPars(01) ;% V initial [#/mL] (cohort 1) +Vi(2) = AllPars(02) ;% V initial [#/mL] (cohort 2) +SNi(1) = AllPars(03) ;% S+N initial [#/uL] (cohort 1) +SNi(2) = AllPars(04) ;% S+N initial [#/uL] (cohort 2) +fS(1) = AllPars(05) ;% initial frequency: S/(S+N) (cohort 1) +fS(2) = AllPars(06) ;% initial frequency: S/(S+N) (cohort 2) +Si = SNi.*fS ;% S initial value [#/uL] [cohort 1,2] +Ni = SNi-Si ;% N initial value [#/uL] [cohort 1,2] + +q = AllPars(07) ;% V growth rate constant [/d] +psi = AllPars(08) ;% gN/gS killing rate constant ratio +V50S = AllPars(09) ;% 50% V saturation for S activation [#/mL] +V50N = AllPars(10) ;% 50% V saturation for N activation [#/mL] +mS = AllPars(11) ;% S reversion rate constant [/d] +mN = AllPars(12) ;% N reversion rate constant [/d] + +S50 = AllPars(13)* fS(1) ;% 50% S saturation for proliferation [#/uL] +N50 = AllPars(13)*(1-fS(1)) ;% 50% N saturation for proliferation [#/uL] +p = AllPars(14) ;% S,N resting proliferation rate constant [/d] +pA = AllPars(15) ;% S active proliferation rate constant [/d] +d = AllPars(16) ;% S,N resting death rate constant [/d] +dA = AllPars(17) ;% S,N acting death rate constant [/d] + +Xi = AllPars(18) ;% X initial condition (at dose) [pmol/kg] +ka = AllPars(19) ;% N803 absorption rate constant [/d] +ke = AllPars(20) ;% N803 elimination rate constant [/d] +vd = AllPars(21) ;% N803 'vol. of distribution'/'bioavailability' [L/kg] +C50 = AllPars(22) ;% 50% N803 saturation (drug effects) [pM] +pm = AllPars(23) ;% S,N induced resting expansion rate [] +alS = AllPars(24) ;% S activation stimulation factor [] +alN = AllPars(25) ;% N activation stimulation factor [] + +dR = AllPars(26) ;% R decay rate constant [/d] +sig = AllPars(27) ;% ratio of regulation generation rates sN/sS +ph = AllPars(28) ;% S,N proliferation regulation factor [] + +%% ------------------------------------------------------------------------ +% Calculate some initial conditions & parameters -------------------------- + +YS = Vi ./ (Vi + V50S) ;% initial V/(V50S+V) [cohort 1,2] +YN = Vi ./ (Vi + V50N) ;% initial V/(V50N+V) [cohort 1,2] + +K = YS + sig*YN ;% collection [cohort 1,2] +Ri = [ 1 K(2)/K(1) ] ;% initial regulation [cohort 1,2] +sS = dR/K(1) ;% R generation rate constant (S activation) [/d] +sN = sig*sS ;% R generation rate constant (N activation) [/d] + +ro = pm/p ;% S,N proliferation stimulation factor [] + +% calculate aS,aN,zeS,zeN +PS = p*S50./(S50+Si)./(1+ph*Ri) ;% initial S0 prolif rates [/d] ([c1,c2]) +PN = p*N50./(N50+Ni)./(1+ph*Ri) ;% initial N0 prolif rates [/d] ([c1,c2]) +US = 2*(2*pA/(pA+dA))^7/(mS+dA) ;% initial S8 / (AS*S0) +UN = 2 /(mN+dA) ;% initial N1 / (AN*N0) +AS = (d-PS) / (mS*US-1) ;% initial S0 activation [/d] ([c1,c2]) +AN = (d-PN) / (mN*UN-1) ;% initial N0 activation [/d] ([c1,c2]) +WS = AS(1)/AS(2)*YS(2)/YS(1) ; +WN = AN(1)/AN(2)*YN(2)/YN(1) ; +zeS = (WS-1) / (Ri(2)-WS) ;% S activation regulation factor [] +zeN = (WN-1) / (Ri(2)-WN) ;% N activation regulation factor [] +aS = AS(1)/YS(1) * (1+zeS) ;% S activation rate constant [/d] +aN = AN(1)/YN(1) * (1+zeN) ;% N activation rate constant [/d] + +% calculate initial S0,N0,SA for each cohort +ZS = US ; +for i = 1:7 + ZS = ZS + 2*(2*pA)^(i-1)/(pA+dA)^i ;% initial SA / (AS*S0) +end +S0 = Si./(1+ZS*AS) ;% initial S0 +N0 = Ni./(1+UN*AN) ;% initial N0 +SA = Si - S0 ;% initial SA (S1+S2...+S8) +NA = Ni - N0 ;% initial NA (N1) + +% calculate g and la +K = SA + psi.*NA ; +la = (K(2) - K(1)) / (K(1)*Ri(2) - K(2)) ;% S,N killing reg factor [] +gS = q * (1+la) / K(1) ;% S killing rate constant [uL/#-d] +gN = psi*gS ;% N killing rate constant [uL/#-d] + +%% Do for each cohort (NOT indenting loop) ================================ +for c = 1:2 + +% solve for initial S1-8 +S = zeros(1,8) ;% initial S1-S8 +S(1) = 2*AS(c)*S0(c) / (dA+pA) ;% S1 +for i = 2:7 + S(i) = 2*pA*S(i-1) / (dA+pA) ;% S2 to S7 +end +S(8) = 2*pA*S(7) / (dA+mS) ;% S8 + +%% ------------------------------------------------------------------------ +% Prepare parameter and initial value vectors and call 'N803_model_2' ----- + +% V S0-8 N0 N1 X C R initial values +Yic = [ Vi(c) S0(c) S N0(c) NA(c) 0 0 Ri(c) Ri(c) ] ; + +Pars(01) = q ;% V growth rate constant [/d] +Pars(02) = gS ;% S killing rate constant [uL/#-d] +Pars(03) = gN ;% N killing rate constant [uL/#-d] +Pars(04) = V50S ;% 50% V saturation for S activation [#/mL] +Pars(05) = V50N ;% 50% V saturation for N activation [#/mL] +Pars(06) = aS ;% S activation rate constant [/d] +Pars(07) = aN ;% N activation rate constant [/d] +Pars(08) = mS ;% S reversion rate constant [/d] +Pars(09) = mN ;% N reversion rate constant [/d] +Pars(10) = S50 ;% 50% S saturation for proliferation [#/uL] +Pars(11) = N50 ;% 50% N saturation for proliferation [#/uL] +Pars(12) = p ;% S,N resting proliferation rate constant [/d] +Pars(13) = pA ;% S active proliferation rate constant [/d] +Pars(14) = d ;% S,N resting death rate constant [/d] +Pars(15) = dA ;% S,N active death rate constant [/d] +Pars(16) = Xi ;% X initial condition (at dose) [pmol/kg] +Pars(17) = ka ;% N803 absorption rate constant [/d] +Pars(18) = ke ;% N803 elimination rate constant [/d] +Pars(19) = vd ;% N803 'volume of distribution'/'bioavailability' [L/kg] +Pars(20) = C50 ;% 50% N803 saturation (drug effects) [pM] +Pars(21) = ro ;% S,N proliferation stimulation factor [] +Pars(22) = alS ;% S activation stimulation factor [] +Pars(23) = alN ;% N activation stimulation factor [] +Pars(24) = sS ;% R generation rate constant (S activation) [/d] +Pars(25) = sN ;% R generation rate constant (N activation) [/d] +Pars(26) = dR ;% R decay rate constant [/d] +Pars(27) = la ;% S,N killing regulation factor [] +Pars(28) = ph ;% S,N proliferation regulation factor [] +Pars(29) = zeS ;% S activation regulation factor [] +Pars(30) = zeN ;% N activation regulation factor [] + +if any( [ Pars Yic ] < 0 ) && noNegs + error('Negative parameters or initial values.') +end + +if iscell(SoluTimes) % if calibrating to Ellis-Connell data + + VirusTimes = SoluTimes{2*(c-1)+1} ;% times for virus data points + CD8TcTimes = SoluTimes{2*(c-1)+2} ;% times for CD8 T cell data points + + X = sort( unique( [ VirusTimes ; CD8TcTimes ; DoseTimes{c} ] ) ) ; + + Y = N803_model_2D( X, DoseTimes{c}, Pars, Yic, varargin ) ; + + Y_OUT{2*(c-1)+1} = interp1( X, Y(:,1), VirusTimes ) ; + Y_OUT{2*(c-1)+2} = interp1( X, Y(:,2), CD8TcTimes ) ; + + PARS = [ PARS ; [ Pars Yic ] ] ; %#ok + +elseif cohort == c % if evaluating for one cohort + + Y_OUT = N803_model_2D( SoluTimes, DoseTimes, Pars, Yic, varargin ) ; + PARS = [ Pars Yic ] ; +end +end +end \ No newline at end of file diff --git a/N803_single_2.m b/N803_single_2.m new file mode 100644 index 0000000..522f25a --- /dev/null +++ b/N803_single_2.m @@ -0,0 +1,220 @@ +%% N803_single_2.m - solves model 2 for 1 cohort +% +% /--------------------------------------------------------------\ +% | Date: 05/06/2023 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Nomenclature: V = SIV virions [#/mL] +% T8 = total CD8+ T cells [#/uL] +% S0 = resting SIV-specific CD8+ T cells [#/uL] +% SA = active SIV-specific CD8+ T cells [#/uL] (S1-S8) +% N0 = resting non-SIV-specific CD8+ T cells [#/uL] +% NA = active non-SIV-specific CD8+ T cells [#/uL] (N1) +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = immune regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% NOTE: The format of 'SoluTimes' input will determine the format of the +% 'Y_OUT' output. +% +% SoluTimes = ascending vector of days at which to evaluate solution +% OR +% SoluTimes = 1x2 cell vector, where each element is a vector of days +% corresponding to a data type (see OUTPUTS) +% +% AllPars = vector of parameters (see list in function) +% +% cohort = scalar to indicate cohort #1 or #2 (determines dose) +% +%% ======================================================================== +% OPTIONS +% ======================================================================== +% NOTE: 'varargin' is used to add optional inputs +% EX: N803_single_2(SoluTimes,AllPars,'Input Name',[Input Value]) +% +% DoseTimes = ascending vector of days at which to administer doses +% DEFAULT: 1x2 cell vector, where each element is a vector of days +% corresponding to doses in Ellis-Connell data +% +% noNegs = logical scalar +% if true, throw an error if parameters are negative +% +% Additional inputs will also be passed to 'N803_model_2.m' and used to +% define options in the same manner (see function for list) +% EX: N803_single_2(SoluTimes,DoseTimes,AllPars,'AbsTol',1e-2} +% will set ode solver absolute tolerance to 1e-2 +% +%% ======================================================================== +% OUTPUTS +% ======================================================================== +% NOTE: The format of 'SoluTimes' input will determine the format of the +% 'Y_OUT' output. +% +% If 'SoluTimes' is a 1x2 cell vector +% Y_OUT{1} = V at points in 'SoluTimes{1}' [log(#/mL)] +% Y_OUT{2} = T8 at points in 'SoluTimes{2}' [#/uL] +% PARS = [ parameters, initials ] (see code) +% +% If 'SoluTimes' is a vector of days +% Y_OUT = model outputs as described in 'N803_model_2.m' +% PARS = [ parameters, initials ] (see code) +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function [Y_OUT,PARS] = N803_single_2(SoluTimes,AllPars,cohort,varargin) + +% Declare optional input default values (then overwrite, if specified) +DoseTimes = { 7*[0:3 5:8 37:40]' ... dose days for cohort 1 + , 7*(0:2:4)' };% dose days for cohort 2 +noNegs = true ;% do not allow negative parameter values + +for n = 1:length(varargin) +if strcmp('DoseTimes', varargin(n)) ; DoseTimes = varargin{n+1} ; end +if strcmp('noNegs', varargin(n)) ; noNegs = varargin{n+1} ; end +end + +Y_OUT = SoluTimes ;% initialize response output + +% Rename inputed parameters ----------------------------------------------- +Vi = AllPars(01) ;% V initial [#/mL] +SNi = AllPars(02) ;% S+N initial [#/uL] +fS = AllPars(03) ;% initial frequency: S/(S+N) +Si = SNi*fS ;% S initial value [#/uL] +Ni = SNi-Si ;% N initial value [#/uL] + +q = AllPars(04) ;% V growth rate constant [/d] +V50S = AllPars(05) ;% 50% V saturation for S activation [#/mL] +V50N = AllPars(06) ;% 50% V saturation for N activation [#/mL] +mS = AllPars(07) ;% S reversion rate constant [/d] +mN = AllPars(08) ;% N reversion rate constant [/d] + +S50 = AllPars(09)* fS(1) ;% 50% S saturation for proliferation [#/uL] +N50 = AllPars(09)*(1-fS(1)) ;% 50% N saturation for proliferation [#/uL] +p = AllPars(10) ;% S,N resting proliferation rate constant [/d] +pA = AllPars(11) ;% S active proliferation rate constant [/d] +d = AllPars(12) ;% S,N resting death rate constant [/d] +dA = AllPars(13) ;% S,N acting death rate constant [/d] + +Xi = AllPars(14) ;% X initial condition (at dose) [pmol/kg] +ka = AllPars(15) ;% N803 absorption rate constant [/d] +ke = AllPars(16) ;% N803 elimination rate constant [/d] +vd = AllPars(17) ;% N803 'vol. of distribution'/'bioavailability' [L/kg] +C50 = AllPars(18) ;% 50% N803 saturation (drug effects) [pM] +pm = AllPars(19) ;% S,N induced resting expansion rate [] +alS = AllPars(20) ;% S activation stimulation factor [] +alN = AllPars(21) ;% N activation stimulation factor [] + +dR = AllPars(22) ;% R decay rate constant [/d] +sig = AllPars(23) ;% ratio of regulation generation rates sN/sS +la = AllPars(24) ;% S killing regulation factor [] +ph = AllPars(25) ;% S,N proliferation regulation factor [] +zeS = AllPars(26) ;% S activation regulation factor [] +zeN = AllPars(27) ;% N activation regulation factor [] + +%% ------------------------------------------------------------------------ +% Calculate some initial conditions & parameters -------------------------- + +YS = Vi / (Vi + V50S) ;% initial V/(V50S+V) +YN = Vi / (Vi + V50N) ;% initial V/(V50N+V) + +sS = dR/(YS + sig*YN) ;% R generation rate constant (S activation) [/d] +sN = sig*sS ;% R generation rate constant (N activation) [/d] + +ro = pm/p ;% S,N proliferation stimulation factor [] + +% calculate aS,aN +PS = p*S50/(S50+Si)/(1+ph) ;% initial S0 prolif rates [/d] +PN = p*N50/(N50+Ni)/(1+ph) ;% initial N0 prolif rates [/d] +US = 2*(2*pA/(pA+dA))^7/(mS+dA) ;% initial S8 / (AS*S0) +UN = 2 /(mN+dA) ;% initial N1 / (AN*N0) +AS = (d-PS) / (mS*US-1) ;% initial S0 activation [/d] +AN = (d-PN) / (mN*UN-1) ;% initial N0 activation [/d] +aS = AS/YS * (1+zeS) ;% S activation rate constant [/d] +aN = AN/YN * (1+zeN) ;% N activation rate constant [/d] + +% calculate initial S0,N0,SA for each cohort +ZS = US ; +for i = 1:7 + ZS = ZS + 2*(2*pA)^(i-1)/(pA+dA)^i ;% initial SA / (AS*S0) +end +S0 = Si/(1+ZS*AS) ;% initial S0 +N0 = Ni/(1+UN*AN) ;% initial N0 +SA = Si - S0 ;% initial SA (S1+S2...+S8) +NA = Ni - N0 ;% initial NA (N1) + +% calculate g +g = q * (1+la) / SA ;% S killing rate constant [uL/#-d] + +% solve for initial S1-8 +S = zeros(1,8) ;% initial S1-S8 +S(1) = 2*AS*S0 / (dA+pA) ;% S1 +for i = 2:7 + S(i) = 2*pA*S(i-1) / (dA+pA) ;% S2 to S7 +end +S(8) = 2*pA*S(7) / (dA+mS) ;% S8 + +%% ------------------------------------------------------------------------ +% Prepare parameter and initial value vectors and call 'N803_model_2' ----- + +% V S0-8 N0 N1 X C R initial values +Yic = [ Vi S0 S N0 NA 0 0 1 1 ] ; + +Pars(01) = q ;% V growth rate constant [/d] +Pars(02) = g ;% S killing rate constant [uL/#-d] +Pars(03) = V50S ;% 50% V saturation for S activation [#/mL] +Pars(04) = V50N ;% 50% V saturation for N activation [#/mL] +Pars(05) = aS ;% S activation rate constant [/d] +Pars(06) = aN ;% N activation rate constant [/d] +Pars(07) = mS ;% S reversion rate constant [/d] +Pars(08) = mN ;% N reversion rate constant [/d] +Pars(09) = S50 ;% 50% S saturation for proliferation [#/uL] +Pars(10) = N50 ;% 50% N saturation for proliferation [#/uL] +Pars(11) = p ;% S,N resting proliferation rate constant [/d] +Pars(12) = pA ;% S active proliferation rate constant [/d] +Pars(13) = d ;% S,N resting death rate constant [/d] +Pars(14) = dA ;% S,N active death rate constant [/d] +Pars(15) = Xi ;% X initial condition (at dose) [pmol/kg] +Pars(16) = ka ;% N803 absorption rate constant [/d] +Pars(17) = ke ;% N803 elimination rate constant [/d] +Pars(18) = vd ;% N803 'volume of distribution'/'bioavailability' [L/kg] +Pars(19) = C50 ;% 50% N803 saturation (drug effects) [pM] +Pars(20) = ro ;% S,N proliferation stimulation factor [] +Pars(21) = alS ;% S activation stimulation factor [] +Pars(22) = alN ;% N activation stimulation factor [] +Pars(23) = sS ;% R generation rate constant (S activation) [/d] +Pars(24) = sN ;% R generation rate constant (N activation) [/d] +Pars(25) = dR ;% R decay rate constant [/d] +Pars(26) = la ;% S killing regulation factor [] +Pars(27) = ph ;% S,N proliferation regulation factor [] +Pars(28) = zeS ;% S activation regulation factor [] +Pars(29) = zeN ;% N activation regulation factor [] + +if any( [ Pars Yic ] < 0 ) && noNegs + error('Negative parameters or initial values.') +end + +if iscell(SoluTimes) % if calibrating to Ellis-Connell data + + VirusTimes = SoluTimes{1} ;% times for virus data points + CD8TcTimes = SoluTimes{2} ;% times for CD8 T cell data points + + X = sort( unique( [ VirusTimes ; CD8TcTimes ; DoseTimes{cohort} ] ) ) ; + + Y = N803_model_2( X, DoseTimes{cohort}, Pars, Yic, varargin ) ; + + Y_OUT{1} = interp1( X, Y(:,1), VirusTimes ) ; + Y_OUT{2} = interp1( X, Y(:,2), CD8TcTimes ) ; + +else + Y_OUT = N803_model_2( SoluTimes, DoseTimes, Pars, Yic, varargin ) ; +end + PARS = [ Pars Yic ] ; +end \ No newline at end of file diff --git a/main_shared_2A_calibration.m b/main_shared_2A_calibration.m new file mode 100644 index 0000000..2572625 --- /dev/null +++ b/main_shared_2A_calibration.m @@ -0,0 +1,139 @@ +%% main.m = calls model analysis functions stored in 'Analyzer' folder + +addpath Analyzer +addpath Data + +%% ======================================================================== +% NOTES +% ======================================================================== + +% Script performs MSLS parameter estimation for an ODE model of N-803 +% treatment of SIV. + +% Training data is from two SIV cohorts: +% Cohort 1 (low VL) = Ellis et. al. J Virol. 2018. (N803_DataSet_1.mat) +% Cohort 2 (high VL) = Ellis et. al. J Virol. 2022. (N803_DataSet_2.mat) + +% Model 2A is a version of model 2 where SIV-specific cells divide 6 times +% after activation, instead of 8. + +% Model is run for both cohorts (from two fixed points), with parameters +% shared between cohorts, fitting to both cohorts simultaneously. +% Parameter and IC calculation is done in 'N803_shared_2A.m'. +% ODEs are solved in 'N803_model_2A.m'. +% Files within 'Analyzer' folder perform MSLS. + +%% ======================================================================== +% INPUTS for 'calibrator.m' +% ======================================================================== + +% GuessBOUNDS = array for boundaries of initial parameter guesses ----- +GuessBOUNDS(:,01) = [1 15]*1e3 ;% V initial [#/mL] (cohort 1) +GuessBOUNDS(:,02) = [1 40]*1e5 ;% V initial [#/mL] (cohort 2) +GuessBOUNDS(:,03) = [200 1000] ;% S+N initial [#/uL] (cohort 1) +GuessBOUNDS(:,04) = [200 1000] ;% S+N initial [#/uL] (cohort 2) +GuessBOUNDS(:,05) = [.03 .3] ;% initial frequency: S/(S+N) (cohort 1) +GuessBOUNDS(:,06) = [.03 .3] ;% initial frequency: S/(S+N) (cohort 2) + +GuessBOUNDS(:,07) = [.01 1] ;% V growth rate constant [/d] +GuessBOUNDS(:,08) = [1e4 1e6] ;% 50% V saturation for S activation [#/mL] +GuessBOUNDS(:,09) = [1e4 1e6] ;% 50% V saturation for N activation [#/mL] +GuessBOUNDS(:,10) = [.004 .04] ;% S reversion rate constant [/d] +GuessBOUNDS(:,11) = [.04 .4] ;% N reversion rate constant [/d] + +GuessBOUNDS(:,12) = [300 3e3] ;% 50% S+N saturation for prolif. [#/uL] +GuessBOUNDS(:,13) = [.01 .1] ;% S,N resting prolif rate constant [/d] +GuessBOUNDS(:,14) = [1 4] ;% S active proliferation rate constant [/d] +GuessBOUNDS(:,15) = [.01 .05] ;% S,N resting death rate constant [/d] +GuessBOUNDS(:,16) = [.1 .5] ;% S,N active death rate constant [/d] + +GuessBOUNDS(:,17) = [1 1]*880 ;% X initial condition (at dose) [pmol/kg] +GuessBOUNDS(:,18) = [1 1]*0.8 ;% N803 absorption rate constant [/d] +GuessBOUNDS(:,19) = [1 1]*2.1 ;% N803 elimination rate constant [/d] +GuessBOUNDS(:,20) = [1 1]*1.3 ;% N803 'vol of dist'/'bioavail' [L/kg] +GuessBOUNDS(:,21) = [1 1e3] ;% 50% N803 saturation (drug effects) [pM] +GuessBOUNDS(:,22) = [.1 2] ;% S,N induced resting expansion rate [] +GuessBOUNDS(:,23) = [0.1 1e3] ;% S activation stimulation factor [] +GuessBOUNDS(:,24) = [0.01 100] ;% N activation stimulation factor [] + +GuessBOUNDS(:,25) = [0.05 1] ;% R decay rate constant [/d] +GuessBOUNDS(:,26) = [0.1 100] ;% ratio of regulation gen rates sN/sS +GuessBOUNDS(:,27) = [0.01 100] ;% S,N proliferation regulation factor + +% FitBOUNDS = array for boundaries of parameter fitting ------------------- +FitBOUNDS = GuessBOUNDS ; + +% X and Y = cell arrays for data x-values and y-values -------------------- +load('N803_DataSet_1') ; DS1 = DataSet ; +load('N803_DataSet_2') ; DS2 = DataSet ; +X = { DS1(1).X_data , DS1(3).X_data , DS2(1).X_data , DS2(3).X_data } ; +Y = { DS1(1).Y_data , DS1(3).Y_data , DS2(1).Y_data , DS2(3).Y_data } ; +% Apply log scaling to virus (see 'prepper' function) +[X,Y,Sid] = prepper(X,Y, [], [], [], [1 0 1 0]) ; + +% modelFunction = handle for calling model function ----------------------- +modelFunction = @(P) N803_shared_2A( X, P,... + 'AbsTol', 1e-3, 'RelTol', 1e-2, 'odeSolver', '23s' ,'timeOut', 1 ) ; + +% costFunction = handle for calling cost function ------------------------- +costFunction = @(P) costFun_NLL( P, modelFunction, Y,... + 'ScalingFacts', [1 .001 1 .001] ) ;% coverting T8 to #/nL + +% costFunction(GuessBOUNDS(1,:)) ;% uncomment to debug your costFunction + +% filePrefix = string to prefix to file names ----------------------------- +filePrefix = 'N803_shared_2A' ; + +% timeLimit = number of minutes allowed for a single optimization run ----- +timeLimit = 30 ; + +% numberGuesses = number of initial parameter guesses --------------------- +numberGuesses = 20000 ; + +% numberWorkers = number of parallel processes (enter [] or 0 for serial) +numberWorkers = 4 ; + +%% ------------------------------------------------------------------------ +% Optional inputs for results summary figures ('calibrator_plot.m') +% ------------------------------------------------------------------------ + +plotSetup.Xdata = X ; +plotSetup.Ydata = Y ; +plotSetup.Fits2Plot = 1:20 ; +plotSetup.nFits2Sum = 20 ; +plotSetup.SubjectID = Sid ; +plotSetup.XLabels = {{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'}} ; +plotSetup.YLabels = {{'Virus','[log CEQ/mL]'} ... + ,{'CD8^{+} T cells','[#/uL]'}... + ,{'Virus','[log CEQ/mL]'} ... + ,{'CD8^{+} T cells','[#/uL]'}} ; +plotSetup.XRANGE = 7*[ -4 41 ; -4 41 ; -4 6 ; -4 6 ] ; +plotSetup.YRANGE = [ 1 8 ; 0 3500 ; 1 8 ; 0 3500 ] ; +plotSetup.PlotGrid = [3 2] ; +plotSetup.PlotOuts = [1 3 2 4] ; +plotSetup.plotPars = 1 ; + +% Set up function handle for plotting +Xeval = X ; +xStep = 1 ; +for n = 1:length(Xeval) + Xeval{n} = ( plotSetup.XRANGE(n,1) : xStep : plotSetup.XRANGE(n,2) )' ; +end +plotSetup.modelFunction = @(P) N803_shared_2A(Xeval,P) ; +plotSetup.Xeval = Xeval ; + +%% ======================================================================== +% PERFORM MSLS +% ======================================================================== + +% Perform calibration using 'calibrator.m' +setup = struct('GuessBOUNDS',GuessBOUNDS,'FitBOUNDS',FitBOUNDS,... + 'costFunction',costFunction,'timeLimit',timeLimit) ; + +Results = calibrator(setup,filePrefix,numberGuesses,numberWorkers) ; + +% Plot results using 'calibrator_plot.m' +calibrator_plot(Results,plotSetup,filePrefix) diff --git a/main_shared_2B_calibration.m b/main_shared_2B_calibration.m new file mode 100644 index 0000000..0fc517f --- /dev/null +++ b/main_shared_2B_calibration.m @@ -0,0 +1,136 @@ +%% main.m = calls model analysis functions stored in 'Analyzer' folder + +addpath Analyzer +addpath Data + +%% ======================================================================== +% NOTES +% ======================================================================== + +% Script performs MSLS parameter estimation for an ODE model of N-803 +% treatment of SIV. + +% Training data is from two SIV cohorts: +% Cohort 1 (low VL) = Ellis et. al. J Virol. 2018. (N803_DataSet_1.mat) +% Cohort 2 (high VL) = Ellis et. al. J Virol. 2022. (N803_DataSet_2.mat) + +% Model 2B is a version of model 2 that models all CD8+ T cells as the +% SIV-specific cells from Model 2, which now only divide 4 times after +% activation. Also, regulation is generated with independent parameters. + +% Model is run for both cohorts (from two fixed points), with parameters +% shared between cohorts, fitting to both cohorts simultaneously. +% Parameter and IC calculation is done in 'N803_shared_2B.m'. +% ODEs are solved in 'N803_model_2B.m'. +% Files within 'Analyzer' folder perform MSLS. + +%% ======================================================================== +% INPUTS for 'calibrator.m' +% ======================================================================== + +% GuessBOUNDS = array for boundaries of initial parameter guesses ----- +GuessBOUNDS(:,01) = [1 15]*1e3 ;% V initial [#/mL] (cohort 1) +GuessBOUNDS(:,02) = [1 40]*1e5 ;% V initial [#/mL] (cohort 2) +GuessBOUNDS(:,03) = [200 1000] ;% E initial [#/uL] (cohort 1) +GuessBOUNDS(:,04) = [200 1000] ;% E initial [#/uL] (cohort 2) + +GuessBOUNDS(:,05) = [.01 1] ;% V growth rate constant [/d] +GuessBOUNDS(:,06) = [1e4 1e6] ;% 50% V saturation for E activation [#/mL] +GuessBOUNDS(:,07) = [1e4 1e6] ;% 50% V saturation for R generation [#/mL] +GuessBOUNDS(:,08) = [.004 .4] ;% E reversion rate constant [/d] + +GuessBOUNDS(:,09) = [300 3e3] ;% 50% E saturation for proliferation [#/uL] +GuessBOUNDS(:,10) = [.01 .1] ;% E resting prolif rate constant [/d] +GuessBOUNDS(:,11) = [.5 4] ;% E active proliferation rate constant [/d] +GuessBOUNDS(:,12) = [.01 .05] ;% E resting death rate constant [/d] +GuessBOUNDS(:,13) = [.1 .5] ;% E active death rate constant [/d] + +GuessBOUNDS(:,14) = [1 1]*880 ;% X initial condition (at dose) [pmol/kg] +GuessBOUNDS(:,15) = [1 1]*0.8 ;% N803 absorption rate constant [/d] +GuessBOUNDS(:,16) = [1 1]*2.1 ;% N803 elimination rate constant [/d] +GuessBOUNDS(:,17) = [1 1]*1.3 ;% N803 'vol of dist'/'bioavail' [L/kg] +GuessBOUNDS(:,18) = [1 1e3] ;% 50% N803 stimulation concentration [pM] +GuessBOUNDS(:,19) = [.1 2] ;% E0 induced resting expansion rate [] +GuessBOUNDS(:,20) = [0.1 1e3] ;% E activation stimulation factor [] +GuessBOUNDS(:,21) = [0.01 100] ;% R generation stimulation factor [] + +GuessBOUNDS(:,22) = [0.05 1] ;% R decay rate constant [/d] +GuessBOUNDS(:,23) = [0.01 100] ;% E proliferation regulation factor + +% FitBOUNDS = array for boundaries of parameter fitting ------------------- +FitBOUNDS = GuessBOUNDS ; + +% X and Y = cell arrays for data x-values and y-values -------------------- +load('N803_DataSet_1') ; DS1 = DataSet ; +load('N803_DataSet_2') ; DS2 = DataSet ; +X = { DS1(1).X_data , DS1(3).X_data , DS2(1).X_data , DS2(3).X_data } ; +Y = { DS1(1).Y_data , DS1(3).Y_data , DS2(1).Y_data , DS2(3).Y_data } ; +% Apply log scaling to virus (see 'prepper' function) +[X,Y,Sid] = prepper(X,Y, [], [], [], [1 0 1 0]) ; + +% modelFunction = handle for calling model function ----------------------- +modelFunction = @(P) N803_shared_2B( X, P,... + 'AbsTol', 1e-3, 'RelTol', 1e-2, 'odeSolver', '23s' ,'timeOut', 1 ) ; + +% costFunction = handle for calling cost function ------------------------- +costFunction = @(P) costFun_NLL( P, modelFunction, Y,... + 'ScalingFacts', [1 .001 1 .001] ) ;% coverting T8 to #/nL + +% costFunction(GuessBOUNDS(1,:)) ;% uncomment to debug your costFunction + +% filePrefix = string to prefix to file names ----------------------------- +filePrefix = 'N803_shared_2B' ; + +% timeLimit = number of minutes allowed for a single optimization run ----- +timeLimit = 30 ; + +% numberGuesses = number of initial parameter guesses --------------------- +numberGuesses = 10000 ; + +% numberWorkers = number of parallel processes (enter [] or 0 for serial) +numberWorkers = 4 ; + +%% ------------------------------------------------------------------------ +% Optional inputs for results summary figures ('calibrator_plot.m') +% ------------------------------------------------------------------------ + +plotSetup.Xdata = X ; +plotSetup.Ydata = Y ; +plotSetup.Fits2Plot = 1:20 ; +plotSetup.nFits2Sum = 20 ; +plotSetup.SubjectID = Sid ; +plotSetup.XLabels = {{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'}} ; +plotSetup.YLabels = {{'Virus','[log CEQ/mL]'} ... + ,{'CD8^{+} T cells','[#/uL]'}... + ,{'Virus','[log CEQ/mL]'} ... + ,{'CD8^{+} T cells','[#/uL]'}} ; +plotSetup.XRANGE = 7*[ -4 41 ; -4 41 ; -4 6 ; -4 6 ] ; +plotSetup.YRANGE = [ 1 8 ; 0 3500 ; 1 8 ; 0 3500 ] ; +plotSetup.PlotGrid = [3 2] ; +plotSetup.PlotOuts = [1 3 2 4] ; +plotSetup.plotPars = 1 ; + +% Set up function handle for plotting +Xeval = X ; +xStep = 1 ; +for n = 1:length(Xeval) + Xeval{n} = ( plotSetup.XRANGE(n,1) : xStep : plotSetup.XRANGE(n,2) )' ; +end +plotSetup.modelFunction = @(P) N803_shared_2B(Xeval,P) ; +plotSetup.Xeval = Xeval ; + +%% ======================================================================== +% PERFORM MSLS +% ======================================================================== + +% Perform calibration using 'calibrator.m' +setup = struct('GuessBOUNDS',GuessBOUNDS,'FitBOUNDS',FitBOUNDS,... + 'costFunction',costFunction,'timeLimit',timeLimit) ; + +Results = calibrator(setup,filePrefix,numberGuesses,numberWorkers) ; + +% Plot results using 'calibrator_plot.m' +calibrator_plot(Results,plotSetup,filePrefix) diff --git a/main_shared_2C_calibration.m b/main_shared_2C_calibration.m new file mode 100644 index 0000000..6cf3f74 --- /dev/null +++ b/main_shared_2C_calibration.m @@ -0,0 +1,139 @@ +%% main.m = calls model analysis functions stored in 'Analyzer' folder + +addpath Analyzer +addpath Data + +%% ======================================================================== +% NOTES +% ======================================================================== + +% Script performs MSLS parameter estimation for an ODE model of N-803 +% treatment of SIV. + +% Training data is from two SIV cohorts: +% Cohort 1 (low VL) = Ellis et. al. J Virol. 2018. (N803_DataSet_1.mat) +% Cohort 2 (high VL) = Ellis et. al. J Virol. 2022. (N803_DataSet_2.mat) + +% Model 2C is a version of model 2 where CD8+ T cells are represented by +% one resting and one active group. Proliferation of both groups is +% promoted by N803, limited by cell density, and inhibited by regulation. +% Also, regulation is generated with independent parameters. + +% Model is run for both cohorts (from two fixed points), with parameters +% shared between cohorts, fitting to both cohorts simultaneously. +% Parameter and IC calculation is done in 'N803_shared_2C.m'. +% ODEs are solved in 'N803_model_2C.m'. +% Files within 'Analyzer' folder perform MSLS. + +%% ======================================================================== +% INPUTS for 'calibrator.m' +% ======================================================================== + +% GuessBOUNDS = array for boundaries of initial parameter guesses ----- +GuessBOUNDS(:,01) = [1 15]*1e3 ;% V initial [#/mL] (cohort 1) +GuessBOUNDS(:,02) = [1 40]*1e5 ;% V initial [#/mL] (cohort 2) +GuessBOUNDS(:,03) = [200 1000] ;% E initial [#/uL] (cohort 1) +GuessBOUNDS(:,04) = [200 1000] ;% E initial [#/uL] (cohort 2) + +GuessBOUNDS(:,05) = [.01 1] ;% V growth rate constant [/d] +GuessBOUNDS(:,06) = [1e4 1e6] ;% 50% V saturation for E activation [#/mL] +GuessBOUNDS(:,07) = [1e4 1e6] ;% 50% V saturation for R generation [#/mL] +GuessBOUNDS(:,08) = [.05 .5] ;% E reversion rate constant [/d] + +GuessBOUNDS(:,09) = [300 3e3] ;% 50% E saturation for proliferation [#/uL] +GuessBOUNDS(:,10) = [.01 .1] ;% E0 prolif rate constant [/d] +GuessBOUNDS(:,11) = [.01 1] ;% normalized EA prolif rate constant [/d] +GuessBOUNDS(:,12) = [.01 .05] ;% E0 death rate constant [/d] +GuessBOUNDS(:,13) = [.1 .5] ;% EA death rate constant [/d] + +GuessBOUNDS(:,14) = [1 1]*880 ;% X initial condition (at dose) [pmol/kg] +GuessBOUNDS(:,15) = [1 1]*0.8 ;% N803 absorption rate constant [/d] +GuessBOUNDS(:,16) = [1 1]*2.1 ;% N803 elimination rate constant [/d] +GuessBOUNDS(:,17) = [1 1]*1.3 ;% N803 'vol of dist'/'bioavail' [L/kg] +GuessBOUNDS(:,18) = [1 1e3] ;% 50% N803 stimulation concentration [pM] +GuessBOUNDS(:,19) = [.02 2] ;% E0 maximum proliferation rate [/d] +GuessBOUNDS(:,20) = [0.4 4] ;% EA maximum proliferation rate [/d] +GuessBOUNDS(:,21) = [.1 1e3] ;% E activation stimulation factor [] +GuessBOUNDS(:,22) = [.1 100] ;% R generation stimulation factor [] + +GuessBOUNDS(:,23) = [0.05 1] ;% R decay rate constant [/d] +GuessBOUNDS(:,24) = [.01 100] ;% E0 proliferation regulation factor +GuessBOUNDS(:,25) = [.01 100] ;% EA proliferation regulation factor + +% FitBOUNDS = array for boundaries of parameter fitting ------------------- +FitBOUNDS = GuessBOUNDS ; + +% X and Y = cell arrays for data x-values and y-values -------------------- +load('N803_DataSet_1') ; DS1 = DataSet ; +load('N803_DataSet_2') ; DS2 = DataSet ; +X = { DS1(1).X_data , DS1(3).X_data , DS2(1).X_data , DS2(3).X_data } ; +Y = { DS1(1).Y_data , DS1(3).Y_data , DS2(1).Y_data , DS2(3).Y_data } ; +% Apply log scaling to virus (see 'prepper' function) +[X,Y,Sid] = prepper(X,Y, [], [], [], [1 0 1 0]) ; + +% modelFunction = handle for calling model function ----------------------- +modelFunction = @(P) N803_shared_2C( X, P,... + 'AbsTol', 1e-3, 'RelTol', 1e-2, 'odeSolver', '23s' ,'timeOut', 1 ) ; + +% costFunction = handle for calling cost function ------------------------- +costFunction = @(P) costFun_NLL( P, modelFunction, Y,... + 'ScalingFacts', [1 .001 1 .001] ) ;% coverting T8 to #/nL + +% costFunction(GuessBOUNDS(1,:)) ;% uncomment to debug your costFunction + +% filePrefix = string to prefix to file names ----------------------------- +filePrefix = 'N803_shared_2C' ; + +% timeLimit = number of minutes allowed for a single optimization run ----- +timeLimit = 30 ; + +% numberGuesses = number of initial parameter guesses --------------------- +numberGuesses = 15000 ; + +% numberWorkers = number of parallel processes (enter [] or 0 for serial) +numberWorkers = 4 ; + +%% ------------------------------------------------------------------------ +% Optional inputs for results summary figures ('calibrator_plot.m') +% ------------------------------------------------------------------------ + +plotSetup.Xdata = X ; +plotSetup.Ydata = Y ; +plotSetup.Fits2Plot = 1:20 ; +plotSetup.nFits2Sum = 20 ; +plotSetup.SubjectID = Sid ; +plotSetup.XLabels = {{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'}} ; +plotSetup.YLabels = {{'Virus','[log CEQ/mL]'} ... + ,{'CD8^{+} T cells','[#/uL]'}... + ,{'Virus','[log CEQ/mL]'} ... + ,{'CD8^{+} T cells','[#/uL]'}} ; +plotSetup.XRANGE = 7*[ -4 41 ; -4 41 ; -4 6 ; -4 6 ] ; +plotSetup.YRANGE = [ 1 8 ; 0 3500 ; 1 8 ; 0 3500 ] ; +plotSetup.PlotGrid = [3 2] ; +plotSetup.PlotOuts = [1 3 2 4] ; +plotSetup.plotPars = 1 ; + +% Set up function handle for plotting +Xeval = X ; +xStep = 1 ; +for n = 1:length(Xeval) + Xeval{n} = ( plotSetup.XRANGE(n,1) : xStep : plotSetup.XRANGE(n,2) )' ; +end +plotSetup.modelFunction = @(P) N803_shared_2C(Xeval,P) ; +plotSetup.Xeval = Xeval ; + +%% ======================================================================== +% PERFORM MSLS +% ======================================================================== + +% Perform calibration using 'calibrator.m' +setup = struct('GuessBOUNDS',GuessBOUNDS,'FitBOUNDS',FitBOUNDS,... + 'costFunction',costFunction,'timeLimit',timeLimit) ; + +Results = calibrator(setup,filePrefix,numberGuesses,numberWorkers) ; + +% Plot results using 'calibrator_plot.m' +calibrator_plot(Results,plotSetup,filePrefix) diff --git a/main_shared_2D_calibration.m b/main_shared_2D_calibration.m new file mode 100644 index 0000000..00921f5 --- /dev/null +++ b/main_shared_2D_calibration.m @@ -0,0 +1,140 @@ +%% main.m = calls model analysis functions stored in 'Analyzer' folder + +addpath Analyzer +addpath Data + +%% ======================================================================== +% NOTES +% ======================================================================== + +% Script performs MSLS parameter estimation for an ODE model of N-803 +% treatment of SIV. + +% Model 2D is a version of model 2 that includes infected cell killing by +% non-SIV-specific CD8+ T cells. + +% Training data is from two SIV cohorts: +% Cohort 1 (low VL) = Ellis et. al. J Virol. 2018. (N803_DataSet_1.mat) +% Cohort 2 (high VL) = Ellis et. al. J Virol. 2022. (N803_DataSet_2.mat) + +% Model is run for both cohorts (from two fixed points), with parameters +% shared between cohorts, fitting to both cohorts simultaneously. +% Parameter and IC calculation is done in 'N803_shared_2D.m'. +% ODEs are solved in 'N803_model_2D.m'. +% Files within 'Analyzer' folder perform MSLS. + +%% ======================================================================== +% INPUTS for 'calibrator.m' +% ======================================================================== + +% GuessBOUNDS = array for boundaries of initial parameter guesses ----- +GuessBOUNDS(:,01) = [1 15]*1e3 ;% V initial [#/mL] (cohort 1) +GuessBOUNDS(:,02) = [1 40]*1e5 ;% V initial [#/mL] (cohort 2) +GuessBOUNDS(:,03) = [200 1000] ;% S+N initial [#/uL] (cohort 1) +GuessBOUNDS(:,04) = [200 1000] ;% S+N initial [#/uL] (cohort 2) +GuessBOUNDS(:,05) = [.03 .3] ;% initial frequency: S/(S+N) (cohort 1) +GuessBOUNDS(:,06) = [.03 .3] ;% initial frequency: S/(S+N) (cohort 2) + +GuessBOUNDS(:,07) = [.01 1] ;% V growth rate constant [/d] +GuessBOUNDS(:,08) = [.01 1] ;% gN/gS killing rate constant ratio +GuessBOUNDS(:,09) = [1e4 1e6] ;% 50% V saturation for S activation [#/mL] +GuessBOUNDS(:,10) = [1e4 1e6] ;% 50% V saturation for N activation [#/mL] +GuessBOUNDS(:,11) = [.004 .04] ;% S reversion rate constant [/d] +GuessBOUNDS(:,12) = [.04 .4] ;% N reversion rate constant [/d] + +GuessBOUNDS(:,13) = [300 3e3] ;% 50% S+N saturation for prolif. [#/uL] +GuessBOUNDS(:,14) = [.01 .1] ;% S,N resting prolif rate constant [/d] +GuessBOUNDS(:,15) = [1 4] ;% S active proliferation rate constant [/d] +GuessBOUNDS(:,16) = [.01 .05] ;% S,N resting death rate constant [/d] +GuessBOUNDS(:,17) = [.1 .5] ;% S,N active death rate constant [/d] + +GuessBOUNDS(:,18) = [1 1]*880 ;% X initial condition (at dose) [pmol/kg] +GuessBOUNDS(:,19) = [1 1]*0.8 ;% N803 absorption rate constant [/d] +GuessBOUNDS(:,20) = [1 1]*2.1 ;% N803 elimination rate constant [/d] +GuessBOUNDS(:,21) = [1 1]*1.3 ;% N803 'vol of dist'/'bioavail' [L/kg] +GuessBOUNDS(:,22) = [1 1e3] ;% 50% N803 saturation (drug effects) [pM] +GuessBOUNDS(:,23) = [.1 2] ;% S,N induced resting expansion rate [] +GuessBOUNDS(:,24) = [0.1 1e3] ;% S activation stimulation factor [] +GuessBOUNDS(:,25) = [0.01 100] ;% N activation stimulation factor [] + +GuessBOUNDS(:,26) = [0.05 1] ;% R decay rate constant [/d] +GuessBOUNDS(:,27) = [0.1 100] ;% ratio of regulation gen rates sN/sS +GuessBOUNDS(:,28) = [0.01 100] ;% S,N proliferation regulation factor + +% FitBOUNDS = array for boundaries of parameter fitting ------------------- +FitBOUNDS = GuessBOUNDS ; + +% X and Y = cell arrays for data x-values and y-values -------------------- +load('N803_DataSet_1') ; DS1 = DataSet ; +load('N803_DataSet_2') ; DS2 = DataSet ; +X = { DS1(1).X_data , DS1(3).X_data , DS2(1).X_data , DS2(3).X_data } ; +Y = { DS1(1).Y_data , DS1(3).Y_data , DS2(1).Y_data , DS2(3).Y_data } ; +% Apply log scaling to virus (see 'prepper' function) +[X,Y,Sid] = prepper(X,Y, [], [], [], [1 0 1 0]) ; + +% modelFunction = handle for calling model function ----------------------- +modelFunction = @(P) N803_shared_2D( X, P,... + 'AbsTol', 1e-3, 'RelTol', 1e-2, 'odeSolver', '23s' ,'timeOut', 1 ) ; + +% costFunction = handle for calling cost function ------------------------- +costFunction = @(P) costFun_NLL( P, modelFunction, Y,... + 'ScalingFacts', [1 .001 1 .001] ) ;% coverting T8 to #/nL + +% costFunction(GuessBOUNDS(1,:)) ;% uncomment to debug your costFunction + +% filePrefix = string to prefix to file names ----------------------------- +filePrefix = 'N803_shared_2D' ; + +% timeLimit = number of minutes allowed for a single optimization run ----- +timeLimit = 30 ; + +% numberGuesses = number of initial parameter guesses --------------------- +numberGuesses = 20000 ; + +% numberWorkers = number of parallel processes (enter [] or 0 for serial) +numberWorkers = 4 ; + +%% ------------------------------------------------------------------------ +% Optional inputs for results summary figures ('calibrator_plot.m') +% ------------------------------------------------------------------------ + +plotSetup.Xdata = X ; +plotSetup.Ydata = Y ; +plotSetup.Fits2Plot = 1:20 ; +plotSetup.nFits2Sum = 20 ; +plotSetup.SubjectID = Sid ; +plotSetup.XLabels = {{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'}} ; +plotSetup.YLabels = {{'Virus','[log CEQ/mL]'} ... + ,{'CD8^{+} T cells','[#/uL]'}... + ,{'Virus','[log CEQ/mL]'} ... + ,{'CD8^{+} T cells','[#/uL]'}} ; +plotSetup.XRANGE = 7*[ -4 41 ; -4 41 ; -4 6 ; -4 6 ] ; +plotSetup.YRANGE = [ 1 8 ; 0 3500 ; 1 8 ; 0 3500 ] ; +plotSetup.PlotGrid = [3 2] ; +plotSetup.PlotOuts = [1 3 2 4] ; +plotSetup.plotPars = 1 ; + +% Set up function handle for plotting +Xeval = X ; +xStep = 1 ; +for n = 1:length(Xeval) + Xeval{n} = ( plotSetup.XRANGE(n,1) : xStep : plotSetup.XRANGE(n,2) )' ; +end +plotSetup.modelFunction = @(P) N803_shared_2D(Xeval,P) ; +plotSetup.Xeval = Xeval ; + +%% ======================================================================== +% PERFORM MSLS +% ======================================================================== + +% Perform calibration using 'calibrator.m' +setup = struct('GuessBOUNDS',GuessBOUNDS,'FitBOUNDS',FitBOUNDS,... + 'costFunction',costFunction,'timeLimit',timeLimit) ; + +Results = calibrator(setup,filePrefix,numberGuesses,numberWorkers) ; + +% Plot results using 'calibrator_plot.m' +calibrator_plot(Results,plotSetup,filePrefix) diff --git a/main_shared_2_calibration.m b/main_shared_2_calibration.m new file mode 100644 index 0000000..d649597 --- /dev/null +++ b/main_shared_2_calibration.m @@ -0,0 +1,152 @@ +%% main.m = calls model analysis functions stored in 'Analyzer' folder + +addpath Analyzer +addpath Data + +%% ======================================================================== +% NOTES +% ======================================================================== + +% Script performs MSLS parameter estimation and Bayesian MCMC uncertainty +% quantification for an ODE model of N-803 treatment of SIV. + +% Training data is from two SIV cohorts: +% Cohort 1 (low VL) = Ellis et. al. J Virol. 2018. (N803_DataSet_1.mat) +% Cohort 2 (high VL) = Ellis et. al. J Virol. 2022. (N803_DataSet_2.mat) + +% Model is run for both cohorts (from two fixed points), with parameters +% shared between cohorts, fitting to both cohorts simultaneously. +% Parameter and IC calculation is done in 'N803_shared_2.m'. +% ODEs are solved in 'N803_model_2.m'. +% Files within 'Analyzer' folder perform MSLS and MCMC. + +%% ======================================================================== +% INPUTS for 'calibrator.m' +% ======================================================================== + +% GuessBOUNDS = array for boundaries of initial parameter guesses ----- +GuessBOUNDS(:,01) = [1 15]*1e3 ;% V initial [#/mL] (cohort 1) +GuessBOUNDS(:,02) = [1 40]*1e5 ;% V initial [#/mL] (cohort 2) +GuessBOUNDS(:,03) = [200 1000] ;% S+N initial [#/uL] (cohort 1) +GuessBOUNDS(:,04) = [200 1000] ;% S+N initial [#/uL] (cohort 2) +GuessBOUNDS(:,05) = [.03 .3] ;% initial frequency: S/(S+N) (cohort 1) +GuessBOUNDS(:,06) = [.03 .3] ;% initial frequency: S/(S+N) (cohort 2) + +GuessBOUNDS(:,07) = [.01 1] ;% V growth rate constant [/d] +GuessBOUNDS(:,08) = [1e4 1e6] ;% 50% V saturation for S activation [#/mL] +GuessBOUNDS(:,09) = [1e4 1e6] ;% 50% V saturation for N activation [#/mL] +GuessBOUNDS(:,10) = [.004 .04] ;% S reversion rate constant [/d] +GuessBOUNDS(:,11) = [.04 .4] ;% N reversion rate constant [/d] + +GuessBOUNDS(:,12) = [300 3e3] ;% 50% S+N saturation for prolif. [#/uL] +GuessBOUNDS(:,13) = [.01 .1] ;% S,N resting prolif rate constant [/d] +GuessBOUNDS(:,14) = [1 4] ;% S active proliferation rate constant [/d] +GuessBOUNDS(:,15) = [.01 .05] ;% S,N resting death rate constant [/d] +GuessBOUNDS(:,16) = [.1 .5] ;% S,N active death rate constant [/d] + +GuessBOUNDS(:,17) = [1 1]*880 ;% X initial condition (at dose) [pmol/kg] +GuessBOUNDS(:,18) = [1 1]*0.8 ;% N803 absorption rate constant [/d] +GuessBOUNDS(:,19) = [1 1]*2.1 ;% N803 elimination rate constant [/d] +GuessBOUNDS(:,20) = [1 1]*1.3 ;% N803 'vol of dist'/'bioavail' [L/kg] +GuessBOUNDS(:,21) = [1 1e3] ;% 50% N803 saturation (drug effects) [pM] +GuessBOUNDS(:,22) = [.1 2] ;% S,N induced resting expansion rate [] +GuessBOUNDS(:,23) = [0.1 1e3] ;% S activation stimulation factor [] +GuessBOUNDS(:,24) = [0.01 100] ;% N activation stimulation factor [] + +GuessBOUNDS(:,25) = [0.05 1] ;% R decay rate constant [/d] +GuessBOUNDS(:,26) = [0.1 100] ;% ratio of regulation gen rates sN/sS +GuessBOUNDS(:,27) = [0.01 100] ;% S,N proliferation regulation factor + +% FitBOUNDS = array for boundaries of parameter fitting ------------------- +FitBOUNDS = GuessBOUNDS ; + +% X and Y = cell arrays for data x-values and y-values -------------------- +load('N803_DataSet_1') ; DS1 = DataSet ; +load('N803_DataSet_2') ; DS2 = DataSet ; +X = { DS1(1).X_data , DS1(3).X_data , DS2(1).X_data , DS2(3).X_data } ; +Y = { DS1(1).Y_data , DS1(3).Y_data , DS2(1).Y_data , DS2(3).Y_data } ; +% Apply log scaling to virus (see 'prepper' function) +[X,Y,Sid] = prepper(X,Y, [], [], [], [1 0 1 0]) ; + +% modelFunction = handle for calling model function ----------------------- +modelFunction = @(P) N803_shared_2( X, P,... + 'AbsTol', 1e-3, 'RelTol', 1e-2, 'odeSolver', '23s' ,'timeOut', 1 ) ; + +% costFunction = handle for calling cost function ------------------------- +costFunction = @(P) costFun_NLL( P, modelFunction, Y,... + 'ScalingFacts', [1 .001 1 .001] ) ;% coverting T8 to #/nL + +% costFunction(GuessBOUNDS(1,:)) ;% uncomment to debug your costFunction + +% filePrefix = string to prefix to file names ----------------------------- +filePrefix = 'N803_shared_2' ; + +% timeLimit = number of minutes allowed for a single optimization run ----- +timeLimit = 30 ; + +% numberGuesses = number of initial parameter guesses --------------------- +numberGuesses = 20000 ; + +% numberWorkers = number of parallel processes (enter [] or 0 for serial) +numberWorkers = 0 ; + +%% ------------------------------------------------------------------------ +% Optional inputs for results summary figures ('calibrator_plot.m') +% ------------------------------------------------------------------------ + +plotSetup.Xdata = X ; +plotSetup.Ydata = Y ; +plotSetup.Fits2Plot = 1:20 ; +plotSetup.nFits2Sum = 20 ; +plotSetup.SubjectID = Sid ; +plotSetup.XLabels = {{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'}} ; +plotSetup.YLabels = {{'Virus','[log CEQ/mL]'} ... + ,{'CD8^{+} T cells','[#/uL]'}... + ,{'Virus','[log CEQ/mL]'} ... + ,{'CD8^{+} T cells','[#/uL]'}} ; +plotSetup.XRANGE = 7*[ -4 41 ; -4 41 ; -4 6 ; -4 6 ] ; +plotSetup.YRANGE = [ 1 8 ; 0 3500 ; 1 8 ; 0 3500 ] ; +plotSetup.PlotGrid = [3 2] ; +plotSetup.PlotOuts = [1 3 2 4] ; +plotSetup.plotPars = 1 ; + +% Set up function handle for plotting +Xeval = X ; +xStep = 1 ; +for n = 1:length(Xeval) + Xeval{n} = ( plotSetup.XRANGE(n,1) : xStep : plotSetup.XRANGE(n,2) )' ; +end +plotSetup.modelFunction = @(P) N803_shared_2(Xeval,P) ; +plotSetup.Xeval = Xeval ; + +%% ======================================================================== +% PERFORM MSLS +% ======================================================================== + +% Perform calibration using 'calibrator.m' +setup = struct('GuessBOUNDS',GuessBOUNDS,'FitBOUNDS',FitBOUNDS,... + 'costFunction',costFunction,'timeLimit',timeLimit) ; + +Results = calibrator(setup,filePrefix,numberGuesses,numberWorkers) ; + +% Plot results using 'calibrator_plot.m' +calibrator_plot(Results,plotSetup,filePrefix) + +%% ======================================================================== +% PERFORM MCMC +% ======================================================================== + +% InitialRanks = ranks (by NLL) of parameter sets to initialize MCMC ------ +InitialRanks = 1:5 ; + +% MCMC sample size -------------------------------------------------------- +numSamples = 2 + 1e6 ; + +% save every 'saveEach' iterations ---------------------------------------- +saveEach = 1000 ; + +% Perform uncertainty quantification using 'MCMCsampler.m' +% MCMCsampler([filePrefix '_Results'],InitialRanks,numSamples,saveEach) ; diff --git a/main_single_2_calibration.m b/main_single_2_calibration.m new file mode 100644 index 0000000..88be245 --- /dev/null +++ b/main_single_2_calibration.m @@ -0,0 +1,147 @@ +%% main.m = calls model analysis functions stored in 'Analyzer' folder + +addpath Analyzer +addpath Data + +%% ======================================================================== +% NOTES +% ======================================================================== + +% Script performs MSLS parameter estimation for an ODE model of N-803 +% treatment of SIV. + +% Training data is from two SIV cohorts: +% Cohort 1 (low VL) = Ellis et. al. J Virol. 2018. (N803_DataSet_1.mat) +% Cohort 2 (high VL) = Ellis et. al. J Virol. 2022. (N803_DataSet_2.mat) + +% Model is run for both cohorts (from two fixed points), with NO parameters +% shared between cohorts, fitting to both cohorts seperately. +% Parameter and IC calculation is done in 'N803_single.m'. +% ODEs are solved in 'N803_model_2.m'. +% Files within 'Analyzer' folder perform MSLS. + +%% ======================================================================== +% INPUTS for 'calibrator.m' +% ======================================================================== + +% GuessBOUNDS = array for boundaries of initial parameter guesses --------- +GuessBOUNDS(:,01) = [1 15]*1e3 ;% V initial [#/mL] (cohort 1) +GuessBOUNDS(:,02) = [200 1000] ;% S+N initial [#/uL] +GuessBOUNDS(:,03) = [.03 .3] ;% initial frequency: S/(S+N) (cohort 1) + +GuessBOUNDS(:,04) = [.01 1] ;% V growth rate constant [/d] +GuessBOUNDS(:,05) = [1e4 1e6] ;% 50% V saturation for S activation [#/mL] +GuessBOUNDS(:,06) = [1e4 1e6] ;% 50% V saturation for N activation [#/mL] +GuessBOUNDS(:,07) = [.004 .04] ;% S reversion rate constant [/d] +GuessBOUNDS(:,08) = [.04 .4] ;% N reversion rate constant [/d] + +GuessBOUNDS(:,09) = [300 3e3] ;% 50% S+N saturation for prolif. [#/uL] +GuessBOUNDS(:,10) = [.01 .1] ;% S,N resting prolif rate constant [/d] +GuessBOUNDS(:,11) = [1 4] ;% S active proliferation rate constant [/d] +GuessBOUNDS(:,12) = [.01 .05] ;% S,N resting death rate constant [/d] +GuessBOUNDS(:,13) = [.1 .5] ;% S,N active death rate constant [/d] + +GuessBOUNDS(:,14) = [1 1]*880 ;% X initial condition (at dose) [pmol/kg] +GuessBOUNDS(:,15) = [1 1]*0.8 ;% N803 absorption rate constant [/d] +GuessBOUNDS(:,16) = [1 1]*2.1 ;% N803 elimination rate constant [/d] +GuessBOUNDS(:,17) = [1 1]*1.3 ;% N803 'vol of dist'/'bioavail' [L/kg] +GuessBOUNDS(:,18) = [1 1e3] ;% 50% N803 saturation (drug effects) [pM] +GuessBOUNDS(:,19) = [.1 2] ;% S,N induced resting expansion rate [] +GuessBOUNDS(:,20) = [0.1 1e3] ;% S activation stimulation factor [] +GuessBOUNDS(:,21) = [0.01 100] ;% N activation stimulation factor [] + +GuessBOUNDS(:,22) = [0.05 1] ;% R decay rate constant [/d] +GuessBOUNDS(:,23) = [0.1 100] ;% ratio of regulation gen rates sN/sS +GuessBOUNDS(:,24) = [0.01 100] ;% S killing regulation factor +GuessBOUNDS(:,25) = [0.01 100] ;% S,N proliferation regulation factor +GuessBOUNDS(:,26) = [0.01 100] ;% S activation regulation factor +GuessBOUNDS(:,27) = [0.01 100] ;% N activation regulation factor + +% timeLimit = number of minutes allowed for a single optimization run ----- +timeLimit = 30 ; + +% numberGuesses = number of initial parameter guesses --------------------- +numberGuesses = 5000 ; + +% numberWorkers = number of parallel processes (enter [] or 0 for serial) +numberWorkers = 4 ; + +%% ------------------------------------------------------------------------ +% Do for each cohort... +for c = 1:2 + +% GuessBOUNDS = array for boundaries of initial parameter guesses --------- +if c == 2 + GuessBOUNDS(:,01) = [1 40]*1e5 ;% V initial [#/mL] (cohort 2) +end + +% FitBOUNDS = array for boundaries of parameter fitting ------------------- +FitBOUNDS = GuessBOUNDS ; + +% X and Y = cell arrays for data x-values and y-values -------------------- +load(['N803_DataSet_' num2str(c)]) +X = { DataSet(1).X_data , DataSet(3).X_data } ; +Y = { DataSet(1).Y_data , DataSet(3).Y_data } ; +% Apply log scaling to virus (see 'prepper' function) +[X,Y,Sid] = prepper(X,Y, [], [], [], [1 0]) ; + +% modelFunction = handle for calling model function ----------------------- +modelFunction = @(P) N803_single_2( X, P, c,... + 'AbsTol', 1e-3, 'RelTol', 1e-2, 'odeSolver', '23s' ,'timeOut', 1 ) ; + +% costFunction = handle for calling cost function ------------------------- +costFunction = @(P) costFun_NLL( P, modelFunction, Y,... + 'ScalingFacts', [1 .001 1 .001] ) ;% coverting T8 to #/nL + +% costFunction(GuessBOUNDS(1,:)) ;% uncomment to debug your costFunction + +% filePrefix = string to prefix to file names ----------------------------- +filePrefix = ['N803_cohort_' num2str(c)]; + +%% ------------------------------------------------------------------------ +% Optional inputs for results summary figures ('calibrator_plot.m') +% ------------------------------------------------------------------------ + +plotSetup.Xdata = X ; +plotSetup.Ydata = Y ; +plotSetup.Fits2Plot = 1:20 ; +plotSetup.nFits2Sum = 20 ; +plotSetup.SubjectID = Sid ; +plotSetup.XLabels = {{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'}} ; +plotSetup.YLabels = {{'Virus','[log CEQ/mL]'} ... + ,{'CD8^{+} T cells','[#/uL]'}} ; +switch c + case 1 + plotSetup.XRANGE = 7*[ -4 41 ; -4 41 ] ; + case 2 + plotSetup.XRANGE = 7*[ -4 6 ; -4 6] ; +end +plotSetup.YRANGE = [ 1 8 ; 0 3500 ] ; +plotSetup.PlotGrid = [2 1] ; +plotSetup.PlotOuts = 1:2 ; +plotSetup.plotPars = 0 ; + +% Set up function handle for plotting +Xeval = X ; +xStep = 1 ; +for n = 1:length(Xeval) + Xeval{n} = ( plotSetup.XRANGE(n,1) : xStep : plotSetup.XRANGE(n,2) )' ; +end +plotSetup.modelFunction = @(P) N803_single_2(Xeval,P,c) ; +plotSetup.Xeval = Xeval ; + +%% ======================================================================== +% FUNCTION CALLS +% ======================================================================== + +% Perform calibration using 'calibrator.m' +setup = struct('GuessBOUNDS',GuessBOUNDS,'FitBOUNDS',FitBOUNDS,... + 'costFunction',costFunction,'timeLimit',timeLimit) ; + +Results = calibrator(setup,filePrefix,numberGuesses,numberWorkers) ; + +% Plot results using 'calibrator_plot.m' +calibrator_plot(Results,plotSetup,filePrefix) + +end \ No newline at end of file