From 63ca8ca5534059df7c556dce55c8c065d587148b Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Sun, 25 Dec 2022 08:59:46 -0600 Subject: [PATCH] Delete HIV_model_1.m --- Legacy Code/HIV_model_1.m | 194 -------------------------------------- 1 file changed, 194 deletions(-) delete mode 100644 Legacy Code/HIV_model_1.m diff --git a/Legacy Code/HIV_model_1.m b/Legacy Code/HIV_model_1.m deleted file mode 100644 index 9961a7b..0000000 --- a/Legacy Code/HIV_model_1.m +++ /dev/null @@ -1,194 +0,0 @@ -%% -% /--------------------------------------------------------------\ -% | Date: 10/28/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 -% 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 inhibition pool (unique for E and K) - -%% ======================================================================== -% INPUTS -% ======================================================================== - -% t = vector of time points at which to evaluate solution [day] -% td = vector of time points at which to administer doses [day] -% P = vector of parameters corresponding to list below -% Q = vector of indices in P (used for equating parameter values) -% S = vector of 0/1 corresponding to P (not used in this version) - -% P(01) = V i.c. [log(/mL)] -% P(02) = E i.c. [/µL] -% P(03) = K i.c. [/µL] -% P(04) = N i.c. [pmol/kg] -% P(05) = C (vol. of dist.)/F [mL/kg] -% P(06) = C absorption rate [/d] -% P(07) = C elimination rate [/d] - -% P(08) = E decay rate [/d] -% P(09) = K decay rate [/d] -% P(10) = E killing rate [µL/d] -% P(11) = K killing rate [µL/d] -% P(12) = E Ce half saturation [nM] -% P(13) = K Ce half saturation [nM] - -% P(14) = E terminal fold change [] -% P(15) = K terminal fold change [] -% P(16) = E Ce killing modifier [] -% P(17) = K Ce killing modifier [] -% P(18) = E i.c. portion due to IL15 [] -% P(19) = K i.c. portion due to IL15 [] - -% P(20) = E Ct max level [] -% P(21) = K Ct max level [] -% P(22) = E Ct decay rate [/d] -% P(23) = K Ct decay rate [/d] -% P(24) = E Cy max level [] -% P(25) = K Cy max level [] -% P(26) = E Cy decay rate [/d] -% P(27) = K Cy decay rate [/d] - -%% ======================================================================== -% OUTPUTS -% ======================================================================== - -% YM(:,1) = V concentration at each time point in 't' [log(#/mL)] -% YM(:,2) = E fold change over baseline at each time point in 't' [] -% YM(:,3) = K fold change over baseline at each time point in 't' [] - -% Y = all variables (see 'system' function) -% YE = extra variables (see 'output' section) - -%% ======================================================================== -% FUNCTION -% ======================================================================== - -function [YM,Y,YE] = HIV_model_1(t,td,P,Q,S) %#ok - -P = P' ;% convert P to column vector - -if ~isempty(Q) - P = P(Q) ;% equate parameters based on indices in Q -% S = S(Q) ;% equate parameter switches based on indices in Q -end - -% Declare calculated parameters & initial conditions ---------------------- - -r = (1-(1-P(18:19))./P(14:15)).*(1+P(20:21)) ;% Ce proliferation modifier -O = P(18:19)./(r-P(18:19).*P(20:21)) ;% Ce initial effect - -if any(O<0)||any(O>1) - error('Initial effect is out of bounds.') -end - -P(01) = 10^(P(01)-3) ;% convert V to [#/µL] -I = [P(01:03) ; 0 ; 0] ;% V,E,K,N,C i.c. -I(6:9) = [P(20:21).*O ; P(24:25).*O] ;% Ct,Cy i.c. - -h = (P(12:13)/P(01)).*O./(1-O) ;% IL15 factor -s = P(08:09).*(1-P(18:19)).*P(02:03) ;% E,K generation rate -p = sum(P(10:11).*G(O,I(6:7),I(8:9)).*P(02:03)) ;% V production rate - -% Solve for each treatment cycle ------------------------------------------ - -if any(td==inf) - D = [1;find(ismember(t,0)) ;length(t)] ;% index in 't' for dose -else - D = [1;find(ismember(t,td));length(t)] ;% index in 't' for dose -end -Y = zeros(length(t),9) ;% percolate output matrix -I = I' ;% convert I to row vector - -for n = 1:length(D)-1 - - ts = tic ;% starting time (for catching a stalled algorithm execution) - - % Solve ODEs (see MATLAB help for ode15s algorithm) - [~,YTEMP] = ode15s(@system,t(D(n):D(n+1)),I) ; - - % For solutions of 2 timepoints, exclude the transition timepoints that - % ode15s adds by default - if length(t(D(n):D(n+1))) == 2 - YTEMP = [YTEMP(1,:);YTEMP(end,:)] ; - end - - Y(D(n):D(n+1),:) = YTEMP ;% add to full solution - I = Y(D(n+1),:) + [0 0 0 P(04) 0 0 0 0 0] ;% i.c. for next cycle -end - -% 'system' function ------------------------------------------------------- - -function dY = system(~,Y) - - if toc(ts) >= 5 - error('ODE solver stalled.') - end - - dY = zeros(9,1) ;% percolate dY/dt - - O = (h*Y(1)+Y(5))./(P(12:13)+h*Y(1)+Y(5)) ;% IL15+N803 effect - c = sum(P(10:11).*G(O,Y(6:7),Y(8:9)).*Y(2:3)) ;% E+K killing rate - d = P(08:09).*(1-r.*F(O,Y(6:7))) ;% E/K decay-prolif rate - - % Gen/Transfer Prolif/Decay Killing - dY(1) = + p*Y(1) - c*Y(1) ;% dV/dt - dY(2:3) = s - d.*Y(2:3) ;% dE,K/dt - dY(4) = - P(06)*Y(4) ;% dN/dt - dY(5) = P(06)*Y(4)/P(05) - P(07)*Y(5) ;% dC/dt - dY(6:7) = P(22:23).*P(20:21).*O - P(22:23).*Y(6:7) ;% dCt/dt - dY(8:9) = P(26:27).*P(24:25).*O - P(26:27).*Y(8:9) ;% dCt/dt - - if any(td==inf) - dY(4) = 0 ;% constant injection site quantity - end -end - -% other functions --------------------------------------------------------- - -% function y = io(in,A,B) -% if S(in) ; y = A ;% value if P(in) is on -% else ; y = B ;% value if P(in) is off -% end -% end - -function y = F(eff,tol) % Ce effect * tolerance - y = eff./(1+tol) ; -end - -function y = G(eff,tol,inh) % Ce killing factor - y = (1+P(16:17).*F(eff,tol))./(1+inh) ; -end - -% preparing outputs ------------------------------------------------------- - -Y = abs(Y) ;% removing negatives resulting from numerical error - -YM(:,1) = log10(Y(:,1))+3 ;% V [log(#/mL)] -YM(:,2) = Y(:,2)/P(02) ;% E fold change [] -YM(:,3) = Y(:,3)/P(03) ;% K fold change [] - -YE = zeros(length(t),11) ;% percolating extra ouput matrix -YE(:,1) = (h(1)*Y(:,1)+Y(:,5))./(P(12)+h(1)*Y(:,1)+Y(:,5)) ;% O for E -YE(:,2) = (h(2)*Y(:,1)+Y(:,5))./(P(13)+h(2)*Y(:,1)+Y(:,5)) ;% O for K -YE(:,3) = 1-1./(1+Y(:,6)) ;% tolerance for E -YE(:,4) = 1-1./(1+Y(:,7)) ;% tolerance for K -YE(:,5) = 1-1./(1+Y(:,8)) ;% inhibition for E -YE(:,6) = 1-1./(1+Y(:,9)) ;% inhibition for K -for n = 1:size(YE,1) % killing rate fold change - YE(n,7:8) = (G(YE(n,1:2)',Y(n,6:7)',Y(n,8:9)')./... - G(YE(1,1:2)',Y(1,6:7)',Y(1,8:9)'))' ; -end - -end