diff --git a/Legacy Code/HIV_treat_V3.m b/Legacy Code/HIV_treat_V3.m deleted file mode 100644 index 45f8533..0000000 --- a/Legacy Code/HIV_treat_V3.m +++ /dev/null @@ -1,245 +0,0 @@ -% -% /-------------------------------------------------------\ -% | Date: 6/24/2018 | -% | Author: Jonathan Cody | -% | Institution: Purdue University | -% | Weldon School of Biomedical Engineering | -% | Pienaar Lab | -% \-------------------------------------------------------/ -% -% Solves host-level treatment model V.3 for ALT-803 project -% -% Nomenclature: T = target CD4+ T cell -% I = infected CD4+ T cell -% R = resting CD8+ T cell -% E = effector CD8+ T cell -% K = natural killer cell -% V = virion -% -% INPUTS ================================================================== -% -% t = time points at which to evaluate solution [day] (inc. initial time) -% td = dosing regimen [day] -% P = parameter vector corresponding to the following list: -% -% R(:,1) = [X X] ;% T+I initial concentration [#/uL] -% R(:,2) = [X X] ;% R+E initial concentration [#/uL] -% R(:,3) = [X X] ;% K initial concentration [#/uL] -% R(:,4) = [X X] ;% V initial concentration [log(#/mL)] -% R(:,5) = [X X] ;% T decay rate constant [/day] -% R(:,6) = [X X] ;% I decay rate constant [*P(5)] -% R(:,7) = [X X] ;% R decay rate constant [/day] -% R(:,8) = [X X] ;% E decay rate constant [*P(7)] -% R(:,9) = [X X] ;% K decay rate constant [/day] -% R(:,10) = [X X] ;% V decay rate constant [/day] -% R(:,11) = [X X] ;% ALT803 initial concentration [nM] -% R(:,12) = [X X] ;% I/T initial concentration proportionality [] -% R(:,13) = [X X] ;% E/R initial concentration proportionality [] -% R(:,14) = [X X] ;% I killing by all E initial rate [/day] -% R(:,15) = [X X] ;% I killing by all K initial rate [*P(14)] -% R(:,16) = [X X] ;% E prolif. by IL15 sat. rate constant [*P(8*7)] -% R(:,17) = [X X] ;% K prolif. by IL15 sat. rate constant [*P(9)] -% R(:,18) = [X X] ;% E prolif. by IL15 half-sat. conc. [#/uL] -% R(:,19) = [X X] ;% K prolif. by IL15 half-sat. conc. [#/uL] -% R(:,20) = [X X] ;% R prolif. by ALT803 sat. rate constant [/day] -% R(:,21) = [X X] ;% E prolif. by ALT803 sat. rate constant [/day] -% R(:,22) = [X X] ;% K prolif. by ALT803 sat. rate constant [/day] -% R(:,23) = [X X] ;% ALT803 decay rate constant [/day] -% R(:,24) = [X X] ;% ALT803 half-saturation concentration [nM] -% R(:,25) = [X X] ;% tolerance production sat. rate constant [/day] -% R(:,26) = [X X] ;% tolerance decay rate constant [/day] -% R(:,27) = [X X] ;% tolerance half-saturation concentration [] -% R(:,28) = [X X] ;% effect production sat. rate constant [/day] -% R(:,29) = [X X] ;% effect decay rate constant [/day] -% R(:,30) = [X X] ;% effect half-saturation concentration [] -% R(:,31) = [X X] ;% induced killing by E sat. rate constant [*cE] -% R(:,32) = [X X] ;% induced killing by K sat. rate constant [*cK] -% -% OUTPUTS ================================================================= -% -% Ym(:,1) = V concentration at each time point in 't' [log(#/mL)] -% Ym(:,2) = T+I concentration at each time point in 't' [#/uL] -% Ym(:,3) = R+E concentration at each time point in 't' [#/uL] -% Ym(:,4) = K concentration at each time point in 't' [#/uL] -% -% Ya = all concentrations and rate terms (see 'preparing outputs') -% -% Pa = all parameter values (see 'preparing outputs') -% -% CALCULATIONS ============================================================ - -function [Ym,Ya,Pa] = HIV_treat_V3(t,td,P) - -% scaling to prevent negative concentrations due to rounding error -------- - -P([1:3,11,18,19,24]) = P([1:3,11,18,19,24])*1e6 ;% [#/L] or [fM] -P(4) = 10.^P(4)*1e3 ;% [#/L] - -% generating constants and ICs by assuming initially steady state --------- - -T0 = P(1)/(1+P(12)) ;% T initial condition -R0 = P(2)/(1+P(13)) ;% R initial condition -I0 = T0*P(12) ;% I initial condition -E0 = R0*P(13) ;% E initial condition -dI = P(6)*P(5) ;% I decay rate constant -dE = P(8)*P(7) ;% E decay rate constant -pE = P(16)*dE ;% E <- IL15 sat. rate constant -pK = P(17)*P(9) ;% K <- IL15 sat. rate constant -cE = P(14)/E0 ;% I killing by E rate constant -cK = P(15)*P(14)/P(3) ;% I killing by K rate constant -b = (dI+P(14)*(1+P(15)))*P(12)/P(4) ;% T infection rate constant -a = (dE/P(4)-pE/(P(18)+P(4)))*P(13) ;% R activation rate constant -sT = (P(5)+b*P(4))*T0 ;% T generation rate -sR = (P(7)+a*P(4))*R0 ;% R generation rate -sK = (P(9)-pK*P(4)/(P(19)+P(4)))*P(3) ;% K generation rate -pV = P(10)*P(4)/I0 ;% V production by I rate constant - -% solving for each treatment cycle ---------------------------------------- - -I = [T0 I0 R0 E0 P(3:4) 0 0 0] ;% initial conditions -d = [1 find(ismember(t,td)) length(t)] ;% index in 't' for dose -Y = zeros(length(t),9) ;% percolating output matrix - -for n = 1:length(d)-1 - % Solving ODEs (see MATLAB documentation for algorithm) - [~,Ytemp] = ode15s(@system,t(d(n):d(n+1)),I) ; - - if length(t(d(n):d(n+1))) == 2 % resolving issue with solver output - Ytemp = [Ytemp(1,:);Ytemp(end,:)] ; - end - if size(Ytemp,1)~=length(t(d(n):d(n+1))) % if solver fails - Ym = 1 ; return - end - Y(d(n):d(n+1),:) = Ytemp ;% building full solution - I = [Y(d(n+1),1:6) P(11) Y(d(n+1),8:9)] ;% I.C. for next cycle -end - -% 'system' function ------------------------------------------------------- - -function dY = system(~,Y) -dY = zeros(9,1) ;% percolating dY/dt -% ALT803 effect -Ef = Y(9)/(P(30)+Y(9))*P(27)/(P(27)+Y(8)) ; -% T rate of change -dY(1) = sT - P(5)*Y(1) - b*Y(6)*Y(1) ; -% I rate of change -dY(2) = b*Y(6)*Y(1) - dI*Y(2) - cE*(1+Ef*P(31))*Y(4)*Y(2) ... - - cK*(1+Ef*P(32))*Y(5)*Y(2) ; -% R rate of change -dY(3) = sR - P(7)*Y(3) - a*Y(6)*Y(3) + P(20)*Ef*Y(3) ; -% E rate of change -dY(4) = a*Y(6)*Y(3) - dE*Y(4) + pE*Y(6)/(P(18)+Y(6))*Y(4) + P(21)*Ef*Y(4) ; -% K rate of change -dY(5) = sK - P(9)*Y(5) + pK*Y(6)/(P(19)+Y(6))*Y(5) + P(22)*Ef*Y(5) ; -% V rate of change -dY(6) = pV*Y(2) - P(10)*Y(6) ; -% ALT803 concentration rate of change -dY(7) = -P(23)*Y(7) ; -% ALT803 tolerance concentration rate of change -dY(8) = P(25)*Y(7)/(P(24)+Y(7))-P(26)*Y(8) ; -% ALT803 effect concentration rate of change -dY(9) = P(28)*Y(7)/(P(24)+Y(7))-P(29)*Y(9) ; -end - -% preparing outputs ------------------------------------------------------- - -Y = abs(Y) ;% required to prevent imaginary numbers - -Ym = [log10(Y(:,6)/1e3) [Y(:,1)+Y(:,2) Y(:,3)+Y(:,4) Y(:,5)]/1e6] ; - -Ef = Y(:,9)./(P(30)+Y(:,9))*P(27)./(P(27)+Y(:,8)) ; - -Ya{1} = {... -'[T]' Y(:,1)/1e6 '\muL^{-1}' ;... -'s_{T}' sT/1e6*(1+0*Y(:,1)) '\muL^{-1} d^{-1}' ;... -'d_{T}[T]' P(5)*Y(:,1)/1e6 '\muL^{-1} d^{-1}' ;... -'b[V][T]' b*Y(:,6).*Y(:,1)/1e6 '\muL^{-1} d^{-1}'}; -Ya{2} = {... -'[I]' Y(:,2)/1e6 '\muL^{-1}' ;... -'b[V][T]' b*Y(:,6).*Y(:,1)/1e6 '\muL^{-1} d^{-1}' ;... -'d_{I}[I]' dI*Y(:,2)/1e6 '\muL^{-1} d^{-1}' ;... -'c_{E}[E][I]' cE*Y(:,4).*Y(:,2)/1e6 '\muL^{-1} d^{-1}' ;... -'c_{K}[K][I]' cK*Y(:,5).*Y(:,2)/1e6 '\muL^{-1} d^{-1}' ;... -'\kappa_{E}[ef][E][I]' cE*P(31)*Ef.*Y(:,4)... - .*Y(:,2)/1e6 '\muL^{-1} d^{-1}' ;... -'\kappa_{K}[ef][K][I]' cK*P(32)*Ef.*Y(:,5)... - .*Y(:,2)/1e6 '\muL^{-1} d^{-1}'}; -Ya{3} = {... -'[R]' Y(:,3)/1e6 '\muL^{-1}' ;... -'s_{R}' sR/1e6*(1+0*Y(:,3)) '\muL^{-1} d^{-1}' ;... -'d_{R}[R]' P(7)*Y(:,3)/1e6 '\muL^{-1} d^{-1}' ;... -'a[V][R]' a*Y(:,6).*Y(:,3)/1e6 '\muL^{-1} d^{-1}' ;... -'\rho_{R}[ef][R]' P(20)*Ef.*Y(:,3)/1e6 '\muL^{-1} d^{-1}'}; -Ya{4} = {... -'[E]' Y(:,4)/1e6 '\muL^{-1}' ;... -'a[V][R]' a*Y(:,6).*Y(:,3)/1e6 '\muL^{-1} d^{-1}' ;... -'d_{E}[E]' dE*Y(:,4)/1e6 '\muL^{-1} d^{-1}' ;... -'p_{E}[E][V]/(h_{E}+[V])' pE*Y(:,4).*Y(:,6)/1e6... - ./(P(18)+Y(:,6)) '\muL^{-1} d^{-1}' ;... -'\rho_{E}[ef][E]' P(21)*Ef.*Y(:,4)/1e6 '\muL^{-1} d^{-1}'}; -Ya{5} = {... -'[K]' Y(:,5)/1e6 '\muL^{-1}' ;... -'s_{K}' sK/1e6*(1+0*Y(:,5)) '\muL^{-1} d^{-1}' ;... -'d_{K}[K]' P(9)*Y(:,5)/1e6 '\muL^{-1} d^{-1}' ;... -'p_{K}[K][V]/(h_{K}+[V])' pK*Y(:,5).*Y(:,6)/1e6... - ./(P(19)+Y(:,6)) '\muL^{-1} d^{-1}' ;... -'\rho_{K}[ef][K]' P(22)*Ef.*Y(:,5)/1e6 '\muL^{-1} d^{-1}'}; -Ya{6} = {... -'[V]' Y(:,6)/1e6 '\muL^{-1}' ;... -'p_{V}[I]' pV*Y(:,2)/1e6 '\muL^{-1} d^{-1}' ;... -'d_{V}[V]' P(10)*Y(:,6)/1e6 '\muL^{-1} d^{-1}'}; -Ya{7} = {... -'Effect' Ef ' ' ;... -'[C]' Y(:,7)/1e6 'nM' ;... -'\delta[C]' P(23)*Y(:,7)/1e6 'nM d^{-1}' }; -Ya{8} = {... -'[C_{/tau}]' Y(:,8) ' ' ;... -'\rho_{\tau}[C]/(\eta+[C])' P(25)*Y(:,7)... - ./(P(24)+Y(:,7)) 'd^{-1}' ;... -'\delta_{\tau}[C_{/tau}]' P(26)*Y(:,8) 'd^{-1}' }; -Ya{9} = {... -'[C_{/epsilon}]' Y(:,9) ' ' ;... -'\rho_{\epsilon}[C]/(\eta+[C])' P(28)*Y(:,7)... - ./(P(24)+Y(:,7)) 'd^{-1}' ;... -'\delta_{\epsilon}[C_{/epsilon}]' P(29)*Y(:,9) 'd^{-1}' }; - -Pa = {'T_{0}' T0/1e6 '\muL^{-1}' ;... - 'I_{0}' I0/1e6 '\muL^{-1}' ;... - 'R_{0}' R0/1e6 '\muL^{-1}' ;... - 'E_{0}' E0/1e6 '\muL^{-1}' ;... - 'K_{0}' P(3)/1e6 '\muL^{-1}' ;... - 'V_{0}' P(4)/1e6 '\muL^{-1}' ;... - 's_{T}' sT/1e6 '\muL^{-1} d^{-1}' ;... - 's_{R}' sR/1e6 '\muL^{-1} d^{-1}' ;... - 's_{K}' sK/1e6 '\muL^{-1} d^{-1}' ;... - 'd_{T}' P(5) 'd^{-1}' ;... - 'd_{I}' dI 'd^{-1}' ;... - 'd_{R}' P(7) 'd^{-1}' ;... - 'd_{E}' dE 'd^{-1}' ;... - 'd_{K}' P(9) 'd^{-1}' ;... - 'd_{V}' P(10) 'd^{-1}' ;... - 'p_{V}' pV 'd^{-1}' ;... - 'b' b*1e6 '\muL d^{-1}' ;... - 'a' a*1e6 '\muL d^{-1}' ;... - 'c_{E}' cE*1e6 '\muL d^{-1}' ;... - 'c_{K}' cK*1e6 '\muL d^{-1}' ;... - 'p_{E}' pE 'd^{-1}' ;... - 'p_{K}' pK 'd^{-1}' ;... - 'h_{E}' P(18)/1e6 '\muL^{-1}' ;... - 'h_{K}' P(19)/1e6 '\muL^{-1}' ;... - '\rho_{R}' P(20) 'd^{-1}' ;... - '\rho_{E}' P(21) 'd^{-1}' ;... - '\rho_{K}' P(22) 'd^{-1}' ;... - 'C_{0}' P(11)/1e6 'nM' ;... - '\delta' P(23) 'd^{-1}' ;... - '\eta' P(24)/1e6 'nM' ;... - '\rho_{\tau}' P(25) 'd^{-1}' ;... - '\delta_{\tau}' P(26) 'd^{-1}' ;... - '\eta_{\tau}' P(27) ' ' ;... - '\rho_{\epsilon}' P(28) 'd^{-1}' ;... - '\delta_{\epsilon}' P(29) 'd^{-1}' ;... - '\eta_{\epsilon}' P(30) ' ' ;... - '\kappa_{E}' cE*P(31)*1e6 '\muL d^{-1}' ;... - '\kappa_{K}' cK*P(32)*1e6 '\muL d^{-1}' }; - -end