diff --git a/Legacy Code/HIV_model_4.m b/Legacy Code/HIV_model_4.m deleted file mode 100644 index 744574d..0000000 --- a/Legacy Code/HIV_model_4.m +++ /dev/null @@ -1,418 +0,0 @@ -%% -% /--------------------------------------------------------------\ -% | Date: 04/16/2819 | -% | 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 = virions (sum of optional strains V1 and V2) -% T = CD4+ T cells -% E = CD8+ T cells -% K = NK cells -% Q = N803 injection site quantity -% C = N803 plasma concentration -% TOL = cellular tolerance (unique for E and K) -% REG = immune regulation (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 -% (elements of 'DoseTimes' must also be in 'SolTimes') - -% P = vector of parameters corresponding to list below -% (parameter pairs are for E and K, respectively) - -% EquatedPars = vector of indices in 'P' used for equating parameter values -% (for example, EquatedPars(09) = 08 will set P(09) = P(08) -% (if empty, EquatedPars is ignored) - -% timeout = scalar time limit [seconds] for solving the ODEs -% (if timeout is exceed, an error is thrown) -% (if empty, timeout = 10) - -% normout = scalar used to change units of Y_MAIN (see OUTPUTS) -% (if empty, normout is ignored) - -% P(01) = V i.c. [log(/mL)] (will convert to [/無]) -% 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) = total initial killing rate [/d] (will convert to [/無-d]) -% P(13) = [*P(12)] (will convert to [/無-d]) -% P(14) = drug killing modifier [] -% P(15) - -% P(16) = drug half saturation concentration [pM] -% P(17) -% P(18) = max drug-induced proliferation rate [/d] -% P(19) - -% P(20) = TOL delay rate constant [/d] -% P(21) -% P(22) = REG delay rate constant [/d] -% P(23) -% P(24) = # of TOL delay chambers [] -% P(25) -% P(26) = # of REG delay chambers [] -% P(27) - -% P(28) = TOL drug effect modifier [] -% P(29) -% P(30) = REG-induced decay rate [*P(18)] (will convert to [/d]) -% P(31) [*P(19)] (will convert to [/d]) -% P(32) = REG killing modifier [] -% P(33) - -% P(34) = T i.c. [/無] -% P(35) = T death rate constant [/d] -% P(36) = T max proliferating concentration [/無] -% P(37) = T infection by V1 rate constant [/無-d] -% P(38) = infected T death rate [*P(35)] (will convert to [/d]) - -% P(39) = V2 initial frequency [] -% P(40) = mutation probability [] -% P(41) = V2 targeting factor [] (if P(41) = 0 , P(41) -> 1) -% P(42) (if P(42) = 0 , P(42) -> 1) - -%% ======================================================================== -% OUTPUTS -% ======================================================================== - -% if normout = 1 -% Y_MAIN(:,1) = V at points in 'SolTimes' [log(#/mL)] [log fold change] -% Y_MAIN(:,2) = E at points in 'SolTimes' [/無] [fold change] -% Y_MAIN(:,3) = K at points in 'SolTimes' [/無] [fold change] -% Y_MAIN(:,4) = T at points in 'SolTimes' [/無] [fold change] -% (see INPUTS) - -% Y_ALL = cell array of extra outputs (see end of script) - -%% ======================================================================== -% FUNCTION -% ======================================================================== - -function [Y_MAIN,Y_ALL] = HIV_model_4(SolTimes,DoseTimes,P,... - EquatedPars,timeout,normout) - -% 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 ... -% [(E,K TOL delay chambers) (E,K REG delay chambers)] -In = 11 + [0 P(24) sum(P(24:25)) sum(P(24:26)) sum(P(24:27))] ; - -% build reference index in 'SolTimes' containing start time, all dose -% times, and end time (used for piecewise solution of ODEs) -PieceIn = [ 1 ; find(ismember(SolTimes,DoseTimes)) ... - ; length(SolTimes) ] ; - -% ------------------------------------------------------------------------- -% Declare initial conditions & calculated parameters ---------------------- - -P(01) = 10^(P(01)-3) ;% convert from [log(/mL)] to [/無] - -Yic = zeros(1,In(5)) ;% initial conditions -Yic(7:11) = [P([02,03,34])' P(01)*[1-P(39) P(39)] ] ;% E,K,T,V1,V2 i.c. - -P(13) = P(13)*P(12) ;% convert from [*P(12)] to [/d] -P(38) = P(38)*P(35) ;% convert from [*P(35)] to [/d] - -if P(41) == 0 ; P(41) = 1 ; end % if V2 targeting factor is zero (off)... -if P(42) == 0 ; P(42) = 1 ; end % then set to 1 (no reduction in targeting) - -% calculate beta [無/d] and b [無/d] -if P(39) == 0 % for single strain - beta = [ (sum(P(12:13)) + P(38)) / P(34) ; 0 ] ; - b = [ P(37) ; 0 ] ; -else % for multi-strain - betaCOEF = [ 1-P(40) P(40)*P(39)/(1-P(39)) ... - ; P(40)*(1-P(39))/P(39) 1-P(40) ]; - betacon = ([ sum( P(12:13)) ... - ; sum(P(41:42).* P(12:13)) ] + P(38) )./ P(34) ; - beta = betaCOEF \ betacon ; - b = P(37) * [ 1 ; beta(2)/beta(1) ] ; -end - -P(30:31) = P(30:31).*P(18:19) ;% convert from [*P(18:19)] to [/d] -P(12:13) = P(12:13)./P(02:03) ;% convert from [/d] to [/無-d] - -% if a maximum proliferating concentration is zero, set it as very large -% (this avoids division by zero and approximately removes its effect) -if P(10)==0 ; P(10) = P(02)*10^30 ; end -if P(11)==0 ; P(11) = P(03)*10^30 ; end -if P(36)==0 ; P(36) = P(34)*10^30 ; end - -% calculate proliferation rate constants [/d] -rEK = P(08:09) .* (P(10:11) + P(02:03))./P(10:11) ; -rT = (P(35)+Yic(10:11)*b) * (P(36) + P(34) ) /P(36) ; - -% ------------------------------------------------------------------------- -% Solve for each solution period ------------------------------------------ - -Y = zeros(length(SolTimes),length(Yic)) ;% percolate output matrix - -for n = 1:length(PieceIn)-1 - - solvertime = tic ;% start timer (for catching a stalled execution) - - % solve ODEs (see MATLAB help for ode15s algorithm) - [~,Y_TEMP] = ode15s(@system,SolTimes(PieceIn(n):PieceIn(n+1)),Yic) ; - - % for solutions of 2 timepoints, exclude the transition timepoints - % that ode15s adds by default - if length(SolTimes(PieceIn(n):PieceIn(n+1))) == 2 - Y_TEMP = [Y_TEMP(1,:);Y_TEMP(end,:)] ; - end - - % add 'Y_TEMP' to full solution 'Y' - Y(PieceIn(n):PieceIn(n+1),:) = Y_TEMP ; - - % update initial condition for next solution window - Yic = Y(PieceIn(n+1),:) + [ P(04) Yic(2:end)*0 ] ; -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 TOL and REG delay chambers ------------ - - Effects = [funO(Y(2)) ; funO(Y(2))] ;% drug effects for E,K,E,K - Generators = Effects ;% terms that generate tolerance and regulation - % (by default these are the drug effects) - - % for E tolerance, K tolerance, E regulation, K regulation - for m = 1:4 - - % if delay chambers are specified - if P(23+m) ~= 0 - - % update delay chamber rates of change - dY(In(m)+1:In(m+1)) = ... - P(19+m) * [ Effects(m) ; Y(In(m)+1:In(m+1)-1) ] ... - - P(19+m) * Y(In(m)+1:In(m+1) ) ; - - % update Generators(m) as last delay chamber value - Generators(m) = Y(In(m+1)) ; - end - end - - % Calculate rates of change for Q, C, TOL, REG, E, K, T, V ------------ - - % E,K modified proliferation rate [/d] - prolifEK = rStar(Y(2),Y(3:4)).* P(10:11)./ (P(10:11) + Y(7:8)) ; - - % T modified proliferation rate [/d] - prolifT = rT * P(36) / (P(36) + Y(9)) ; - - % V killing by E,K rate [/d] - Killing = cStar(Y(2),Y(3:4),Y(5:6)).* Y(7:8) ; - - % 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(20:21).*Generators(1:2) - P(20:21) .*Y(3:4) ;% dTOL/dt - dY(5:6) = P(22:23).*Generators(3:4) - P(22:23) .*Y(5:6) ;% dREG/dt - dY(7:8) = prolifEK.*Y(7:8) - dStar(Y(5:6)).*Y(7:8) ;% dE:K/dt - dY(9) = prolifT *Y(9) - P(35) *Y(9) ... - ... Infection/Mutation Decay/Killing - - sum( b.*Y(10:11) ) *Y(9) ;% dT/dt - dY(10) = beta(1)*Y(10)*(1-P(40))*Y(9) - P(38) *Y(10) ... - + beta(2)*Y(11)* P(40) *Y(9) - sum(Killing)*Y(10) ;% dV1/dt - - dY(11) = beta(2)*Y(11)*(1-P(40))*Y(9) - P(38) *Y(11) ... - + beta(1)*Y(10)* P(40) *Y(9) - sum(Killing ... - .*P(41:42))*Y(11) ;% dV2/dt -end - -% ------------------------------------------------------------------------- -% other functions --------------------------------------------------------- - -% function theta (equation for N803 effect) -function y = funO(C) - y = C./(P(16:17)+C) ; -end - -% function xi (applies tolerance to N803 effect) -function y = funE(C,TOL) - y = funO(C)./ (1 + P(28:29).*TOL) ; -end - -% function rStar (base proliferation + N803-induced proliferation) -function y = rStar(C,TOL) - y = rEK + P(18:19).* funE(C,TOL) ; -end - -% function dStar (base death + regulation-induced decay) -function y = dStar(REG) - y = P(08:09) + P(30:31).*REG ; -end - -% function cStar (base killing rate with N803 and regulation modifiers) -function y = cStar(C,TOL,REG) - y = P(12:13).* (1 + P(14:15).* funE(C,TOL))./ (1 + P(32:33).*REG) ; -end - -% ------------------------------------------------------------------------- -% Prepare main outputs 'Y_MAIN' ------------------------------------------- - -Y = abs(Y) ;% removing negatives resulting from numerical error - -if normout == 1 % units are [log fold change] and [fold change] - Y_MAIN(:,1) = log10( (Y(:,10)+Y(:,11)) / P(01) ) ;% V - Y_MAIN(:,2:4) = Y(:,7:9) ./ P([02,03,34])' ;% E,K,T -else % units are [log(#/mL)] and [/無] - Y_MAIN(:,1) = log10( Y(:,10)+Y(:,11) ) + 3 ;% V - Y_MAIN(:,2:4) = Y(:,7:9) ;% E,K,T -end - -% ------------------------------------------------------------------------- -% Prepare extra outputs 'Y_ALL' ------------------------------------------- - -% if only one output is requested ('Y_MAIN'), exit function -if nargout == 1 ; return ;end - -% generate placeholder for extra outputs -Z = SolTimes * 0 ; - -% percolate outpute variable and declare associated line type, legend -% entry, and y-axis label -Y_ALL{01} = {Z, '-' , 'V' , 'log fold \Delta' }; - -Y_ALL{02} = {Z, '-' , 'c_{K}^{*}K(\chi_{K})' , 'frequency' ;... - Z, ':' , '1' , 'frequency' ;... - Z, ':' , '0' , 'frequency' }; - -Y_ALL{03} = {Z, ':' , 'c_{E}^{*}E' , 'log fold \Delta' ;... - Z, '-' , 'c_{E}^{*}E(\chi_{E})' , 'log fold \Delta' }; - -Y_ALL{04} = {Z, ':' , 'c_{K}^{*}K' , 'log fold \Delta' ;... - Z, '-' , 'c_{K}^{*}K(\chi_{K})' , 'log fold \Delta' }; - -Y_ALL{05} = {Z, '-' , 'V2' , 'frequency' ;... - Z, ':' , '1' , 'frequency' ;... - Z, ':' , '0' , 'frequency' }; - -Y_ALL{06} = {Z, '-' , 'E' , 'fold \Delta' }; - -Y_ALL{07} = {Z, ':' , 'r_{E}(E_{p})-d_{E}' , 'd^{-1}' ;... - Z, '-' , 'r_{E}^{*}(E_{p})-d_{E}^{*}' , 'd^{-1}' }; - -Y_ALL{08} = {Z, ':' , 'c_{E}^{*}(-REG)' , 'log fold \Delta' ;... - Z, '-' , 'c_{E}^{*}' , 'log fold \Delta' }; - -Y_ALL{09} = {Z, ':' , '\theta_{E}' , ' ' ;... - Z, '-' , '\xi_{E}' , ' ' }; - -Y_ALL{10} = {Z, ':' , 'TOL_{E}' , ' ' ;... - Z, '-' , 'REG_{E}' , ' ' }; - -Y_ALL{11} = {Z, '-' , 'K' , 'fold \Delta' }; - -Y_ALL{12} = {Z, ':' , 'r_{K}(K_{p})-d_{K}' , 'd^{-1}' ;... - Z, '-' , 'r_{K}^{*}(K_{p})-d_{K}^{*}' , 'd^{-1}' }; - -Y_ALL{13} = {Z, ':' , 'c_{K}^{*}(-REG)' , 'log fold \Delta' ;... - Z, '-' , 'c_{K}^{*}' , 'log fold \Delta' }; - -Y_ALL{14} = {Z, ':' , '\theta_{K}' , ' ' ;... - Z, '-' , '\xi_{K}' , ' ' }; - -Y_ALL{15} = {Z, ':' , 'TOL_{K}' , ' ' ;... - Z, '-' , 'REG_{K}' , ' ' }; - -% calculate extra output variables -for n = 1:length(Z) - - Funs(1:2) = funO (Y(n,2) ) ;% theta at time n - Funs(3:4) = funE (Y(n,2),Y(n,3:4)' ) ;% xi at time n - Funs(5:6) = rStar(Y(n,2),Y(n,3:4)' ) ;% r* at time n - Funs(7:8) = dStar( Y(n,5:6)') ;% d* at time n - Funs(9:10) = cStar(Y(n,2),Y(n,3:4)',Y(n,5:6)') ;% c* at time n - Funs(11:12) = ( Y(n,10) + P(41:42)*Y(n,11) ) ... - / ( Y(n,10) + Y(n,11) ) ;% escape at time n - Funs(13:14) = Funs(11:12).* Funs(9:10).* Y(n,7:8) ;% killing at time n - - if n == 1 ; Funs0 = Funs ; end % saving initial values - - Y_ALL{01}{1,1}(n) = log10( (Y(n,10)+Y(n,11)) / P(01) ) ; - - Y_ALL{02}{1,1}(n) = Funs(14) / sum(Funs(13:14)) ; - Y_ALL{02}{2,1}(n) = 1 ; - Y_ALL{02}{3,1}(n) = 0 ; - - Y_ALL{03}{1,1}(n) = log10( Funs(9)*Y(n,7) / Funs0(9)/Y(1,7) ) ; - Y_ALL{03}{2,1}(n) = log10( Funs(13) / Funs0(13) ) ; - - Y_ALL{04}{1,1}(n) = log10( Funs(10)*Y(n,8) / Funs0(10)/Y(1,8) ) ; - Y_ALL{04}{2,1}(n) = log10( Funs(14) / Funs0(14) ) ; - - Y_ALL{05}{1,1}(n) = Y(n,11) / (Y(n,10)+Y(n,11)) ; - Y_ALL{05}{2,1}(n) = 1 ; - Y_ALL{05}{3,1}(n) = 0 ; - - Y_ALL{06}{1,1}(n) = Y(n,7) / Y(1,7) ; - - Y_ALL{07}{1,1}(n) = rEK(1) * P(10)/(P(10)+Y(n,7)) - P(08) ; - Y_ALL{07}{2,1}(n) = Funs(5) * P(10)/(P(10)+Y(n,7)) - Funs(7) ; - - Y_ALL{08}{1,1}(n) = log10( 1 + P(14) * Funs(3) ) ; - Y_ALL{08}{2,1}(n) = log10( Funs(9) / P(12) ) ; - - Y_ALL{09}{1,1}(n) = Funs(1) ; - Y_ALL{09}{2,1}(n) = Funs(3) ; - - Y_ALL{10}{1,1}(n) = Y(n,3) ; - Y_ALL{10}{2,1}(n) = Y(n,5) ; - - Y_ALL{11}{1,1}(n) = Y(n,8) / Y(1,8) ; - - Y_ALL{12}{1,1}(n) = rEK(2) * P(11)/(P(11)+Y(n,8)) - P(09) ; - Y_ALL{12}{2,1}(n) = Funs(6) * P(11)/(P(11)+Y(n,8)) - Funs(8) ; - - Y_ALL{13}{1,1}(n) = log10( 1 + P(15) * Funs(4) ) ; - Y_ALL{13}{2,1}(n) = log10( Funs(10) / P(13) ) ; - - Y_ALL{14}{1,1}(n) = Funs(2) ; - Y_ALL{14}{2,1}(n) = Funs(4) ; - - Y_ALL{15}{1,1}(n) = Y(n,4) ; - Y_ALL{15}{2,1}(n) = Y(n,6) ; -end - -end