From a4f9c07e8713ac4832a93298cbac7c400872438e Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Sun, 25 Dec 2022 09:00:08 -0600 Subject: [PATCH] Delete HIV_model_5.m --- Legacy Code/HIV_model_5.m | 470 -------------------------------------- 1 file changed, 470 deletions(-) delete mode 100644 Legacy Code/HIV_model_5.m diff --git a/Legacy Code/HIV_model_5.m b/Legacy Code/HIV_model_5.m deleted file mode 100644 index aba9672..0000000 --- a/Legacy Code/HIV_model_5.m +++ /dev/null @@ -1,470 +0,0 @@ -%% -% /--------------------------------------------------------------\ -% | Date: 03/29/2019 | -% | Author: Jonathan Cody | -% | Affiliation: Purdue University | -% | Weldon School of Biomedical Engineering | -% | Pienaar Computational Systems Pharmacology Lab | -% \--------------------------------------------------------------/ - -% Solves host-level model for virus and cell concentrations during a -% subcutaneously administered regimen of N803 - -% Nomenclature: V = virion -% T = CD4+ T cell -% R = resting CD8+ T cell (naive or memory) -% A = active-proliferating CD8+ T cell -% E = active-effector CD8+ T cell -% K = NK cell -% Q = N803 injection site quantity -% C = N803 plasma concentration -% Rc = receptor regulation due to N803 -% Sc = immune suppression due to N803 - -%% ======================================================================== -% INPUTS -% ======================================================================== - -% SolTimes = ascending vector of days at which to evaluate solution -% DoseTimes = ascending vector of days at which to administer doses -% (all elements of 'DoseTimes' must also be in 'SolTimes') -% P = vector of parameters corresponding to list below -% EquatedPars = vector of indices in 'P' used for equating parameter values -% (if empty, EquatedPars is ignored) -% timeout = exit a stalled solution after this many seconds -% (if empty, timeout = 10 [seconds]) - -% Initial Conditions -% P(01) = virion i.c. [log(/mL)] -% P(02) = CD8+ T cell (R+A+E) i.c. [/無] -% P(03) = NK cell i.c. [/無] - -% Pharmacokinetic Parameters -% P(04) = N803 injection site i.c. [pmol/kg] -% P(05) = absorption rate constant [/d] -% P(06) = P(05) * 'bioavailability'/'volume of distribution' [kg/L-d] -% P(07) = elimination rate constant [/d] - -% Pharmacodynamic Parameters -% P(08) = drug half saturation [pM] -% P(09) = Rc 1/(delay time constant) [/d] -% P(10) = Sc " -% P(11) = Rc number of delay chambers [] -% P(12) = Sc " -% P(13) = R receptor effect modifier [] -% P(14) = E " -% P(15) = K " - -% Parameters governing R,K proliferation and death -% P(16) = R death rate constant [/d] -% P(17) = K " -% P(18) = R max proliferating concentration [/無] -% P(19) = K " -% P(20) = R drug-induced proliferation factor [/d] -% P(21) = K " -% P(22) = R population suppression factor [/d] -% P(23) = K " - -% Parameters governing R activation & E,K killing -% P(24) = R activation rate constant [/d] -% P(25) = E initial killing rate (a_E * E_0) [/d] -% P(26) = K " -% P(27) = R drug-induced activation factor [/d] -% P(28) = E drug-induced killing factor [無/d] -% P(29) = K " -% P(30) = R activation suppression factor [] -% P(31) = E killing suppression factor [] -% P(32) = K " -% P(33) = R half-activation virus concentration [/無] - -% Parameters governing CD8+ T cell clonal expansion and retraction -% P(34) = A death rate constant [/d] -% P(35) = A division rate constant [/d] -% P(36) = A number of divisions [] -% P(37) = E death rate constant [/d] -% P(38) = E deactivation rate constant [/d] - -% Parameters governing CD4+ T cells -% P(39) = CD4+ T cell i.c. [/無] -% P(40) = T death rate constant [/d] -% P(41) = T max proliferating concentration [/無] -% P(42) = T infection rate constant [/無-d] -% P(43) = infected T decay rate constant [*(P(40)] - -%% ======================================================================== -% OUTPUTS -% ======================================================================== - -% Y_MAIN = model output at each time point in 'SolTimes' -% Y_MAIN(:,1) = virion [log(#/mL)] -% Y_MAIN(:,2) = CD8+ T cell (R+A+E) [/無] -% Y_MAIN(:,3) = NK cell [/無] -% Y_MAIN(:,4) = CD4+ T cell [/無] - -% Y_ALL = cell array of extra outputs (see 'output' section) - -%% ======================================================================== -% FUNCTION -% ======================================================================== - -function [Y_MAIN,Y_ALL] = HIV_model_5(SolTimes,DoseTimes,... - P,EquatedPars,timeout) - -% convert inputs to column vectors -if size(SolTimes ,1)==1 ; SolTimes = SolTimes' ; end -if size(DoseTimes ,1)==1 ; DoseTimes = DoseTimes' ; end -if size(P ,1)==1 ; P = P' ; end - -% equate parameters in 'P' based on indices in 'EquatedPars' -if ~isempty(EquatedPars) ; P = P(EquatedPars) ; end - -% set default 'timeout' to 10 seconds -if isempty(timeout) ; timeout = 10 ; end - -% declare last index in Y (solution matrix) for ... -% [(A generations) (Rc delay chambers) (Sc delay chambers)] -in = 9 + [ P(36) sum(P([11,36])) sum(P([11,12,36])) ] ; - -% build reference index in 'SolTimes' containing start time, all dose -% times, and end time (used for piecewise solution of ODEs) -PiecewiseIndex = [ 1 ; find(ismember(SolTimes,DoseTimes)) ; ... - length(SolTimes) ] ; - -% ------------------------------------------------------------------------- -% Declare initial conditions & calculated parameters ---------------------- - -V0 = 10^(P(01)-3) ;% convert V i.c. to [/無] - -% calculate i.c. for CD8+ T cell subsets (R,A,E) -CD8TCOEFS = [ ones(1,P(36)+2) ; zeros(P(36)+1,P(36)+2) ] ; -CD8TCOEFS(2,1) = P(24) * ( V0 / ( P(33) + V0 ) ) ; -CD8TCOEFS(end,end) = - P(37) - P(38) * ( P(33) / ( P(33) + V0 ) ) ; -for n = 1:P(36) - CD8TCOEFS((n+1):(n+2),n+1) = [ - P(34) - P(35) ; 2*P(35) ] ; -end -CD8Tic = CD8TCOEFS \ [ P(02) ; zeros(P(36)+1,1) ] ; - -% declare initial conditions -Yic = [ 0 ; 0 ... Q,C i.c. (pre-treatment phase) - ; 0 ; 0 ... Rc,Sc i.c. - ; CD8Tic([1,end]) ... R,E i.c. - ; P([3,39]) ... K,T i.c. - ; V0 ... V i.c. - ; CD8Tic(2:end-1) ... A i.c. - ; zeros(sum(P([11,12])),1) ];% Rc,Sc delay chambers i.c. - -% redefine P(25:26) as E,K per-capita killing rate constants -P(25:26) = P(25:26)./ Yic(6:7) ; - -% (infection rate constant) * (infected cell to virus conversion) [/無-d] -beta = ( P(40)*P(43) + sum(P(25:26).*Yic(6:7)) ) / Yic(8) ; - -% calculate [R;K;T] proliferation rate [/d] -Prolif = CD8TCOEFS(2,1) + (CD8TCOEFS(end,end)+P(37)) * Yic(6)/Yic(5) ; -Prolif = [ ( P(16) + Prolif ) * ( P(18) + Yic(5) ) / P(18) ... - ; ( P(17) ) * ( P(19) + Yic(7) ) / P(19) ... - ; ( P(40) + P(42)*V0 ) * ( P(41) + Yic(8) ) / P(41) ]; - -% ------------------------------------------------------------------------- -% Solve for each solution period ------------------------------------------ - -Y = zeros(length(SolTimes),length(Yic)) ;% percolate output matrix -Yic = Yic' ;% convert 'Yic' to row vector - -for n = 1:length(PiecewiseIndex)-1 - - solvertime = tic ;% start timer (for catching a stalled execution) - - % if a dose time is reached, add dose to N803 injection site quantity - if ismember(SolTimes(PiecewiseIndex(n)),DoseTimes) - Yic = Yic + [P(04) Yic(2:end)*0 ] ; - end - - % solve ODEs (see MATLAB help for ode15s algorithm) - [~,Y_TEMP] = ode15s(@system,... - SolTimes(PiecewiseIndex(n):PiecewiseIndex(n+1)),... - Yic) ; - - % for solutions of 2 timepoints, exclude the transition timepoints - % that ode15s adds by default - if length(SolTimes(PiecewiseIndex(n):PiecewiseIndex(n+1))) == 2 - Y_TEMP = [Y_TEMP(1,:);Y_TEMP(end,:)] ; - end - - % add 'Y_TEMP' to full solution 'Y' - Y(PiecewiseIndex(n):PiecewiseIndex(n+1),:) = Y_TEMP ; - - % update initial condition for next solution window - Yic = Y(PiecewiseIndex(n+1),:) ; -end - -% ------------------------------------------------------------------------- -% 'system' function ------------------------------------------------------- - -function dY = system(~,Y) - - % throw error if 'solvertime' exceeds 'timeout' limit - if toc(solvertime) >= timeout - error('ODE solver stalled.') - end - - dY = Y*0 ;% percolate dY/dt - - % Calculate rates of change for Rc and Sc delay chambers -------------- - - Generators = funO(Y(2))*[1;1] ;% terms that generate Rc and Sc - % (by default these are the drug effect) - % for Rc and Sc - for m = 1:2 - - % if delay chambers are specified - if P(10+m) ~= 0 - - % update delay chamber rates of change - dY(in(m)+1:in(1+m)) = ... - P(08+m) * [ funO(Y(2)) ; Y(in(m)+1:in(1+m)-1) ] ... - - P(08+m) * Y(in(m)+1:in(1+m) ) ; - - % update Generators(m) as last delay chamber value - Generators(m) = Y(in(1+m)) ; - end - end - - % Calculate rates of change for Q, C, Rc, Sc -------------------------- - - % Transfer Transfer/Decay - dY(1) = - P(05)*Y(1) ;% dQ/dt - dY(2) = P(06)*Y(1) - P(07)*Y(2) ;% dC/dt - dY(3) = P(09)*Generators(1) - P(09)*Y(3) ;% dRc/dt - dY(4) = P(10)*Generators(2) - P(10)*Y(4) ;% dSc/dt - - % Calculate rates of change for R, E, K, T, V, A ---------------------- - - % Proliferation Death - dY(5) = funF(Y,1)*P(18)/(P(18)+Y(5))*Y(5) - funG(Y,1) *Y(5) ;% dR/dt - dY(6) = - P(37) *Y(6) ;% dE/dt - dY(7) = funF(Y,2)*P(19)/(P(19)+Y(7))*Y(7) - funG(Y,2) *Y(7) ;% dK/dt - dY(8) = Prolif(3)*P(41)/(P(41)+Y(8))*Y(8) - P(40) *Y(8) ;% dT/dt - dY(9) = - P(40)*P(43)*Y(9) ;% dV/dt - dY(10:in(1)) = - P(34)*Y(10:in(1)) ;% dA/dt - - % Division - dY(10) = dY(10) - P(35)*Y(10) ;% dA1/dt - dY(11:in(1)) = dY(11:in(1)) + 2*P(35)*Y(10:in(1)-1) ... - - P(35)*Y(11:in(1)) ;% dAn/dt - dY(6) = dY(6) + 2*P(35)*Y( in(1)) ;% dE/dt - - % Activation/Deactivation - dY(5) = dY(5) - funH(Y,1)*Y(9) /(P(33)+Y(9))*Y(5) ;% dR/dt - dY(10) = dY(10) + funH(Y,1)*Y(9) /(P(33)+Y(9))*Y(5) ;% dA1/dt - dY(6) = dY(6) - P(38) *P(33)/(P(33)+Y(9))*Y(6) ;% dE/dt - dY(5) = dY(5) + P(38) *P(33)/(P(33)+Y(9))*Y(6) ;% dR/dt - - % Infection Killing - dY(8) = dY(8) - P(42)*Y(9)*Y(8) ;% dT/dt - dY(9) = dY(9) + beta *Y(9)*Y(8) - funH(Y,2)*Y(6)*Y(9) ... - - funH(Y,3)*Y(7)*Y(9) ;% dV/dt - -end - -% ------------------------------------------------------------------------- -% other functions --------------------------------------------------------- - -% function theta (equation for N803 effect) -function y = funO(C) - y = C / ( P(08) + C ) ; -end - -% function E (applies receptor regulation to N803 effect) -% s = 1 for R , s = 2 for E , s = 3 for K -function y = funE(Y,s) - y = funO(Y(2)) / ( 1 + P(12+s) * Y(3) ) ; -end - -% function F (R,K proliferation w/N803 modification) -% s = 1 for R , s = 2 for K -function y = funF(Y,s) - y = Prolif(s) + P(19+s) * funE(Y,s+(s==2)) ; -end - -% function G (R,K death w/immune suppression) -% s = 1 for R , s = 2 for K -function y = funG(Y,s) - y = P(15+s) + P(21+s) * Y(4) ; -end - -% function H (R activation, E,K killing w/immune suppression) -% s = 1 for R , s = 2 for E , s = 3 for K -function y = funH(Y,s) - y = ( P(23+s) + P(26+s) * funE(Y,s) ) / ( 1 + P(29+s) * Y(4) ) ; -end - -% ------------------------------------------------------------------------- -% Prepare main outputs 'Y_MAIN' ------------------------------------------- - -Y = abs(Y) ;% removing negatives resulting from numerical error - -Y_MAIN(:,1) = log10(Y(:,9))+3 ;% virion [log(#/mL)] -Y_MAIN(:,2) = sum(Y(:,[5,6,10:in(1)]),2) ;% CD8+ T cell (R+A+E) [/無] -Y_MAIN(:,3) = Y(:,7) ;% NK cell [/無] -Y_MAIN(:,4) = Y(:,8) ;% CD4+ T cell [/無] - -% ------------------------------------------------------------------------- -% Prepare extra outputs 'Y_ALL' ----------------------------------------- - -% If only one output is requested ('Y_MAIN'), exit function -if nargout == 1 ; return ;end - -% percolate outpute variable and declare associated line type, legend -% entry, and y-axis label -Y_ALL{01} = {Y(:,1), '-' , 'V' , 'log(mL^{-1})' }; - -Y_ALL{02} = {Y(:,1), '-' , 'H_{E}E' , 'frequency' ;... - Y(:,1), ':' , '1' , 'frequency' ;... - Y(:,1), ':' , '0' , 'frequency' }; - -Y_ALL{03} = {Y(:,1), '-' , 'H_{E}E' , 'log fold change' ;... - Y(:,1), '--', 'E' , 'log fold change' ;... - Y(:,1), ':' , 'H_{E}' , 'log fold change' }; - -Y_ALL{04} = {Y(:,1), '-' , 'H_{K}K' , 'log fold change' ;... - Y(:,1), '--', 'K' , 'log fold change' ;... - Y(:,1), ':' , 'H_{K}' , 'log fold change' }; - -Y_ALL{05} = {Y(:,1), '-' , 'T' , 'fold change' }; - -Y_ALL{06} = {Y(:,1), '- ', 'R' , 'fold change' }; - -Y_ALL{07} = {Y(:,1), '-' , 'F_{R}R^{*}-G_{R}' , 'd^{-1}';... - Y(:,1), '--', 'p_{R}R^{*}-d_{R}' , 'd^{-1}';... - Y(:,1), ':' , '\rho_{R}\xi_{R}R^{*}-\delta_{R}[IS]','d^{-1}'}; - -Y_ALL{08} = {Y(:,1), '-' , 'H_{R}' , 'log fold change';... - Y(:,1), '--', 'a_{R}+\alpha_{R}\xi_{R}', 'log fold change';... - Y(:,1), ':' , '1/(1+\sigma_{R}[IS])' , 'log fold change'}; - -Y_ALL{09} = {Y(:,1), '-' , '\xi_{R}' , ' ' ;... - Y(:,1), '--', '1/(1+\varsigma_{R}[RR])', ' ' ;... - Y(:,1), ':' , 'C/(C+\eta)' , ' ' }; - -Y_ALL{10} = {Y(:,1), '-' , 'C/(C+\eta)' , ' ' ;... - Y(:,1), '--', '[RR]' , ' ' ;... - Y(:,1), ':' , '[IS]' , ' ' }; - -Y_ALL{11} = {Y(:,1), '- ', 'A_{0}' , 'fold change' ;... - Y(:,1), '--', 'E' , 'fold change' }; - -Y_ALL{12} = {Y(:,1), '-' , 'A+E' , 'frequency' ;... - Y(:,1), '--', 'E' , 'frequency' ;... - Y(:,1), ':' , '1' , ' ' ;... - Y(:,1), ':' , '0' , ' ' }; - -Y_ALL{13} = {Y(:,1), '-' , 'H_{E}' , 'log fold change';... - Y(:,1), '--', 'a_{E}+\alpha_{E}\xi_{E}', 'log fold change';... - Y(:,1), ':' , '1/(1+\sigma_{E}[IS])' , 'log fold change'}; - -Y_ALL{14} = {Y(:,1), '-' , '\xi_{E}' , ' ' ;... - Y(:,1), '--', '1/(1+\varsigma_{E}[RR])', ' ' ;... - Y(:,1), ':' , 'C/(C+\eta)' , ' ' }; - -Y_ALL{15} = {Y(:,1), '-' , 'C/(C+\eta)' , ' ' ;... - Y(:,1), '--', '[RR]' , ' ' ;... - Y(:,1), ':' , '[IS]' , ' ' }; - -Y_ALL{16} = {Y(:,1), '- ', 'K' , 'fold change' }; - -Y_ALL{17} = {Y(:,1), '-' , 'F_{K}K^{*}-G_{K}' , 'd^{-1}';... - Y(:,1), '--', 'p_{K}K^{*}-d_{K}' , 'd^{-1}';... - Y(:,1), ':' , '\rho_{K}\xi_{K}K^{*}-\delta_{K}[IS]','d^{-1}'}; - -Y_ALL{18} = {Y(:,1), '-' , 'H_{K}' , 'log fold change';... - Y(:,1), '--', 'a_{K}+\alpha_{K}\xi_{K}', 'log fold change';... - Y(:,1), ':' , '1/(1+\sigma_{K}[IS])' , 'log fold change'}; - -Y_ALL{19} = {Y(:,1), '-' , '\xi_{K}' , ' ' ;... - Y(:,1), '--', '1/(1+\varsigma_{K}[RR])', ' ' ;... - Y(:,1), ':' , 'C/(C+\eta)' , ' ' }; - -Y_ALL{20} = {Y(:,1), '-' , 'C/(C+\eta)' , ' ' ;... - Y(:,1), '--', '[RR]' , ' ' ;... - Y(:,1), ':' , '[IS]' , ' ' }; - -% calculate extra output variables -for n = 1:length(SolTimes) - Yn = Y(n,:) ; - Killing = [ funH(Yn,2) * Y(n,6) , funH(Yn,3) * Y(n,7) ] ; - - Y_ALL{01}{1,1}(n) = log10(Y(n,9))+3 ; - - Y_ALL{02}{1,1}(n) = Killing(1) / sum(Killing) ; - Y_ALL{02}{2,1}(n) = 1 ; - Y_ALL{02}{3,1}(n) = 0 ; - - Y_ALL{03}{1,1}(n) = log10( Killing(1) / Y(1,6) / P(25) ) ; - Y_ALL{03}{2,1}(n) = log10( Y(n,6) / Y(1,6) ) ; - Y_ALL{03}{3,1}(n) = log10( funH(Yn,2) / P(25) ) ; - - Y_ALL{04}{1,1}(n) = log10( Killing(2) / Y(1,7) / P(26) ) ; - Y_ALL{04}{2,1}(n) = log10( Y(n,7) / Y(1,7) ) ; - Y_ALL{04}{3,1}(n) = log10( funH(Yn,3) / P(26) ) ; - - Y_ALL{05}{1,1}(n) = Y(n,8) / Y(1,8) ; - - Y_ALL{06}{1,1}(n) = Y(n,5) / Y(1,5) ; - - Y_ALL{07}{1,1}(n) = funF(Yn,1) * P(18)/(P(18)+Yn(5)) - funG(Yn,1) ; - Y_ALL{07}{2,1}(n) = Prolif(1) * P(18)/(P(18)+Yn(5)) - P(16) ; - Y_ALL{07}{3,1}(n) = Y_ALL{07}{1,1}(n) - Y_ALL{07}{2,1}(n) ; - - Y_ALL{08}{1,1}(n) = log10( funH(Yn,1) / P(24) ) ; - Y_ALL{08}{2,1}(n) = log10( 1 + P(27)*funE(Yn,1) / P(24) ) ; - Y_ALL{08}{3,1}(n) = log10( 1 / (1 + P(30)*Yn(4)) ) ; - - Y_ALL{09}{1,1}(n) = funE(Yn,1) ; - Y_ALL{09}{2,1}(n) = 1 / (1 + P(13)*Yn(3)) ; - Y_ALL{09}{3,1}(n) = funO(Yn(2)) ; - - Y_ALL{10}{1,1}(n) = funO(Yn(2)) ; - Y_ALL{10}{2,1}(n) = Yn(3) ; - Y_ALL{10}{3,1}(n) = Yn(4) ; - - Y_ALL{11}{1,1}(n) = Y(n,10) / Y(1,10) ; - Y_ALL{11}{2,1}(n) = Y(n,6) / Y(1,6) ; - - Y_ALL{12}{1,1}(n) = 1 - Y(n,5) / sum(Yn([5,6,10:in(1)])) ; - Y_ALL{12}{2,1}(n) = Y(n,6) / sum(Yn([5,6,10:in(1)])) ; - Y_ALL{12}{3,1}(n) = 1 ; - Y_ALL{12}{4,1}(n) = 0 ; - - Y_ALL{13}{1,1}(n) = log10( funH(Yn,2) / P(25) ) ; - Y_ALL{13}{2,1}(n) = log10( 1 + P(28)*funE(Yn,2) / P(25) ) ; - Y_ALL{13}{3,1}(n) = log10( 1 / (1 + P(31)*Yn(4)) ) ; - - Y_ALL{14}{1,1}(n) = funE(Yn,2) ; - Y_ALL{14}{2,1}(n) = 1 / (1 + P(14)*Yn(3)) ; - Y_ALL{14}{3,1}(n) = funO(Yn(2)) ; - - Y_ALL{15}{1,1}(n) = funO(Yn(2)) ; - Y_ALL{15}{2,1}(n) = Yn(3) ; - Y_ALL{15}{3,1}(n) = Yn(4) ; - - Y_ALL{16}{1,1}(n) = Y(n,7) / Y(1,7) ; - - Y_ALL{17}{1,1}(n) = funF(Yn,2) * P(19)/(P(19)+Yn(7)) - funG(Yn,2) ; - Y_ALL{17}{2,1}(n) = Prolif(2) * P(19)/(P(19)+Yn(7)) - P(17) ; - Y_ALL{17}{3,1}(n) = Y_ALL{17}{1,1}(n) - Y_ALL{17}{2,1}(n) ; - - Y_ALL{18}{1,1}(n) = log10( funH(Yn,3) / P(26) ) ; - Y_ALL{18}{2,1}(n) = log10( 1 + P(29)*funE(Yn,3) / P(26) ) ; - Y_ALL{18}{3,1}(n) = log10( 1 / (1 + P(32)*Yn(4)) ) ; - - Y_ALL{19}{1,1}(n) = funE(Yn,3) ; - Y_ALL{19}{2,1}(n) = 1 / (1 + P(15)*Yn(3)) ; - Y_ALL{19}{3,1}(n) = funO(Yn(2)) ; - - Y_ALL{20}{1,1}(n) = funO(Yn(2)) ; - Y_ALL{20}{2,1}(n) = Yn(3) ; - Y_ALL{20}{3,1}(n) = Yn(4) ; -end - -end