From 7ab9e8396213dc80d39474c85bf4fa5a75f4754f Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Sun, 25 Dec 2022 08:59:51 -0600 Subject: [PATCH] Delete HIV_model_2.m --- Legacy Code/HIV_model_2.m | 388 -------------------------------------- 1 file changed, 388 deletions(-) delete mode 100644 Legacy Code/HIV_model_2.m diff --git a/Legacy Code/HIV_model_2.m b/Legacy Code/HIV_model_2.m deleted file mode 100644 index 88b7264..0000000 --- a/Legacy Code/HIV_model_2.m +++ /dev/null @@ -1,388 +0,0 @@ -%% -% /--------------------------------------------------------------\ -% | Date: 01/23/2018 | -% | Author: Jonathan Cody | -% | Affiliation: Purdue University | -% | Weldon School of Biomedical Engineering | -% | Pienaar Computational Systems Pharmacology Lab | -% \--------------------------------------------------------------/ - -% Solves host-level model for virus and effector cell concentrations during -% a subcutaneously administered regimen of N803 - -% Nomenclature: V = virion -% T = CD4+ T cell -% E = CD8+ T cell -% K = NK cell -% N = N803 injection site quantity -% C = N803 plasma concentration -% Ce = IL15+N803 effective pool -% Ct = IL15+N803 tolerance pool (unique for E and K) -% Cy = IL15+N803 exhaustion pool (unique for E and K) - -%% ======================================================================== -% INPUTS -% ======================================================================== - -% SolTimes = ascending vector of days at which to evaluate solution -% DoseTimes = ascending vector of days at which to administer doses -% MutanTimes = ascending vector of days at which new mutations emerge -% -> NOTE: all elements of 'DoseTimes' and 'MutanTimes' must also be -% elements of 'SolTimes' - -% P = vector of parameters corresponding to list below -% -> NOTE: parameter pairs are for E and K, respectively - -% P(01) = V i.c. [log(/mL)] -% P(02) = E i.c. [/無] -% P(03) = K i.c. [/無] - -% P(04) = N i.c. [pmol/kg] -% P(05) = C absorption rate const. [/d] -% P(06) = P(05) * 'bioavailability'/'volume of distribution' [kg/L-d] -% P(07) = C elimination rate const. [/d] - -% P(08) = death rate constant [/d] -% P(09) -% P(10) = maximum concentration [/無] -% P(11) - -% P(12) = killing rate constant [無/d] -% P(13) -% P(14) = Ce killing modifier [] -% P(15) - -% P(16) = Ce half saturation [pM] -% P(17) -% P(18) = constant IL15 factor [pM] -% P(19) -% P(20) = V-dep IL15 factor [pM*無] -% P(21) - -% P(22) = Ct max level [] -% P(23) -% P(24) = Ct delay rate constant [/d] -% P(25) - -% P(26) = Cy max level [] -% P(27) -% P(28) = Cy delay rate constant [/d] -% P(29) - -% P(30) = escape mutation rate [] -% P(31) = initial escape proportion [] -% P(32) = response activation threshold [/無] -% P(33) = response prolif. rate const. [/d] - -% P(34) = T i.c. [/無] -% P(35) = T death rate constant [/d] -% P(36) = T maximum concentration [/無] -% P(37) = T infection rate constant [/無-d] -% P(38) = infected T decay rate constant [*(P(35)] - -% EquatedPars = vector of indices in 'P' used for equating parameter values -% (if empty, EquatedPars is ignored) -% timeout = time out a stalled solution after this many seconds -% (if empty, timeout = 10 [seconds]) - -%% ======================================================================== -% OUTPUTS -% ======================================================================== - -% Y_MAIN = model output at each time point in 'SolTimes' -% Y_MAIN(:,1) = V [log(#/mL)] -% Y_MAIN(:,2) = E [/無] -% Y_MAIN(:,3) = K [/無] -% Y_MAIN(:,4) = T [/無] - -% Y_EXTRA = cell array of extra outputs (see 'output' section) - -%% ======================================================================== -% FUNCTION -% ======================================================================== - -function [Y_MAIN,Y_EXTRA] = HIV_model_2... - (SolTimes,DoseTimes,MutanTimes,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(MutanTimes,1)==1 ; MutanTimes = MutanTimes' ; 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 vectors to store strain mutation rates and response proliferation -% rates (both initally 0) and declare number of mutant strains -MutanRates = MutanTimes*0 ; -RespoRates = MutanTimes*0 ; -mutanum = length(MutanTimes) ; - -% build reference index in 'SolTimes' containing start time, all dose and -% mutation times, and end time (used for piecewise solution of ODEs) -PiecewiseIndex = [ 1 ; ... - find(ismember(SolTimes,unique([DoseTimes;MutanTimes]))) ; ... - length(SolTimes) ] ; - -% ------------------------------------------------------------------------- -% Declare initial conditions & calculated parameters ---------------------- - -V = 10^(P(01)-3) ;% converting V i.c. to [/無] - -V = V *10^(3*mutanum) ;% scaling up V to avoid rounding errors -P(20:21) = P(20:21)/10^(3*mutanum) ;% scaling parameters to match -P(32) = P(32) *10^(3*mutanum) ; -P(37) = P(37) /10^(3*mutanum) ; - -Yinit = [ 0 ; 0 ... N:C i.c. (pre-treatment phase) - ; P(22:23).*funO(V,0) ... Ct for E:K i.c. - ; P(26:27).*funO(V,0) ... Cy for E:K i.c. - ; P([02:03,34]) ... E:K,T i.c. - ; V ... V i.c. - ; MutanTimes*0 ... V mutant i.c. - ; MutanTimes*0 + P(31) ];% V mutant escape proportion i.c. - -% E:K maximum proliferation rate constant [/d] (from initial steady-state) -maxprolifEK = P(08:09)./ funF(V,0,Yinit(3:4))./ (1-Yinit(7:8)./P(10:11)) ; - -% T maximum proliferation rate constant [/d] (from initial steady-state) -maxprolifT = (P(35)+P(37)*V)/ (1-Yinit(9)/P(36)) ; - -% V 'production' by T rate constant [/d] (from initial steady-state) -betaV = sum(P(12:13).* funG(V,0,Yinit(3:4),Yinit(5:6)).* Yinit(7:8)) ; -betaV = (betaV + P(35)*P(38))/Yinit(9) ; - -% ------------------------------------------------------------------------- -% Solve for each solution period ------------------------------------------ - -Y = zeros(length(SolTimes),length(Yinit)) ;% percolate output matrix -Yinit = Yinit' ;% convert 'Yinit' to row vector - -for n = 1:length(PiecewiseIndex)-1 - - solvertime = tic ;% start timer (for catching a stalled execution) - - % update 'MutanRates' if a mutation start time is reached - MutanRates( MutanTimes <= SolTimes(PiecewiseIndex(n)) ) = P(30) ; - - % if a dose time is reached, add dose to N803 injection site quantity - if ismember(SolTimes(PiecewiseIndex(n)),DoseTimes) - Yinit = Yinit + [P(04) Yinit(2:end)*0 ] ; - end - - % solve ODEs (see MATLAB help for ode15s algorithm) - [~,Y_TEMP] = ode15s(@system,... - SolTimes(PiecewiseIndex(n):PiecewiseIndex(n+1)),... - Yinit) ; - - % 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 - Yinit = 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 - - V = sum(Y(10:10+mutanum)) ;% total virions - - % Calculate rates of change for N, C, Ct, Cy, E, K , T ---------------- - - % E:K modified proliferation rate [/d] - modprolifEK = maxprolifEK... - .* (1 - Y(7:8)./P(10:11)).* funF(V,Y(2),Y(3:4)) ; - - % T modified proliferation rate [/d] - modprolifT = maxprolifT * (1 - Y(9)/P(36)) ; - - % Transfer/Gen./Prolif. Transfer/Decay - dY(1) = - P(05)*Y(1) ;% dN/dt - dY(2) = P(06)*Y(1) - P(07)*Y(2) ;% dC/dt - dY(3:4) = P(24:25).*P(22:23).*funO(V,Y(2)) - P(24:25).*Y(3:4) ;% dCt/dt - dY(5:6) = P(28:29).*P(26:27).*funO(V,Y(2)) - P(28:29).*Y(5:6) ;% dCy/dt - dY(7:8) = modprolifEK.*Y(7:8) - P(08:09).*Y(7:8) ;%dE:K/dt - dY(9) = modprolifT.*Y(9) - P(35)*Y(9) ... - ... Infection - - P(37)*V*Y(9) ;% dT/dt - - % Calculate rates of change for V, mutants, and escape factors -------- - - % V killing by E:K rate constants [/d] - killingEK = P(12:13).* funG(V,Y(2),Y(3:4),Y(5:6)).* Y(7:8) ; - - % dV/dt:dVm1/dt:dVm2/dt:etc - dY(10:10+mutanum) = ... - ... Proliferation - betaV*Y(9)*(1-[MutanRates;0]-[0;MutanRates]).*Y(10:10+mutanum) ... - ... Mutation - + betaV*Y(9)*[0;MutanRates].*[0;Y(10:9+mutanum)] ... - + betaV*Y(9)*[MutanRates;0].*[Y(11:10+mutanum);0] ... - ... Decay - - P(35)*P(38)*Y(10:10+mutanum) ... - ... Killing - - ( killingEK(1)*(1-[0;Y(11+mutanum:10+2*mutanum)]) ... - + killingEK(2) ).*Y(10:10+mutanum) ; - - % update 'RespoRates' if mutant exceeds response activation threshold - RespoRates( Y(11:10+mutanum) >= P(32) ) = P(33) ; - - % decay of mutant escape proportions (i.e. proliferation of response) - dY(11+mutanum:10+2*mutanum) = -RespoRates.*Y(11+mutanum:10+2*mutanum) ; -end - -% ------------------------------------------------------------------------- -% other functions --------------------------------------------------------- - -% function theta (equation for IL15+N803 effect) -function y = funO(V,C) - Ce = P(18:19)+P(20:21)*V+C ;% Ce for E:K - y = Ce./(P(16:17)+Ce) ;% Ce effect for E:K -end - -% function F (applies tolerance for IL15+N803 effect) -function y = funF(V,C,Ct) - y = funO(V,C)./(1+Ct) ;% Ce effect * tolerance for E:K -end - -% function G (IL15+N803 killing rate modifier, including exhaustion) -function y = funG(V,C,Ct,Cy) - y = (1+P(14:15).*funF(V,C,Ct))./(1+Cy) ;% Ce killing factor for E:K -end - -% ------------------------------------------------------------------------- -% Prepare main outputs 'Y_MAIN' ------------------------------------------- - -Y(:,10:10+mutanum) = Y(:,10:10+mutanum)/10^(3*mutanum) ;% scaling back V -P(20:21) = P(20:21) *10^(3*mutanum) ;% scaling to match -P(32) = P(32) /10^(3*mutanum) ; -P(37) = P(37) *10^(3*mutanum) ; - -Y = abs(Y) ;% removing negatives resulting from numerical error - -Y_MAIN(:,1) = log10(sum(Y(:,10:10+mutanum),2))+3 ;% V [log(#/mL)] -Y_MAIN(:,2:4) = Y(:,7:9) ;% E:K:T [/無] - -% ------------------------------------------------------------------------- -% Prepare extra outputs 'Y_EXTRA' ----------------------------------------- - -% percolate outpute variable and declare associated line type, legend -% entry, and y-axis label -Y_EXTRA{01} = {SolTimes,'-' ,'V' ,'log(mL^{-1})'}; -Y_EXTRA{02} = {SolTimes,'--' ,'c_{E}HG_{E}E' ,'frequency' ;... - SolTimes,':' ,'c_{K}G_{K}K' ,'frequency' }; -Y_EXTRA{03} = {SolTimes,'-' ,'E' ,'fold change' ;... - SolTimes,'--' ,'HG_{E}' ,'fold change' }; -Y_EXTRA{04} = {SolTimes,'-' ,'K' ,'fold change' ;... - SolTimes,'--' ,'G_{K}' ,'fold change' }; -Y_EXTRA{05} = {SolTimes,'-' ,'H' ,' ' }; -Y_EXTRA{06} = {SolTimes,'-' ,'E' ,'fold change' }; -Y_EXTRA{07} = {SolTimes,':' ,'E_{m}-E' ,'fold change' ;... - SolTimes,'--' ,'F_{E}' ,'fold change' ;... - SolTimes,'-' ,'product' ,'fold change' }; -Y_EXTRA{08} = {SolTimes,'-' ,'G_{E}' ,'fold change' ;... - SolTimes,'--' ,'(1+\varsigma_{E}F_{E})' ,'fold change' ;... - SolTimes,':' ,'(1+C_{\gamma,E})^{-1}' ,'fold change' }; -Y_EXTRA{09} = {SolTimes,'-' ,'F_{E}' ,'fold change' ;... - SolTimes,'--' ,'\theta_{E}' ,'fold change' ;... - SolTimes,':' ,'(1+C_{\tau,E})^{-1}' ,'fold change' }; -Y_EXTRA{10} = {SolTimes,'--' ,'\Lambda_{E}' ,'effect' ;... - SolTimes,':' ,'+ \lambda_{E}V' ,'effect' ;... - SolTimes,'-' ,'+ C' ,'effect' }; -Y_EXTRA{11} = {SolTimes,'-' ,'K' ,'fold change' }; -Y_EXTRA{12} = {SolTimes,':' ,'K_{m}-K' ,'fold change' ;... - SolTimes,'--' ,'F_{K}' ,'fold change' ;... - SolTimes,'-' ,'product' ,'fold change' }; -Y_EXTRA{13} = {SolTimes,'-' ,'G_{K}' ,'fold change' ;... - SolTimes,'--' ,'(1+\varsigma_{K}F_{K})' ,'fold change' ;... - SolTimes,':' ,'(1+C_{\gamma,K})^{-1}' ,'fold change' }; -Y_EXTRA{14} = {SolTimes,'-' ,'F_{K}' ,'fold change' ;... - SolTimes,'--' ,'\theta_{K}' ,'fold change' ;... - SolTimes,':' ,'(1+C_{\tau,K})^{-1}' ,'fold change' }; -Y_EXTRA{15} = {SolTimes,'--' ,'\Lambda_{K}' ,'effect' ;... - SolTimes,':' ,'+ \lambda_{K}V' ,'effect' ;... - SolTimes,'-' ,'+ C' ,'effect' }; - -% calculate extra output variables -funO0 = funO(Y(1,9),Y(1,2)) ; -funF0 = funF(Y(1,9),Y(1,2),Y(1,3:4)) ; -funG0 = funG(Y(1,9),Y(1,2),Y(1,3:4),Y(1,5:6)) ; -for n = 1:length(SolTimes) - - V = sum(Y(n,10:10+mutanum)) ; - funOn = funO(V,Y(n,2)) ; - funFn = funF(V,Y(n,2),Y(n,3:4)) ; - funGn = funG(V,Y(n,2),Y(n,3:4),Y(n,5:6)) ; - H = sum((1-[0,Y(n,11+mutanum:10+2*mutanum)]).* Y(n,10:10+mutanum)/ V) ; - killingEK = P(12:13) .* funGn .* Y(n,7:8) ; - - Y_EXTRA{01}{1,1}(n) = log10(V)+3 ; - - Y_EXTRA{02}{1,1}(n) = H*killingEK(1)/ (H*killingEK(1) + killingEK(2)) ; - Y_EXTRA{02}{2,1}(n) = killingEK(2)/ (H*killingEK(1) + killingEK(2)) ; - - Y_EXTRA{03}{1,1}(n) = Y(n,7) / Y(1,7) ; - Y_EXTRA{03}{2,1}(n) = H * funGn(1) / funG0(1) ; - - Y_EXTRA{04}{1,1}(n) = Y(n,8) / Y(1,8) ; - Y_EXTRA{04}{2,1}(n) = funGn(2) / funG0(2) ; - - Y_EXTRA{05}{1,1}(n) = H ; - - Y_EXTRA{06}{1,1}(n) = Y(n,7) / Y(1,7) ; - - Y_EXTRA{07}{1,1}(n) = (P(10) - Y(n,7)) / (P(10) - Y(1,7)) ; - Y_EXTRA{07}{2,1}(n) = funFn(1) / funF0(1) ; - Y_EXTRA{07}{3,1}(n) = Y_EXTRA{07}{1,1}(n) * Y_EXTRA{07}{2,1}(n) ; - - Y_EXTRA{08}{1,1}(n) = funGn(1) / funG0(1) ; - Y_EXTRA{08}{2,1}(n) = (1 + P(14) * funFn(1)) / (1 + P(14) * funF0(1)) ; - Y_EXTRA{08}{3,1}(n) = (1 + Y(1,5)) / (1 + Y(n,5)) ; - - Y_EXTRA{09}{1,1}(n) = funFn(1) / funF0(1) ; - Y_EXTRA{09}{2,1}(n) = funOn(1) / funO0(1) ; - Y_EXTRA{09}{3,1}(n) = (1 + Y(1,3)) / (1 + Y(n,3)) ; - - Y_EXTRA{10}{1,1}(n) = P(18) / (P(16) + P(18)) ; - Y_EXTRA{10}{2,1}(n) = (P(18) + P(20)*V) / (P(16) + P(18) + P(20)*V) ; - Y_EXTRA{10}{3,1}(n) = funOn(1) ; - - Y_EXTRA{11}{1,1}(n) = Y(n,8) / Y(1,8) ; - - Y_EXTRA{12}{1,1}(n) = (P(11) - Y(n,8)) / (P(11) - Y(1,8)) ; - Y_EXTRA{12}{2,1}(n) = funFn(2) / funF0(2) ; - Y_EXTRA{12}{3,1}(n) = Y_EXTRA{12}{1,1}(n) * Y_EXTRA{12}{2,1}(n) ; - - Y_EXTRA{13}{1,1}(n) = funGn(2) / funG0(2) ; - Y_EXTRA{13}{2,1}(n) = (1 + P(15) * funFn(2)) / (1 + P(15) * funF0(2)) ; - Y_EXTRA{13}{3,1}(n) = (1 + Y(1,6)) / (1 + Y(n,6)) ; - - Y_EXTRA{14}{1,1}(n) = funFn(2) / funF0(2) ; - Y_EXTRA{14}{2,1}(n) = funOn(2) / funO0(2) ; - Y_EXTRA{14}{3,1}(n) = (1 + Y(1,4)) / (1 + Y(n,4)) ; - - Y_EXTRA{15}{1,1}(n) = P(19) / (P(17) + P(19)) ; - Y_EXTRA{15}{2,1}(n) = (P(19) + P(21)*V) / (P(17) + P(19) + P(21)*V) ; - Y_EXTRA{15}{3,1}(n) = funOn(2) ; -end - -end