diff --git a/Legacy Code/HIV_model_3.m b/Legacy Code/HIV_model_3.m deleted file mode 100644 index 4c275ac..0000000 --- a/Legacy Code/HIV_model_3.m +++ /dev/null @@ -1,421 +0,0 @@ -%% -% /--------------------------------------------------------------\ -% | Date: 02/19/2019 | -% | 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 -% Q = N803 injection site quantity -% C = N803 plasma concentration -% Ct = N803 tolerance (unique for E and K) -% Cy = N803 exhaustion (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) = Q i.c. [pmol/kg] -% P(05) = C absorption rate constant [/d] -% P(06) = P(05) * 'bioavailability'/'volume of distribution' [kg/L-d] -% P(07) = C elimination rate constant [/d] - -% P(08) = death rate constant [/d] -% P(09) -% P(10) = max proliferating concentration [/無] -% P(11) - -% P(12) = killing rate constant [無/d] -% P(13) -% P(14) = drug killing modifier -% P(15) - -% P(16) = drug half saturation [pM] -% P(17) -% P(18) = drug max prolif rate [/d] -% P(19) - -% P(20) = tolerance modifier [] -% P(21) -% P(22) = Ct delay rate constant [/d] -% P(23) - -% P(24) = exhaustion modifier [] -% P(25) -% P(26) = Cy delay rate constant [/d] -% P(27) - -% P(28) = escape mutation rate [] -% P(29) = initial escape proportion [] -% P(30) = response activation threshold [/無] -% P(31) = response proliferation rate constant [/d] - -% P(32) = T i.c. [/無] -% P(33) = T death rate constant [/d] -% P(34) = T max proliferating concentration [/無] -% P(35) = T infection rate constant [/無-d] -% P(36) = infected T decay rate constant [*(P(33)] - -% Delays = number of delay chambers for: -% [E tolerance, K tolerance, E exhaustion, K exhaustion] -% (if empty, Delays = [0 0 0 0]) -% 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_3... - (SolTimes,DoseTimes,MutanTimes,... - P,Delays,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 - -% set default 'Delays' to [0 0 0 0] -if isempty(Delays) ; Delays = [0 0 0 0] ; 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) ; - -% declare last index in Y (solution) for: -mutend = 10 + mutanum ;% mutant strains -resend = mutend + mutanum ;% responses to mutant strains -DelEnd = [1 1 1 1 1]*resend ;% responses to mutant strains and ... -for n = 2:5 - DelEnd(n:5) = DelEnd(n:5) + Delays(n-1) ;% delay chambers -end - -% store delays rates for E:K tolerance and exhaustion -DelRates = [P(22) P(23) P(26) P(27)] ; - -% 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(30) = P(30) *10^(3*mutanum) ; -P(35) = P(35) /10^(3*mutanum) ; - -Yinit = [ 0 ; 0 ... Q:C i.c. (pre-treatment phase) - ; 0 ; 0 ; 0 ; 0 ... Ct,Cy for E:K i.c. - ; P([02:03,32]) ... E:K,T i.c. - ; V ... V i.c. - ; MutanTimes*0 ... V mutant i.c. - ; MutanTimes*0 + P(29) ... V mutant escape proportion i.c. - ; zeros(sum(Delays),1) ];% E:K tol:exh delay chambers i.c. - -% V 'production' by T rate constant [/d] (from initial steady-state) -betaV = (sum(P(12:13).* Yinit(7:8)) + P(33)*P(36)) /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(28) ; - - % 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 - - % Calculate rates of change for Ct and Cy delay chambers -------------- - - Effects = funO(Y(2)) ;% drug effects for E:K - Effects = [Effects(1) Effects(2) Effects(1) Effects(2)] ; - Generators = Effects' ;% terms that generate tolerance and exhaustion - % (by default these are the drug effects) - - % for E tolerance, K tolerance, E exhaustion, K exhaustion - for m = 1:4 - - % if delay chambers are specified - if Delays(m) ~= 0 - - % update delay chamber rates of change - dY(DelEnd(m)+1:DelEnd(m+1)) = ... - DelRates(m)*[Effects(m);Y(DelEnd(m)+1:DelEnd(m+1)-1)] ... - - DelRates(m)*Y(DelEnd(m)+1:DelEnd(m+1)) ; - - % update Generator(m) as last delay chamber value - Generators(m) = Y(DelEnd(m+1)) ; - end - end - - % Calculate rates of change for Q, C, Ct, Cy, E, K , T ---------------- - - % E:K modified proliferation rate [/d] - prolifEK = P(08:09).*(P(10:11)+P(02:03))./(P(10:11) + Y(7:8)) ... - + P(18:19).*funF(Y(2),Y(3:4)) ; - - % T modified proliferation rate [/d] - prolifT = (P(33)+P(35)*V)*(P(34)+P(32))/(P(34) + Y(9)) ; - - % Transfer/Prolif. Transfer/Decay - dY(1) = - P(05)*Y(1) ;% dQ/dt - dY(2) = P(06)*Y(1) - P(07)*Y(2) ;% dC/dt - dY(3:4) = P(22:23).*Generators(1:2) - P(22:23).*Y(3:4) ;% dCt/dt - dY(5:6) = P(26:27).*Generators(3:4) - P(26:27).*Y(5:6) ;% dCy/dt - dY(7:8) = prolifEK.*Y(7:8) - P(08:09).*Y(7:8) ;%dE:K/dt - dY(9) = prolifT.*Y(9) - P(33)*Y(9) ... - ... Infection - - P(35)*sum(Y(10:mutend))*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(Y(2),Y(3:4),Y(5:6)).* Y(7:8) ; - - % dV/dt:dVm1/dt:dVm2/dt:etc - dY(10:mutend) = ... - ... Proliferation - betaV*Y(9)*(1-[MutanRates;0]-[0;MutanRates]).*Y(10:mutend) ... - ... Mutation - + betaV*Y(9)*[0;MutanRates].*[0;Y(10:mutend-1)] ... - + betaV*Y(9)*[MutanRates;0].*[Y(11:mutend);0] ... - ... Decay - - P(33)*P(36)*Y(10:mutend) ... - ... Killing - - ( killingEK(1)*(1-[0;Y(mutend+1:resend)]) ... - + killingEK(2) ).*Y(10:mutend) ; - - % update 'RespoRates' if mutant exceeds response activation threshold - RespoRates( Y(11:mutend) >= P(30) ) = P(31) ; - - % decay of mutant escape proportions (i.e. proliferation of response) - dY(mutend+1:resend) = -RespoRates.*Y(mutend+1:resend) ; - -end - -% ------------------------------------------------------------------------- -% other functions --------------------------------------------------------- - -% function theta (equation for N803 effect) -function y = funO(C) - y = C./(P(16:17)+C) ; -end - -% function F (applies tolerance to N803 effect) -function y = funF(C,Ct) - y = funO(C)-Ct.*P(20:21) ; -end - -% function G (N803 killing rate modifier, including exhaustion) -function y = funG(C,Ct,Cy) - f = funO(C)./(1+Ct.*P(20:21)) ; - y = (1+P(14:15).*f)./(1+Cy.*P(24:25)) ; -end - -% ------------------------------------------------------------------------- -% Prepare main outputs 'Y_MAIN' ------------------------------------------- - -Y(:,10:10+mutanum) = Y(:,10:10+mutanum)/10^(3*mutanum) ;% scaling back V -P(30) = P(30) /10^(3*mutanum) ; -P(35) = P(35) *10^(3*mutanum) ; - -Y = abs(Y) ;% removing negatives resulting from numerical error - -Y_MAIN(:,1) = log10(sum(Y(:,10:mutend),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' ,'log fold change' ;... - SolTimes,':' ,'G_{E}H' ,'log fold change' ;... - SolTimes,'-' ,'EG_{E}H' ,'log fold change' }; -Y_EXTRA{04} = {SolTimes,'--' ,'K' ,'log fold change' ;... - SolTimes,':' ,'G_{K}' ,'log fold change' ;... - SolTimes,'-' ,'KG_{K}' ,'log fold change' }; -Y_EXTRA{05} = {SolTimes,'--' ,'G_{E}' ,'log fold change' ;... - SolTimes,':' ,'H' ,'log fold change' ;... - SolTimes,'-' ,'G_{E}H' ,'log fold change' }; -Y_EXTRA{06} = {SolTimes,'-' ,'E' ,'fold change' }; -Y_EXTRA{07} = {SolTimes,'--' ,'r_{E}*E_{m}/(E_{m}+E)-d_{E}' ,' ' ;... - SolTimes,':' ,'\rho_{E}*F_{E}' ,' ' ;... - SolTimes,'-' ,'sum' ,' ' }; -Y_EXTRA{08} = {SolTimes,'--' ,'\theta_{E}' ,' ' ;... - SolTimes,':' ,'\sigma_{\tau,E}C_{\tau,E}' ,' ' ;... - SolTimes,'-' ,'F_{E}' ,' ' }; -Y_EXTRA{09} = {SolTimes,'--' ,'G_{E} (no exh)' ,'log' ;... - SolTimes,':' ,'1/(1+\sigma_{\gamma,E}*C_{\gamma,E})',' ';... - SolTimes,'-' ,'G_{E}' ,'log' }; -Y_EXTRA{10} = {SolTimes,'-' ,'\theta_{E}' ,' ' ;... - SolTimes,'--' ,'C_{\tau,E}' ,' ' ;... - SolTimes,':' ,'C_{\gamma,E}' ,' ' }; -Y_EXTRA{11} = {SolTimes,'-' ,'K' ,'fold change' }; -Y_EXTRA{12} = {SolTimes,'--' ,'r_{K}*K_{m}/(K_{m}+K)-d_{K}' ,' ' ;... - SolTimes,':' ,'\rho_{K}*F_{K}' ,' ' ;... - SolTimes,'-' ,'sum' ,' ' }; -Y_EXTRA{13} = {SolTimes,'--' ,'\theta_{K}' ,' ' ;... - SolTimes,':' ,'\sigma_{\tau,K}C_{\tau,K}' ,' ' ;... - SolTimes,'-' ,'F_{K}' ,' ' }; -Y_EXTRA{14} = {SolTimes,'--' ,'G_{K} (no exh)' ,'log' ;... - SolTimes,':' ,'1/(1+\sigma_{\gamma,K}*C_{\gamma,K})',' ';... - SolTimes,'-' ,'G_{K}' ,'log' }; -Y_EXTRA{15} = {SolTimes,'-' ,'\theta_{K}' ,' ' ;... - SolTimes,'--' ,'C_{\tau,K}' ,' ' ;... - SolTimes,':' ,'C_{\gamma,K}' ,' ' }; - -% calculate extra output variables -for n = 1:length(SolTimes) - - V = sum(Y(n,10:mutend)) ; - funOn = funO(Y(n,2)) ; - funFn = funF(Y(n,2),Y(n,3:4)) ; - funGn = funG(Y(n,2),Y(n,3:4),Y(n,5:6)) ; - H = sum((1-[0,Y(n,mutend+1:resend)]).* Y(n,10:mutend)/ 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) = log10( Y(n,7) / Y(1,7) ) ; - Y_EXTRA{03}{2,1}(n) = log10( H * funGn(1) ) ; - Y_EXTRA{03}{3,1}(n) = log10( H * funGn(1) * Y(n,7) / Y(1,7) ) ; - - Y_EXTRA{04}{1,1}(n) = log10( Y(n,8) / Y(1,8) ) ; - Y_EXTRA{04}{2,1}(n) = log10( funGn(2) ) ; - Y_EXTRA{04}{3,1}(n) = log10( funGn(2) * Y(n,8) / Y(1,8) ) ; - - Y_EXTRA{05}{1,1}(n) = log10( funGn(1) ) ; - Y_EXTRA{05}{2,1}(n) = log10( H ) ; - Y_EXTRA{05}{3,1}(n) = log10( H * funGn(1) ) ; - - Y_EXTRA{06}{1,1}(n) = Y(n,7) / Y(1,7) ; - - Y_EXTRA{07}{1,1}(n) = P(8) * ( (P(10)+P(02))/(P(10)+Y(n,7)) - 1 ); - Y_EXTRA{07}{2,1}(n) = P(18) * funFn(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) = funOn(1) ; - Y_EXTRA{08}{2,1}(n) = P(20)*Y(n,3) ; - Y_EXTRA{08}{3,1}(n) = funFn(1) ; - - Y_EXTRA{09}{1,1}(n) = log10( 1 + P(14)*funOn(1) / (1+P(20)*Y(n,3)) ) ; - Y_EXTRA{09}{2,1}(n) = log10( 1 / (1 + P(24)*Y(n,5)) ) ; - Y_EXTRA{09}{3,1}(n) = log10( funGn(1) ) ; - - Y_EXTRA{10}{1,1}(n) = funOn(1) ; - Y_EXTRA{10}{2,1}(n) = Y(n,3) ; - Y_EXTRA{10}{3,1}(n) = Y(n,5) ; - - Y_EXTRA{11}{1,1}(n) = Y(n,8) / Y(1,8) ; - - Y_EXTRA{12}{1,1}(n) = P(9) * ( (P(11)+P(03))/(P(11)+Y(n,8)) - 1 ); - Y_EXTRA{12}{2,1}(n) = P(19) * funFn(1) ; - Y_EXTRA{12}{3,1}(n) = Y_EXTRA{12}{1,1}(n) + Y_EXTRA{12}{2,1}(n) ; - - Y_EXTRA{13}{1,1}(n) = funOn(2) ; - Y_EXTRA{13}{2,1}(n) = P(21)*Y(n,4) ; - Y_EXTRA{13}{3,1}(n) = funFn(2) ; - - Y_EXTRA{14}{1,1}(n) = log10( 1 + P(15)*funOn(2) / (1+P(21)*Y(n,4)) ) ; - Y_EXTRA{14}{2,1}(n) = log10( 1 / (1 + P(25)*Y(n,6)) ) ; - Y_EXTRA{14}{3,1}(n) = log10( funGn(2) ) ; - - Y_EXTRA{15}{1,1}(n) = funOn(2) ; - Y_EXTRA{15}{2,1}(n) = Y(n,4) ; - Y_EXTRA{15}{3,1}(n) = Y(n,6) ; -end - -end