From 70d61bf678bde2b33e9397ffe6d21ea9eb24f037 Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Tue, 27 Dec 2022 11:00:25 -0600 Subject: [PATCH] Add files via upload --- N803_model_2S.m | 353 +++++++++++++++++++++++ N803_shared_noN.m | 254 +++++++--------- main_shared_calibration_2022_08_18.m | 158 ++++++++++ main_shared_noN_calibration_2022_11_02.m | 155 ++++++++++ main_singles_calibration_2022_09_25.m | 163 +++++++++++ 5 files changed, 936 insertions(+), 147 deletions(-) create mode 100644 N803_model_2S.m create mode 100644 main_shared_calibration_2022_08_18.m create mode 100644 main_shared_noN_calibration_2022_11_02.m create mode 100644 main_singles_calibration_2022_09_25.m diff --git a/N803_model_2S.m b/N803_model_2S.m new file mode 100644 index 0000000..4c0e754 --- /dev/null +++ b/N803_model_2S.m @@ -0,0 +1,353 @@ +%% N803_model_2S.m - solves ODE model for N-803 treatment of SIV +% +% /--------------------------------------------------------------\ +% | Date: 12/27/2022 | +% | Author: Jonathan Cody | +% | Affiliation: Purdue University | +% | Weldon School of Biomedical Engineering | +% | Pienaar Computational Systems Pharmacology Lab | +% \--------------------------------------------------------------/ +% +% Model 2S is a simplified version of Model 2 that combines SIV-specific +% and non-SIV-specific CD8+ T cells +% +% Nomenclature: V = SIV virions [#/無] +% T8 = total CD8+ T cells [#/無] +% E0 = resting CD8+ T cells [#/無] +% EA = active CD8+ T cells [#/無] +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = regulation [] (dimensionless quantity) +% +%% ======================================================================== +% INPUTS +% ======================================================================== +% +% 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') +% +% Pars = vector of parameters (see list on line 93) +% +% Yic = vector of initial conditions (see list on line 230) +% +%% ======================================================================== +% 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) = V at points in 'SoluTimes' [log fold change] +% Y_OUT(:,2) = T8 at points in 'SoluTimes' [fold change] +% +%% ======================================================================== +% FUNCTION +% ======================================================================== +function Y_OUT = N803_model_2S(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 = [] ;% if non-empty, more outputs are added (see line 170) +timeOut = 60 ;% scalar time limit [seconds] for ODE solver + % if exceeded, an error is thrown +warnOff = [] ;% if non-empty, 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('EqualPars', Opts(n)) ; EqualPars = Opts{n+1} ; end +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 inputs 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 + +% 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 ~isempty(warnOff) ; warning off ; end + +%% ------------------------------------------------------------------------ +% Rename inputed parameters ----------------------------------------------- + +q = Pars(01) ;% V growth rate (if E+B were absent) [/d] +g0 = Pars(02) ;% EA killing rate constant [無/#-d] +V50E = Pars(03) ;% 50% viral stimulation saturation for E [#/無] +V50R = Pars(04) ;% 50% viral stimulation saturation for R [#/無] +a0 = Pars(05) ;% E0 activation rate constant [/d] +m = Pars(06) ;% EA reversion rate constant [/d] + +E50 = Pars(07) ;% 50% E0 proliferation saturation [#/無] +p0 = 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 [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 stimulation concentration [pM] +p1 = Pars(17) ;% E0 proliferation stimulation factor [] +a1 = Pars(18) ;% E0 activation stimulation factor [] +s1 = Pars(19) ;% R generation stimulation factor [] + +s0 = Pars(20) ;% R generation maximum base rate [/d] +dR = Pars(21) ;% R decay rate constant [/d] +g2 = Pars(22) ;% EA killing regulation factor [] +p2 = Pars(23) ;% E0 proliferation regulation factor [] +a2 = Pars(24) ;% E0 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(11) = Yic(11) + Xi ;% add drug dose as initial condition +end + +%% ------------------------------------------------------------------------ +% Prepare main outputs 'Y_OUT' -------------------------------------------- + +Y = abs(Y) ;% removing negatives resulting from numerical error + +Y_OUT = [ log10( Y(:,1) / Y(1,1) )... V [log fold change] + , sum(Y(:,2:10),2) / sum(Y(1,2:10)) ];% E+B [fold change] + +%% ------------------------------------------------------------------------ +% Add columns to 'Y_OUT' (if requested) ----------------------------------- + +if isempty(fullOut) ; return ; end + +Y_OUT = [ Y_OUT, Y, zeros(Pieces(end),27) ] ; + +V = Y(:,01) ;% SIV virions [#/無] +E0 = Y(:,02) ;% resting specific CD8+ T cells [#/無] +EA = sum( Y(:,03:10) ,2) ;% active specific CD8+ T cells [#/無] +Eaf = Y(:,10) ;% last generation of EA +Eap = EA - Eaf ;% EA that is proliferating +F = terms(Y') ;% modifiers for 'p0','a0','s0','g0' + +p = F(11) ;% E0 proliferation rate [/d] +a = F(12) ;% E0 activation rate [/d] +s = F(13) ;% regulation generation rate [/d] +g = F(14) ;% EA killing rate [#/無-d] + +Y_OUT(:,17) = F(:,01) ;% N803 effect [] + +Y_OUT(:,18) = (E0+EA) ;% CD8+ T cells [#/無] +Y_OUT(:,19) = EA ;% EA +Y_OUT(:,20) = (EA)./(E0+EA) ;% activated frequency in E +Y_OUT(:,21) = log10( EA ) ;% EA (log) + +Y_OUT(:,22) = p.*E0 ;% E0 proliferation term [#/無/d] +Y_OUT(:,23) = d.*E0 ;% E0 death term [#/無/d] +Y_OUT(:,24) = p ;% E0 proliferation rate [/d] +Y_OUT(:,25) = F(:,02) ;% C multiplier for p* (log) +Y_OUT(:,26) = F(:,05) ;% R multiplier for p* (log) +Y_OUT(:,27) = F(:,08) ;% density-dep multiplier for p* (log) + +Y_OUT(:,28) = a.*E0 ;% E0 activation term [#/無/d] +Y_OUT(:,29) = pA*Eap ;% EA proliferation term [#/無/d] +Y_OUT(:,30) = dA*EA ;% EA death term [#/無/d] +Y_OUT(:,31) = m*Eaf ;% EA reversion term [#/無/d] +Y_OUT(:,32) = a ;% E0 activation rate [/d] +Y_OUT(:,33) = F(:,03) ;% C multiplier for a* +Y_OUT(:,34) = F(:,06) ;% R multiplier for a* +Y_OUT(:,35) = F(:,09) ;% V multiplier for a* + +Y_OUT(:,36) = q*V ;% V growth term [#/無/d] +Y_OUT(:,37) = g.*EA.*V ;% E killing term [#/無/d] +Y_OUT(:,38) = g.*EA ;% E killing rate [/d] +Y_OUT(:,39) = F(:,07) ;% R multiplier for g* + +Y_OUT(:,40) = s ;% R generation +Y_OUT(:,41) = 1 + F(:,04)./F(:,10) ;% C multiplier for s* +Y_OUT(:,42) = F(:,10) ;% V multiplier for s* + +% Converting most outputs to log-value +logID = [ 3:12 , 15:16 , 22:42 ] ; +Y_OUT(:,logID) = log10( Y_OUT(:,logID) ) ; + +% Converting virus to #/mL +Y_OUT(:,3) = Y_OUT(:,3) + 3 ; + +%% ======================================================================== +% NESTED FUNCTIONS +% ======================================================================== + +function dY = system(~,Y) + + if toc(solvertime) >= timeOut + error('ODE solver stalled.') + end + + V = Y(01) ;% SIV virions [#/無] + E0 = Y(02) ;% resting specific CD8+ T cells [#/無] + EA = Y(03:10) ;% active specific CD8+ T cells [#/無] (E1 to E8) + 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') ;% modifiers for 'p0','a0','s0','g0' + + p = F(11) ;% E0 proliferation rate [/d] + a = F(12) ;% E0 activation rate [/d] + s = F(13) ;% regulation generation rate [/d] + g = F(14) ;% EA killing rate [#/無-d] + + dY = Y ;% dY/dt + + % calculate rates of change for V,E,B,X,C,R + dY(01) = q*V - g*sum(EA)*V ;% dV/dt + dY(02) = p*E0 - a*E0 - d*E0 + m*EA(8) ;% dE0/dt + dY(03) = + 2*a*E0 - dA*EA(1) - pA*EA(1) ;% dE1/dt + dY(4:9)= 2*pA*EA(1:6) - dA*EA(2:7) - pA*EA(2:7) ;% dE2-7/dt + dY(10) = 2*pA*EA(7) - dA*EA(8) - m*EA(8) ;% dE8/dt + dY(11) = - ka*X ;% dX/dt + dY(12) = ka*X/vd - ke*C ;% dC/dt + dY(13) = s - dR*R(1) ;% dR1/dt + dY(14) = dR*R(1) - dR*R(2) ;% dR2/dt + +end + +function J = jacob(~,Y) + + V = Y(01) ;% SIV virions [#/無] + E0 = Y(02) ;% resting specific CD8+ T cells [#/無] + EA = sum(Y(03:10)) ;% active specific CD8+ T cells [#/無] (E1 to E8) + C = Y(12) ;% N803 plasma concentration [pM] + F = terms(Y') ;% modifiers for 'p0','a0','s0','g0' + + p = F(11) ;% E0 proliferation rate [/d] + a = F(12) ;% E0 activation rate [/d] + g = F(14) ;% EA killing rate [#/無-d] + + J = zeros(14) ; + + % Partial derivatives of g,p,a,s w.r.t. V,E,C,R + dpdR = - p * p2 * F(05) ; + dadR = - a * a2 * F(06) ; + dgdR = - g * g2 * F(07) ; + dpdE = - p * F(08) / E50 ; + dadV = a0 * V50E / (V50E+V)^2 * F(03) * F(06) ; + dsdV = s0 * V50R / (V50R+V)^2 ; + dOdC = C50 / (C50+C)^2 ; + dpdC = p0 * F(08) * p1 * dOdC * F(05) ; + dadC = a0 * F(09) * a1 * dOdC * F(06) ; + dsdC = s0 * s1 * dOdC ; + + % Partial derivatives of f1 = dV/dt + J(01,01 ) = q - g*EA ; + J(01,02:10) = - g * V ; + J(01,14 ) = - dgdR*EA * V ; + + % Partial derivatives of f2 = dE0/dt + J(02,01 ) = - dadV * E0 ; + J(02,02 ) = p + dpdE*E0 - d - a ; + J(02,03:09) = dpdE*E0 ; + J(02,10 ) = dpdE*E0 + m ; + J(02,12 ) = ( dpdC - dadC ) * E0 ; + J(02,14 ) = ( dpdR - dadR ) * E0 ; + + % Partial derivatives of f3 = dE1/dt (save f3/dE1) + J(03,01 ) = 2*dadV * E0 ; + J(03,02 ) = 2*a ; + J(03,12 ) = 2*dadC * E0 ; + J(03,14 ) = 2*dadR * E0 ; + + % Partial derivatives of f4-10 = dE2-8/dt + J(31:15:121) = - dA - pA ; + J(32:15:122) = 2*pA ; + J(10,10) = - dA - m ; + + % Partial derivatives of f11 = dX/dt + J(11,11) = -ka ; + + % Partial derivatives of f12 = dC/dt + J(12,11) = ka/vd ; + J(12,12) = -ke ; + + % Partial derivatives of f13 = dR1/dt + J(13,01) = dsdV ; + J(13,12) = dsdC ; + J(13,13) = -dR ; + + % Partial derivatives of f14 = dR2/dt + J(14,13) = dR ; + J(14,14) = -dR ; +end + +function F = terms(Y) + + V = Y(:,01) ;% SIV virions [#/無] + E = sum(Y(:,02:10),2) ;% CD8+ T cells [#/無] + C = Y(:,12) ;% N803 plasma concentration [pM] + R = Y(:,14) ;% immune regulation [] + + % percolate vector for storing terms + F = zeros(size(Y,1),14) ; + + F(:,01) = C./ (C50+C) ;% N803 effect [] + + F(:,02) = 1 + p1*F(:,1) ;% E0 proliferation stimulation [] + F(:,03) = 1 + a1*F(:,1) ;% E0 activation stimulation [] + F(:,04) = s1*F(:,1) ;% N803-induced R generation [] + + F(:,05) = 1./ (1+p2*R) ;% E0 proliferation regulation [] + F(:,06) = 1./ (1+a2*R) ;% E0 activation regulation [] + F(:,07) = 1./ (1+g2*R) ;% EA killing regulation [] + + F(:,08) = E50./(E50+E) ;% E0 proliferation density dependence [] + F(:,09) = V./ (V50E+V) ;% E0 activation antigen dependence [] + F(:,10) = V./ (V50R+V) ;% R generation antigen dependence [] + + F(:,11) = p0* F(:,08).* F(:,02).* F(:,05) ;% E0 proliferation [/d] + F(:,12) = a0* F(:,09).* F(:,03).* F(:,06) ;% E0 activation rate [/d] + F(:,13) = s0* (F(:,10)+F(:,04)) ;% R generation rate [/d] + F(:,14) = g0* F(:,07) ;% EA killing rate [#/無-d] + +end + +end \ No newline at end of file diff --git a/N803_shared_noN.m b/N803_shared_noN.m index 595f4dc..d30f462 100644 --- a/N803_shared_noN.m +++ b/N803_shared_noN.m @@ -1,23 +1,20 @@ -%% N803_shared_noN.m - solves model for 2 cohorts with shared parameters -% -% Edits of 'N803_shared.m' to zero out non-SIV-specific CD8+ T cells noted !! -% -%#ok<*NASGU> suppressing 'variable is not used' +%% N803_shared_noN.m - solves model 2S for 2 cohorts with shared parameters % % /--------------------------------------------------------------\ -% | Date: 12/25/2022 | +% | Date: 12/27/2022 | % | Author: Jonathan Cody | % | Affiliation: Purdue University | % | Weldon School of Biomedical Engineering | % | Pienaar Computational Systems Pharmacology Lab | % \--------------------------------------------------------------/ % +% Model 2S is a simplified version of Model 2 that combines SIV-specific +% and non-SIV-specific CD8+ T cells +% % Nomenclature: V = SIV virions [#/無] % T8 = total CD8+ T cells [#/無] -% S0 = resting SIV-specific CD8+ T cells [#/無] -% Sa = active SIV-specific CD8+ T cells [#/無] -% N0 = resting non-SIV-specific CD8+ T cells [#/無] -% Na = active non-SIV-specific CD8+ T cells [#/無] +% E0 = resting CD8+ T cells [#/無] +% EA = active CD8+ T cells [#/無] % X = N803 at absorption site [pmol/kg] % C = N803 plasma concentration [pM] % R = regulation [] (dimensionless quantity) @@ -77,159 +74,122 @@ P_Cohort = cell(1,2) ;% cell for storing parmeters % Rename inputed parameters ----------------------------------------------- -Vi(1) = AllPars(01) ;% V initial value [log(#/mL)] (cohort 1) -Vi(2) = AllPars(02) ;% V initial value [log(#/mL)] (cohort 2) -SNi(1) = AllPars(03) ;% S+N initial value [#/無] (cohort 1) -SNi(2) = AllPars(04) ;% S+N initial value [#/無] (cohort 2) -fS(1) = AllPars(05) ;% initial frequency: S/(S+N) (cohort 1) -fSn = AllPars(06) ;% normalized S/(S+N) (co 2) -fS(2) = fS(1) ;% +(0.3-fS(1))*fSn ;% initial frequency: S/(S+N) (cohort 2) !! - -q = AllPars(07) ;% V growth rate (if S+N were absent) [/d] -psi = AllPars(08) ;% ratio of base killing rates gN0/gS0 -V50S = AllPars(09) ;% 50% viral stimulation saturation for S [#/mL] -V50N = AllPars(10) ;% 50% viral stimulation saturation for N [#/mL] -mSn = AllPars(11) ;% normalized Sa reversion rate constant [] -mNn = AllPars(12) ;% normalized Na reversion rate constant [] - -SN50 = AllPars(13) ;% 50% S+N proliferation saturation [#/無] -p0n = AllPars(14) ;% nomalized S0/N0 proliferation rate constant [/d] -pS = AllPars(15) ;% Sa proliferation rate constant [/d] -pN = AllPars(16) ;% Na proliferation rate constant [/d] -d = AllPars(17) ;% S0/N0 death rate constant [/d] -dA = AllPars(18) ;% Sa/Na death rate constant [/d] - -Xi = AllPars(19) ;% X initial condition [pmol/kg] -ka = AllPars(20) ;% N803 absorption rate constant [/d] -ke = AllPars(21) ;% N803 elimination rate constant [/d] -vd = AllPars(22) ;% N803 'volume of distribution'/'bioavailability' [L/kg] -C50 = AllPars(23) ;% 50% N803 stimulation concentration [pM] (Cohort 1) -pm = AllPars(24) ;% S0/N0 maximum proliferation rate [] -aS1 = AllPars(25) ;% S activation stimulation factor [] -aN1 = AllPars(26) ;% N activation stimulation factor [] - -dR = AllPars(27) ;% R decay rate constant [/d] -sig = AllPars(28) ;% ratio of initial regulation due to N/S -p2 = AllPars(29) ;% S0/N0 proliferation regulation factor [] -gN2 = AllPars(30) ;% N killing regulation factor [] (cohort 1) +Vi(1) = AllPars(01) ;% V initial value [log(#/mL)] (cohort 1) +Vi(2) = AllPars(02) ;% V initial value [log(#/mL)] (cohort 2) +Ei(1) = AllPars(03) ;% E+B initial value [#/無] (cohort 1) +Ei(2) = AllPars(04) ;% E+B initial value [#/無] (cohort 2) + +q = AllPars(05) ;% V growth rate (if E were absent) [/d] +V50E = AllPars(06) ;% 50% viral stimulation saturation for E [#/mL] +V50R = AllPars(07) ;% 50% viral stimulation saturation for B [#/mL] +m = AllPars(08) ;% EA reversion rate constant [/d] + +E50 = AllPars(09) ;% 50% E proliferation saturation [#/無] +p0n = AllPars(10) ;% E0 nomalized proliferation rate constant [/d] +pA = AllPars(11) ;% 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 [pmol/kg] +ka = AllPars(15) ;% N803 absorption rate constant [/d] +ke = AllPars(16) ;% N803 elimination rate constant [/d] +vd = AllPars(17) ;% N803 'volume of distribution'/'bioavailability' [L/kg] +C50 = AllPars(18) ;% 50% N803 stimulation concentration [pM] (Cohort 1) +pm = AllPars(19) ;% E0 maximum proliferation rate [] +a1 = AllPars(20) ;% E activation stimulation factor [] +s1 = AllPars(21) ;% R generation stimulation factor [] + +dR = AllPars(22) ;% R decay rate constant [/d] +p2 = AllPars(23) ;% E0 proliferation regulation factor [] +g2 = AllPars(24) ;% EA killing regulation factor [] %% ------------------------------------------------------------------------ % Calculate some initial conditions & parameters -------------------------- Vi = 10.^(Vi - 3) ;% converting V initial value [#/無] -V50S = V50S/1000 ;% 50% viral stimulation saturation for S [#/無] -V50N = V50N/1000 ;% 50% viral stimulation saturation for N [#/無] -YS = Vi ./ (Vi+ V50S) ;% initial V/(V50S+V) [cohort 1,2] -YN = Vi ./ (Vi+ V50N) ;% initial V/(V50N+V) [cohort 1,2] - -% calculate initial R, and sS,sN -z = YS + sig*YN ; -Ri = [ 1 , (z(2)/z(1)) ] ;% initial regulation [cohort 1,2] -R = Ri(2) ;% initial regulation (cohort 2) -sS = dR / ( YS(1) + sig*YN(1) ) ;% R generation due to S0 activation [/d] -sN = sig*sS ;% R generation due to N0 activation [/d] - -% restrict mS and mN such that initial activation aS and aN are positive -US = 2*(2*pS/(pS+dA))^7 ; -UN = 2* 2*pN/(pN+dA) ; -mS = mSn*dA/(US-1) ;% Sa reversion rate constant [/d] -mN = mNn*dA/(UN-1) ;% Na reversion rate constant [/d] - -% restrict p0 to ensure same sign for S0/SA and for N0/NA -p0_max = min( d*(1+p2*Ri).*(SN50+SNi)/SN50 ) ;% maximum value for p0 -p0 = p0n * p0_max ;% S0/N0 prolif rate constant [/d] -p1 = pm/p0 ;% S0/N0 proliferation stimulation factor -pi = p0*SN50./(SN50+SNi)./(1+p2*Ri) ;% initial S0/N0 prolif rate [/d] - -% calculate aS0,aN0 and aS2,aN2 -aSi = (d-pi)/( mS*US/(dA+mS) - 1 ) ;% initial S0 activation rate [/d] -aNi = (d-pi)/( mN*UN/(dA+mN) - 1 ) ;% initial N0 activation rate [/d] -z = aSi./YS ; -aS2 = ( z(1)-z(2) ) / ( R*z(2)-z(1) ) ;% S0 activation rate constant [/d] -aS0 = aSi(1) / YS(1) * (1+aS2) ;% S activation regulation factor [] -z = aNi./YN ; -aN2 = ( z(1)-z(2) ) / ( R*z(2)-z(1) ) ;% N0 activation rate constant [/d] -aN0 = aNi(1) / YN(1) * (1+aN2) ;% N activation regulation factor [] - -% solve for initial ratios below (based on active steady-state) -ZS = US/(mS+dA) ; +V50E = V50E/1000 ;% 50% viral stimulation saturation for E [#/無] +V50R = V50R/1000 ;% 50% viral stimulation saturation for B [#/無] +YE = Vi ./ (Vi+ V50E) ;% initial V/(V50E+V) [cohort 1,2] +YR = Vi ./ (Vi+ V50R) ;% initial V/(V50B+V) [cohort 1,2] + +% calculate initial R (normalized to Cohort 1) and s +s0 = dR / YR(1) ;% R generation rate [/d] +Ri = [ 1 , s0*YR(2)/dR ] ;% initial regulation [cohort 1,2] + +% restrict p0 to ensure same sign for E0/EA and for B0/BA +p0_max = min( d*(1+p2*Ri).*(E50+Ei)/E50 ) ;% maximum value for p0 +p0 = p0n * p0_max ;% E0/B0 prolif rate constant [/d] +p1 = pm/p0 ;% E0/B0 proliferation stimulation factor +pi = p0*E50./(E50+Ei)./(1+p2*Ri) ;% initial E0/B0 prolif rate [/d] + +% calculate initial ratio: E8/E0/a +U = 2/(m+dA) * (2*pA/(pA+dA))^7 ; + +% calculate a0, a2 +ai = (d-pi)/( m*U - 1 ) ;% initial E0 activation rate [/d] +F = ai./YE ; +a2 = ( F(1)-F(2) ) / ( Ri(2)*F(2)-F(1) );% E0 activation rate constant [/d] +a0 = ai(1) / YE(1) * (1+a2) ;% E0 activation regulation factor [] + +% calculate initial ratio of EA/E8 +W = 1 ; for i = 1:7 - ZS = ZS + 2*(2*pS)^(i-1)/(pS+dA)^i ;% SAi/aSi/S0i + W = W + (pA+dA)^(i-1) / (2*pA)^i * (m+dA) ;% EAi/E8i end -ZN = UN/(mN+dA) + 2/(pN+dA) ;% NAi/aNi/N0i - -% solve for initial values of S0,Sa,N0,Na each cohort -Si = SNi.*fS ;% initial S -Ni = SNi.*(1-fS) ;% initial N -S0 = Si./(1+ZS*aSi) ;% initial S0 -N0 = Ni./(1+ZN*aNi) ;% initial N0 -SA = S0.*(ZS*aSi) ;% initial SA -NA = N0.*(ZN*aNi) ;% initial NA - -% calculate gS0,gN0 and gS2 -beta = psi * NA ./ (1+gN2*Ri) ; z = beta(1) - beta(2) ; -a = z*R ; b = SA(1)*R - SA(2) + z*(1+R) ; c = SA(1) - SA(2) + z ; -gS2 = -c/b ;% ( -b + sqrt(b^2 - 4*a*c) ) / (2*a) ;% S killing regulation factor [] !! -gS0 = q / ( SA(1)/(1+gS2) + beta(1) ) ;% Sa killing rate constant [無/#-d] -gN0 = psi * gS0 ;% Na killing rate constant [無/#-d] + +% calculate initial values of E0,EA for each cohort +E0 = Ei.*(m/W)./ ( (m/W) - pi + d + ai ) ;% initial E0 +EA = Ei - E0 ;% initial EA + +% calculate g0 +g0 = q * (1+g2) / EA(1) ;% EA killing rate constant [無/#-d] %% Do for each cohort (NOT indenting loop) ================================ for c = 1:2 -% solve for initial S1-8 and N1-2 -S = zeros(1,8) ;% initial S1-S8 -S(1) = 2*aSi(c)*S0(c) / (dA+pS) ;% S1 +% solve for initial E1-8 and B1-2 +E = zeros(1,8) ;% initial E1-E8 +E(1) = 2*ai(c)*E0(c) / (pA+dA) ;% E1 for i = 2:7 - S(i) = 2*pS*S(i-1) / (dA+pS) ;% S2 to S7 + E(i) = 2*pA*E(i-1) / (pA+dA) ;% E2 to E7 end -S(8) = 2*pS*S(7) / (dA+mS) ;% S8 - -N(1) = 2*aNi(c)*N0(c) / (dA+pN) ;% N1 -N(2) = 2*pN*N(1) / (dA+mN) ;% N2 +E(8) = 2*pA*E(7) / (m+dA) ;% E8 %% ------------------------------------------------------------------------ % Prepare parameter and initial value vectors and call 'N803_model_2' ----- -Pars(01) = q ;% V growth rate (if S+N were absent) [/d] -Pars(02) = gS0 ;% Sa killing rate constant [無/#-d] -Pars(03) = gN0 ;% Na killing rate constant [無/#-d] - -Pars(04) = V50S ;% 50% viral stimulation saturation for S [#/無] -Pars(05) = V50N ;% 50% viral stimulation saturation for N [#/無] -Pars(06) = aS0 ;% S0 activation rate constant [/d] -Pars(07) = aN0 ;% N0 activation rate constant [/d] -Pars(08) = mS ;% Sa reversion rate constant [/d] -Pars(09) = mN ;% Na reversion rate constant [/d] - -Pars(10) = SN50 ;% 50% S+N proliferation saturation [#/無] -Pars(11) = p0 ;% S0/N0 proliferation rate constant [/d] -Pars(12) = pS ;% Sa proliferation rate constant [/d] -Pars(13) = pN ;% Na proliferation rate constant [/d] -Pars(14) = d ;% S0/N0 death rate constant [/d] -Pars(15) = dA ;% Sa/Na death rate constant [/d] - -Pars(16) = Xi ;% X initial condition [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 stimulation concentration [pM] -Pars(21) = p1 ;% S0/N0 proliferation stimulation factor [] -Pars(22) = aS1 ;% S activation stimulation factor [] -Pars(23) = aN1 ;% N activation stimulation factor [] - -Pars(24) = sS ;% R generation due to S0 activation [/d] -Pars(25) = sN ;% R generation due to N0 activation [/d] -Pars(26) = dR ;% R decay rate constant [/d] -Pars(27) = gS2 ;% S killing regulation factor [] -Pars(28) = gN2 ;% N killing regulation factor [] -Pars(29) = p2 ;% S0/N0 proliferation regulation factor [] -Pars(30) = aS2 ;% S activation regulation factor [] -Pars(31) = aN2 ;% N activation regulation factor [] - -% V S0-8 N0-2 X C R initial values -Yic = [ Vi(c) S0(c) S N0(c) N 0 0 Ri(c) Ri(c) ] ; - -% if any( [ Pars([1:6 8 10:30 ]) Yic ] < 0 ) % !! +Pars(01) = q ;% V growth rate (if E+B were absent) [/d] +Pars(02) = g0 ;% EA killing rate constant [無/#-d] +Pars(03) = V50E ;% 50% viral stimulation saturation for E [#/無] +Pars(04) = V50R ;% 50% viral stimulation saturation for R [#/無] +Pars(05) = a0 ;% E0 activation rate constant [/d] +Pars(06) = m ;% EA reversion rate constant [/d] + +Pars(07) = E50 ;% 50% E proliferation saturation [#/無] +Pars(08) = p0 ;% 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 [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 stimulation concentration [pM] +Pars(17) = p1 ;% E0 proliferation stimulation factor [] +Pars(18) = a1 ;% E0 activation stimulation factor [] +Pars(19) = s1 ;% R generation stimulation factor [] + +Pars(20) = s0 ;% R generation maximum base rate [/d] +Pars(21) = dR ;% R decay rate constant [/d] +Pars(22) = g2 ;% EA killing regulation factor [] +Pars(23) = p2 ;% E0 proliferation regulation factor [] +Pars(24) = a2 ;% E0 activation regulation factor [] + +% V E0-8 X C R initial values +Yic = [ Vi(c) E0(c) E 0 0 Ri(c) Ri(c) ] ; + +% if any( [ Pars Yic ] < 0 ) % error('Negative parameters or initial values.') % end @@ -241,7 +201,7 @@ idSol = ~ ( idLo | idHi ) ;% index of times in 'SoluTimes' to solve if RunCohort(c) == 1 - Y_COH = N803_model_2(SoluTimes(idSol),DoseTimes{c},Pars,Yic,varargin) ; + Y_COH = N803_model_2S(SoluTimes(idSol),DoseTimes{c},Pars,Yic,varargin); Y_LO = ones(sum(idLo),1)*Y_COH(1 ,:) ;% constant Y for early times Y_HI = ones(sum(idHi),1)*Y_COH(end,:) ;% constant Y for later times Y_COH = [ Y_LO ; Y_COH ; Y_HI ] ;%#ok % total 'solution' matrix diff --git a/main_shared_calibration_2022_08_18.m b/main_shared_calibration_2022_08_18.m new file mode 100644 index 0000000..fd3a098 --- /dev/null +++ b/main_shared_calibration_2022_08_18.m @@ -0,0 +1,158 @@ +%% 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.m'. +% ODEs are solved in 'N803_model_2.m'. +% Files within 'Analyzer' folder perform MSLS and MCMC. + +%% ======================================================================== +% INPUTS for 'calibrator.m' +% ======================================================================== + +% 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 censoring and normalization to data (see 'prepper' function) +[Xp,Yp,Sid,Cid,IVs] = ... + prepper(X,Y, {2,[],2,[]}, {2 ,[],2 ,[]}, {0,0,0,0}, [1 0 1 0]) ; +[X,Y] = prepper(X,Y, {2,[],2,[]}, {[],[],[],[]}, {0,0,0,0}, [1 0 1 0]) ; + +% GuessBOUNDS = array for boundaries of initial parameter guesses ----- +GuessBOUNDS(:,01) = [1 1]*mean(cat(1,IVs{1}{:})) ;% V iv [log(#/mL)] (co 1) +GuessBOUNDS(:,02) = [1 1]*mean(cat(1,IVs{3}{:})) ;% V iv [log(#/mL)] (co 2) +GuessBOUNDS(:,03) = [1 1]*mean(cat(1,IVs{2}{:})) ;% S+N iv [#/無] (co 1) +GuessBOUNDS(:,04) = [1 1]*mean(cat(1,IVs{4}{:})) ;% S+N iv [#/無] (co 2) +GuessBOUNDS(:,05) = [.03 .3] ;% initial frequency: S/(S+N) (co 1) +GuessBOUNDS(:,06) = [.01 1] ;% normalized S/(S+N) (co 2) + +GuessBOUNDS(:,07) = [.01 1] ;% V growth rate (if S+N were absent) [/d] +GuessBOUNDS(:,08) = [.01 1] ;% Na/Sa killing rate ratio [gN0/gS0] +GuessBOUNDS(:,09) = [1e4 1e6] ;% 50% viral stim saturation for S [#/mL] +GuessBOUNDS(:,10) = [1e4 1e6] ;% 50% viral stim saturation for N [#/mL] +GuessBOUNDS(:,11) = [1 10] ;% normalized Sa reversion rate constant [] +GuessBOUNDS(:,12) = [1 10] ;% normalized Na reversion rate constant [] + +GuessBOUNDS(:,13) = [300 3e3] ;% 50% S+N proliferation saturation [#/無] +GuessBOUNDS(:,14) = [0.1 1] ;% S0/N0 normalized prolif rate constant [/d] +GuessBOUNDS(:,15) = [1 4] ;% Sa proliferation rate constant [/d] +GuessBOUNDS(:,16) = [0.5 2] ;% Na proliferation rate constant [/d] +GuessBOUNDS(:,17) = [.01 .05] ;% S0/N0 death rate constant [/d] +GuessBOUNDS(:,18) = [.1 .5] ;% Sa/Na death rate constant [/d] + +GuessBOUNDS(:,19) = [1 1]*880 ;% X initial condition [pmol/kg] +GuessBOUNDS(:,20) = [1 1]*0.8 ;% N803 absorption rate constant [/d] +GuessBOUNDS(:,21) = [1 1]*2.1 ;% N803 elimination rate constant [/d] +GuessBOUNDS(:,22) = [1 1]*1.3 ;% N803 'vol of dist'/'bioavail' [L/kg] +GuessBOUNDS(:,23) = [1 1e3] ;% 50% N803 stimulation concentration [pM] +GuessBOUNDS(:,24) = [.1 2] ;% S0/N0 maximum prolif rate [] +GuessBOUNDS(:,25) = [0.1 1e3] ;% S activation stimulation factor [] +GuessBOUNDS(:,26) = [0.01 100] ;% N activation stimulation factor [] + +GuessBOUNDS(:,27) = [0.05 1] ;% R decay rate constant [/d] +GuessBOUNDS(:,28) = [0.1 100] ;% sN/sS regulation generation ratio +GuessBOUNDS(:,29) = [0.01 100] ;% S0/N0 prolif reg factor [] +GuessBOUNDS(:,30) = [0.01 100] ;% N killing regulation factor [] + +% Inputs for model function and cost function ----------------------------- +Xeval = min(cat(1,X{:})) : 1 : max(cat(1,X{:})) ;% x-values (days) +Doses = { 7*[0:3 5:8 37:40] ... x-values for doses (Cohort 1) + , 7*(0:2:4) };% x-values for doses (Cohort 2) +Skips = { [-7 max(cat(1,X{1:2}))] ... [min max] time to solve (Cohort 1) + , [-7 max(cat(1,X{3:4}))] };% [min max] time to solve (Cohort 2) + +% modelFunction = handle for calling model function ------------------- +modelFunction = @(t,P) N803_shared(t,Doses,P,Skips,[],... +'warnOff',1,'AbsTol',1e-3,'RelTol',1e-2,'odeSolver','23s') ; + +% FitBOUNDS = array for boundaries of parameter fitting ------------------- +FitBOUNDS = GuessBOUNDS ; + +% costFunction = handle for calling cost function ------------------------- +costFunction = @(P) costFun_NLL(P,modelFunction,X,Y,'Xevaluate',Xeval) ; + +% costFunction(GuessBOUNDS(1,:)) ;% uncomment to debug your costFunction + +% filePrefix = string to prefix to file names ----------------------------- +filePrefix = 'N803_shared' ; + +% timeLimit = number of minutes allowed for a single optimization run ----- +timeLimit = 30 ; + +% numberGuesses = number of initial parameter guesses --------------------- +numberGuesses = 2000 ; + +% numberWorkers = number of parallel processes (enter [] or 0 for serial) +numberWorkers = 0 ; + +%% ------------------------------------------------------------------------ +% Optional inputs for results summary figures ('calibrator_plot.m') +% ------------------------------------------------------------------------ + +plotSetup.Xdata = Xp ; +plotSetup.Ydata = Yp ; +plotSetup.CensorID = Cid ; +plotSetup.Fits2Plot = 1:82 ; +plotSetup.modelFunction = modelFunction ; +plotSetup.nFits2Sum = 82 ; +plotSetup.SubjectID = Sid ; +plotSetup.XLabels = {{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'}} ; +plotSetup.YLabels = {{'vRNA copy','log fold \Delta'} ... + ,{'CD8^{+} T cell','fold \Delta'}... + ,{'vRNA copy','log fold \Delta'} ... + ,{'CD8^{+} T cell','fold \Delta'}} ; +plotSetup.XRANGE = 7*[ -4 41 ; -4 41 ; -4 6 ; -4 6 ] ; +plotSetup.YRANGE = [ -3 1 ; 0 6 ; -2 2 ; 0 6 ] ; +plotSetup.PlotGrid = [2 2] ; +plotSetup.PlotOuts = [1 3 2 4] ; +plotSetup.plotPars = 0 ; +plotSetup.xPrecise = 1 ; + +%% ======================================================================== +% 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:9:82 ; + +% MCMC sample size -------------------------------------------------------- +numSamples = 1e6 ; + +% save every 'saveEach' iterations ---------------------------------------- +saveEach = 1000 ; + +% Perform uncertainty quantification using 'MCMCsampler.m' +ResultsFile = [filePrefix '_Results'] ; +MCMCsampler(ResultsFile,InitialRanks,numSamples,saveEach) ; diff --git a/main_shared_noN_calibration_2022_11_02.m b/main_shared_noN_calibration_2022_11_02.m new file mode 100644 index 0000000..c84a889 --- /dev/null +++ b/main_shared_noN_calibration_2022_11_02.m @@ -0,0 +1,155 @@ +%% 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) + +% This is a repeat of 'main_shared_calibration_2022_08_18.m' without +% non-SIV-specific CD8+ T cells + +% 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_noN.m'. +% ODEs are solved in 'N803_model_2S.m'. +% Files within 'Analyzer' folder perform MSLS and MCMC. + +%% ======================================================================== +% INPUTS for 'calibrator.m' +% ======================================================================== + +% 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 censoring and normalization to data (see 'prepper' function) +[Xp,Yp,Sid,Cid,IVs] = ... + prepper(X,Y, {2,[],2,[]}, {2 ,[],2 ,[]}, {0,0,0,0}, [1 0 1 0]) ; +[X,Y] = prepper(X,Y, {2,[],2,[]}, {[],[],[],[]}, {0,0,0,0}, [1 0 1 0]) ; + +% GuessBOUNDS = array for boundaries of initial parameter guesses ----- +GuessBOUNDS(:,01) = [1 1]*mean(cat(1,IVs{1}{:})) ;% V iv [log(#/mL)] (co 1) +GuessBOUNDS(:,02) = [1 1]*mean(cat(1,IVs{3}{:})) ;% V iv [log(#/mL)] (co 2) +GuessBOUNDS(:,03) = [1 1]*mean(cat(1,IVs{2}{:})) ;% E iv [#/無] (co 1) +GuessBOUNDS(:,04) = [1 1]*mean(cat(1,IVs{4}{:})) ;% E iv [#/無] (co 2) + +GuessBOUNDS(:,05) = [.01 1] ;% V growth rate (if E were absent) [/d] +GuessBOUNDS(:,06) = [1e4 1e6] ;% 50% viral stim saturation for E [#/mL] +GuessBOUNDS(:,07) = [1e4 1e6] ;% 50% viral stim saturation for R [#/mL] +GuessBOUNDS(:,08) = [.005 .05] ;% EA reversion rate constant [] + +GuessBOUNDS(:,09) = [300 3e3] ;% 50% E proliferation saturation [#/無] +GuessBOUNDS(:,10) = [0.1 1] ;% E0 normalized prolif rate constant [/d] +GuessBOUNDS(:,11) = [1 4] ;% EA proliferation 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 [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 maximum prolif 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] ;% E0 proliferation regulation factor [] +GuessBOUNDS(:,24) = [0.01 100] ;% EA killing regulation factor [] + +% Inputs for model function and cost function ----------------------------- +Xeval = min(cat(1,X{:})) : 1 : max(cat(1,X{:})) ;% x-values (days) +Doses = { 7*[0:3 5:8 37:40] ... x-values for doses (Cohort 1) + , 7*(0:2:4) };% x-values for doses (Cohort 2) +Skips = { [-7 max(cat(1,X{1:2}))] ... [min max] time to solve (Cohort 1) + , [-7 max(cat(1,X{3:4}))] };% [min max] time to solve (Cohort 2) + +% modelFunction = handle for calling model function ------------------- +modelFunction = @(t,P) N803_shared_noN(t,Doses,P,Skips,[],... +'warnOff',1,'AbsTol',1e-3,'RelTol',1e-2,'odeSolver','23s') ; + +% FitBOUNDS = array for boundaries of parameter fitting ------------------- +FitBOUNDS = GuessBOUNDS ; + +% costFunction = handle for calling cost function ------------------------- +costFunction = @(P) costFun_NLL(P,modelFunction,X,Y,'Xevaluate',Xeval) ; + +% costFunction(GuessBOUNDS(1,:)) ;% uncomment to debug your costFunction + +% filePrefix = string to prefix to file names ----------------------------- +filePrefix = 'N803_shared_2S' ; + +% timeLimit = number of minutes allowed for a single optimization run ----- +timeLimit = 30 ; + +% numberGuesses = number of initial parameter guesses --------------------- +numberGuesses = 2000 ; + +% numberWorkers = number of parallel processes (enter [] or 0 for serial) +numberWorkers = 12 ; + +%% ------------------------------------------------------------------------ +% Optional inputs for results summary figures ('calibrator_plot.m') +% ------------------------------------------------------------------------ + +plotSetup.Xdata = Xp ; +plotSetup.Ydata = Yp ; +plotSetup.CensorID = Cid ; +plotSetup.Fits2Plot = 1:82 ; +plotSetup.modelFunction = modelFunction ; +plotSetup.nFits2Sum = 82 ; +plotSetup.SubjectID = Sid ; +plotSetup.XLabels = {{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'}} ; +plotSetup.YLabels = {{'vRNA copy','log fold \Delta'} ... + ,{'CD8^{+} T cell','fold \Delta'}... + ,{'vRNA copy','log fold \Delta'} ... + ,{'CD8^{+} T cell','fold \Delta'}} ; +plotSetup.XRANGE = 7*[ -4 41 ; -4 41 ; -4 6 ; -4 6 ] ; +plotSetup.YRANGE = [ -3 1 ; 0 6 ; -2 2 ; 0 6 ] ; +plotSetup.PlotGrid = [2 2] ; +plotSetup.PlotOuts = [1 3 2 4] ; +plotSetup.plotPars = 0 ; +plotSetup.xPrecise = 1 ; + +%% ======================================================================== +% 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:9:82 ; + +% MCMC sample size -------------------------------------------------------- +numSamples = 1e6 ; + +% save every 'saveEach' iterations ---------------------------------------- +saveEach = 1000 ; + +% Perform uncertainty quantification using 'MCMCsampler.m' +ResultsFile = [filePrefix '_Results'] ; +MCMCsampler(ResultsFile,InitialRanks,numSamples,saveEach) ; diff --git a/main_singles_calibration_2022_09_25.m b/main_singles_calibration_2022_09_25.m new file mode 100644 index 0000000..c4eab89 --- /dev/null +++ b/main_singles_calibration_2022_09_25.m @@ -0,0 +1,163 @@ +%% 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 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 and MCMC. + +%% ======================================================================== +% INPUTS for 'calibrator.m' +% ======================================================================== + +% GuessBOUNDS = array for boundaries of initial parameter guesses --------- +GuessBOUNDS(:,03) = [.03 .3] ;% initial frequency: S/(S+N) +GuessBOUNDS(:,04) = [.009 .9] ;% initial frequency: Sa/S +GuessBOUNDS(:,05) = [.01 1] ;% V growth rate (if S+N were absent) [/d] +GuessBOUNDS(:,06) = [.01 1] ;% Na/Sa killing rate ratio [gN0/gS0] +GuessBOUNDS(:,07) = [1e4 1e6] ;% 50% viral stim saturation for S [#/mL] +GuessBOUNDS(:,08) = [1e4 1e6] ;% 50% viral stim saturation for N [#/mL] +GuessBOUNDS(:,09) = [1 10] ;% normalized Sa reversion rate constant [] +GuessBOUNDS(:,10) = [1 10] ;% normalized Na reversion rate constant [] +GuessBOUNDS(:,11) = [300 3e3] ;% 50% S+N proliferation saturation [#/無] +GuessBOUNDS(:,12) = [1 4] ;% Sa proliferation rate constant [/d] +GuessBOUNDS(:,13) = [0.5 2] ;% Na proliferation rate constant [/d] +GuessBOUNDS(:,14) = [.01 .05] ;% S0/N0 death rate constant [/d] +GuessBOUNDS(:,15) = [.1 .5] ;% Sa/Na death rate constant [/d] +GuessBOUNDS(:,16) = [1 1]*880 ;% X initial condition [pmol/kg] +GuessBOUNDS(:,17) = [1 1]*0.8 ;% N803 absorption rate constant [/d] +GuessBOUNDS(:,18) = [1 1]*2.1 ;% N803 elimination rate constant [/d] +GuessBOUNDS(:,19) = [1 1]*1.3 ;% N803 'vol of dist'/'bioavail' [L/kg] +GuessBOUNDS(:,20) = [1 1e3] ;% 50% N803 stimulation concentration [pM] +GuessBOUNDS(:,21) = [.1 2] ;% S0/N0 maximum proliferation rate [] +GuessBOUNDS(:,22) = [0.1 1e3] ;% S activation stimulation factor [] +GuessBOUNDS(:,23) = [0.01 100] ;% N activation stimulation factor [] +GuessBOUNDS(:,24) = [0.1 100] ;% sN/sS regulation generation rate ratio +GuessBOUNDS(:,25) = [0.05 1] ;% R decay rate constant [/d] +GuessBOUNDS(:,26) = [0.01 100] ;% initial S killing regulation [] +GuessBOUNDS(:,27) = [0.01 100] ;% initial N killing regulation [] +GuessBOUNDS(:,28) = [0.01 100] ;% initial S0/N0 proliferation regulation [] +GuessBOUNDS(:,29) = [0.01 100] ;% initial S activation regulation [] +GuessBOUNDS(:,30) = [0.01 100] ;% initial N activation regulation [] + +% timeLimit = number of minutes allowed for a single optimization run ----- +timeLimit = 30 ; + +% numberGuesses = number of initial parameter guesses --------------------- +numberGuesses = 500 ; + +% numberWorkers = number of parallel processes (enter [] or 0 for serial) +numberWorkers = 0 ; + +%% ------------------------------------------------------------------------ +% Do for each cohort... +for c = 1:2 + +% 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 censoring and normalization to data (see 'prepper' function) +[Xp,Yp,Sid,Cid,IVs] = prepper(X,Y,{2,[]},{ 2,[]},{0,0},[1 0]) ; +[X,Y] = prepper(X,Y,{2,[]},{[],[]},{0,0},[1 0]) ; + +% GuessBOUNDS = array for boundaries of initial parameter guesses --------- +GuessBOUNDS(:,01) = [1 1]*mean(cat(1,IVs{1}{:})) ;% V iv [log(#/mL)] +GuessBOUNDS(:,02) = [1 1]*mean(cat(1,IVs{2}{:})) ;% S0+Sa iv [#/無] + +% FitBOUNDS = array for boundaries of parameter fitting ------------------- +FitBOUNDS = GuessBOUNDS ; + +% Inputs for model function and cost function ----------------------------- +Xeval = min(cat(1,X{:})) : 1 : max(cat(1,X{:})) ;% x-values (days) +Doses = { 7*[0:3 5:8 37:40] ... x-values for doses (Cohort 1) + , 7*(0:2:4) };% x-values for doses (Cohort 2) +Skips = [-7 inf] ;% [min max] time point beyond which to skip solving + +% modelFunction = handle for calling model function ----------------------- +modelFunction = @(t,P) N803_single(t,Doses{c},P,Skips,... + 'warnOff',1,'AbsTol',1e-3,'RelTol',1e-2,'odeSolver','23s') ; + +% costFunction = handle for calling cost function ------------------------- +costFunction = @(P) costFun_NLL(P,modelFunction,X,Y,'Xevaluate',Xeval) ; + +% 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 = Xp ; +plotSetup.Ydata = Yp ; +plotSetup.modelFunction = modelFunction ; +plotSetup.CensorID = Cid ; +plotSetup.Fits2Plot = 1:101 ; +plotSetup.nFits2Sum = 101 ; +plotSetup.SubjectID = Sid ; +plotSetup.XLabels = {{'Days Post-Treatment'} ... + ,{'Days Post-Treatment'}} ; +plotSetup.YLabels = {{'vRNA copy','log fold \Delta'} ... + ,{'CD8^{+} T cell','fold \Delta'}} ; +switch c + case 1 + plotSetup.XRANGE = 7*[ -4 41 ; -4 41 ] ; + plotSetup.YRANGE = [ -5 3 ; 0 6 ] ; + case 2 + plotSetup.XRANGE = 7*[ -4 6 ; -4 6] ; + plotSetup.YRANGE = [ -4 4 ; 0 6 ] ; +end +plotSetup.PlotGrid = [2 1] ; +plotSetup.PlotOuts = 1:2 ; +plotSetup.plotPars = 0 ; + +%% ======================================================================== +% 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 + +%% ======================================================================== +% PERFORM MCMC +% ======================================================================== + +% InitialRanks = ranks (by NLL) of parameter sets to initialize MCMC ------ +InitialRanks = 1:25:101 ; + +% MCMC sample size -------------------------------------------------------- +numSamples = 1e6 ; + +% save every 'saveEach' iterations ---------------------------------------- +saveEach = 1000 ; + +parfor c = 1:2 + + % Perform uncertainty quantification using 'MCMCsampler.m' + ResultsFile = ['N803_cohort_' num2str(c) '_Results'] ; + MCMCsampler(ResultsFile,InitialRanks,numSamples,saveEach) ; +end \ No newline at end of file