From f2e491b24f91eaea4a9c40067466f7411950c834 Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Thu, 28 Oct 2021 20:15:02 -0400 Subject: [PATCH] Add files via upload --- N803_model_2.m | 284 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 284 insertions(+) create mode 100644 N803_model_2.m diff --git a/N803_model_2.m b/N803_model_2.m new file mode 100644 index 0000000..50252c8 --- /dev/null +++ b/N803_model_2.m @@ -0,0 +1,284 @@ +% N803_model_2.m - solves host-level model for SIV virions, CD8+ T cells, +% and NK cells during a subcutaneously administered regimen of N803 + +% /--------------------------------------------------------------\ +% | Date: 07/12/2021 | +% | Author: Jonathan Cody | +% | Affiliation: Pienaar Computational Systems Pharmacology Lab | +% | Weldon School of Biomedical Engineering | +% | Purdue University | +% \--------------------------------------------------------------/ +% +% Version 2A +% +% Nomenclature: V = SIV virions [#/無] +% K = NK cells [#/無] +% E = CD8+ T cells [#/無] (E0 resting, EA active) +% X = N803 at absorption site [pmol/kg] +% C = N803 plasma concentration [pM] +% R = regulation [] (dimensionless quantity) +% T = tolerance [] (dimensionless quantity) + +%% ======================================================================== +% INPUTS +% ======================================================================== + +% SoluTimes = 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 'SoluTimes') + +% P = vector of parameters (see list in function) + +% EquatedPars = vector of indices in 'P' used for equating parameter values +% (for example, EquatedPars(18) = 17 will set P(18) = P(17) +% (if empty, EquatedPars is ignored) + +% timeOut = scalar time limit [seconds] for ODE solver 'ode15s' +% (if exceeded during 'ode15s' execution, an error is thrown) +% (if empty, timeOut = 10) + +% fullOut = switch to determine output variables (see end of code) +% (if empty, fullOut is ignored) + +%% ======================================================================== +% OUTPUTS +% ======================================================================== + +% Y_OUT(:,1) = V at points in 'SoluTimes' [log fold change] +% Y_OUT(:,2) = E at points in 'SoluTimes' [fold change] +% Y_OUT(:,3) = K at points in 'SoluTimes' [fold change] + +% If 'fullOut' is non-empty, more columns are included (see end of code) + +%% ======================================================================== +% FUNCTION +% ======================================================================== + +function Y_OUT = N803_model_2(SoluTimes,DoseTimes,Pars,EquatedPars,fullOut) + +% convert inputs to column vectors +if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end +if size(DoseTimes ,1)==1 ; DoseTimes = DoseTimes' ; end +if size(Pars ,1)==1 ; Pars = Pars' ; end + +% declare index in 'SoluTimes' for start, dose, and end times +Pieces = [ 1 ; find(ismember(SoluTimes,DoseTimes)) ; length(SoluTimes) ] ; + +% equate parameters in 'P' based on indices in 'EquatedPars' +if ~isempty(EquatedPars) ; Pars = Pars(EquatedPars) ; end + +%% ------------------------------------------------------------------------ +% Rename inputed parameters ----------------------------------------------- + +Vi = Pars(01) ;% V initial condition [log(#/mL)] +Ki = Pars(02) ;% K initial condition [#/無] +Ei = Pars(03) ;% E initial condition [#/無] +fEA = Pars(04) ;% active E initial frequency + +gK_Ki = Pars(05)*Pars(06) ;% K initial killing rate [/d] +gE_EAi = Pars(06) ;% E initial killing rate [/d] +Km = Pars(07) ;% K max proliferating concentration [#/無] +Em = Pars(08) ;% E max proliferating concentration [#/無] +V50 = Pars(09) ;% V half-sat concentration for E activation [#/uL] +p = Pars(10) ;% EA proliferation rate [/d] +fm = Pars(11) ;% fraction of EA that revert to memory +dK = Pars(12) ;% K death rate constant [/d] +d0 = Pars(13) ;% E0 death rate constant [/d] +dA = Pars(14) ;% EA death rate constant [/d] + +Xi = Pars(15) ;% X initial condition [pmol/kg] +ka = Pars(16) ;% N803 absorption rate constant [/d] +ke = Pars(17) ;% N803 elimination rate constant [/d] +vd = Pars(18) ;% N803 'volume of distribution'/'bioavailability' [L/kg] + +C50 = Pars(19) ;% N803 half-saturation concentration [pM] +rK1_dK = Pars(20) ;% K maximum growth rate [/d] +rE1_d0 = Pars(21) ;% E0 maximum growth rate [/d] +a1 = Pars(22) ;% E activation stimulation factor [] + +dR = Pars(23) ;% regulation decay rate constant [/d] +rK2 = Pars(24) ;% K proliferation regulation factor [] +rE2 = Pars(25) ;% E0 proliferation regulation factor [] +a2 = Pars(26) ;% E activation regulation factor [] + +N = Pars(27) ;% # of tolerance variables (delays + acting) +del = Pars(28) ;% tolerance rate constant [/d] +tau = Pars(29) ;% tolerance recovery [] +eta = Pars(30) ;% tolerance effect factor [] + +%% ------------------------------------------------------------------------ +% Declare initial conditions & calculated parameters ---------------------- + +% some parameters are solved for by assuming an initially steady state + +Vi = 10^(Vi-3) ;% V initial condition [#/無] +E0i = Ei*(1-fEA) ;% E0 initial condition [#/無] +EAi = Ei*( fEA) ;% EA initial condition [#/無] (sum of E1 to E8) + +F5i = Vi/(Vi+V50) / (1+a2) ;% initial E activation rate modifiers [] + +% solve for activation rate a [/d] and E1 to E8 initial conditions [#/無] +MAT = diag( [0 -p*ones(1,7) -dA] ) ; +MAT(1,2:9) = ones(1,8) ; +MAT(2,1) = F5i * E0i ; +MAT(3:9,2:8) = MAT(3:9,2:8) + diag( 2*p*ones(1,7) ) ; +vec = [ EAi ; zeros(8,1) ] ; +sol = MAT\vec ; + +a = sol(1) ;% E activation rate constant [/d] +pV = gE_EAi + gK_Ki ;% V proliferation rate constant [/d] +gE = gE_EAi/EAi ;% E killing rate constant [] +gK = gK_Ki/Ki ;% K killing rate constant [] +rK = (1+rK2) * (Km+Ki)/Km * dK ;% K proliferation rate constant [/d] +rE = (1+rE2) * (Em+Ei)/Em * ( d0 + a*F5i - fm*dA*sol(9)/E0i ) ;% E0 prolif +rK1 = rK1_dK/dK ;% K proliferation stimulation factor [] +rE1 = rE1_d0/d0 ;% E0 proliferation stimulation factor [] +psi = dR/F5i ;% R generation rate constant [/d] + +% declare initial conditions for Y (see 'system' function) +% V K E0 E1-8 X C R tolerance +Yic = [ Vi Ki E0i sol(2:9)' 0 0 1 zeros(1,N) ] ; + +%% ------------------------------------------------------------------------ +% Solve for each solution period ------------------------------------------ + +Y = zeros(Pieces(end),length(Yic)) ;% percolate output matrix + +for n = 1:length(Pieces)-1 + + solvertime = tic ;% timer (for catching a stalled execution) + + Piece = Pieces(n):Pieces(n+1) ;% index for solution piece + + [~,Y_TEMP] = ode15s(@system,SoluTimes(Piece),Yic) ;% solution piece + + if length(Piece) == 2 % exclude timepoints added by ode15s + Y_TEMP = [ Y_TEMP(1,:) ; Y_TEMP(end,:) ] ; + end + + Y(Piece,:) = Y_TEMP ; % add 'Y_TEMP' to full solution 'Y' + + Yic = Y_TEMP(end,:) ;% update ICs for next solution piece + Yic(12) = Xi ; +end + +%% ------------------------------------------------------------------------ +% 'system' function ------------------------------------------------------- + +function dY = system(~,Y) + + if toc(solvertime) >= 10 + error('ODE solver stalled.') + end + + V = Y(01) ;% SIV virions [#/無] + K = Y(02) ;% NK cells [#/無] + E0 = Y(03) ;% resting CD8+ T cells [#/無] + EA = Y(04:11) ;% active CD8+ T cells [#/無] + X = Y(12) ;% N803 at absorption site [pmol/kg] + C = Y(13) ;% N803 plasma concentration [pM] + R = Y(14) ;% immune regulation [] + F = terms(Y') ;% modifiers for 'rK','rE','a' + + dY = Y ;% dY/dt + + % calculate rates of change for V,K,E,X,C,R + dY(01) = pV*V - gK*K*V - gE*sum(EA)*V ;% dV/dt + dY(02) = rK*K *F(13) - dK*K ;% dK/dt + dY(03) = rE*E0*F(14) - d0*E0 + fm*dA*EA(8) - a*E0*F(15) ;% dE0/dt + dY(04) = - p*EA(1) + a*E0*F(15) ;% dE1/dt + dY(05:10) = 2*p*EA(1:6) - p*EA(2:7) ;% dE2/dt to dE7/dt + dY(11) = 2*p*EA(7) - dA*EA(8) ;% dE8/dt + dY(12) = - ka*X ;% dX/dt + dY(13) = ka*X/vd - ke*C ;% dC/dt + dY(14) = psi*F(15) - dR*R ;% dR/dt + + % calculate rates of change for tolerance + if N ~= 0 + dY(15) = del * ( F(1) - Y(15) ) ; + dY(16:end-2) = del * ( Y(15:end-3) - Y(16:end-2) ) ; + dY(end-1) = del * ( Y( end-2) - (1+tau)*Y( end-1) ) ; + dY(end) = del * ( tau*Y( end-1) - tau *Y( end ) ) ; + end +end + +%% ------------------------------------------------------------------------ +% 'terms' function -------------------------------------------------------- + +function F = terms(Y) + + V = Y(:,01) ;% SIV virions [#/無] + K = Y(:,02) ;% NK cells [#/無] + E = sum(Y(:,03:11),2) ;% CD8+ T cells [#/無] + C = Y(:,13) ;% N803 plasma concentration [pM] + R = Y(:,14) ;% immune regulation [] + T = Y(:,13+N)+Y(:,14+N) ;% drug tolerance [] + + % percolate vector for storing terms + F = zeros(size(Y,1),13) ; + + F(:,01) = C ./ (C50+C) ;% N803 effect [] + F(:,02) = 1 ./ (1+eta*T) ;% tolerance effect [] + F(:,03) = F(:,1) .* F(:,2) ;% N803 effect (with tolerance) [] + F(:,04) = 1 + rK1*F(:,3) ;% K proliferation stimulation [] + F(:,05) = 1 + rE1*F(:,3) ;% E0 proliferation stimulation [] + F(:,06) = 1 + a1*F(:,3) ;% activation stimulation [] + F(:,07) = 1 ./ (1+rK2*R) ;% K proliferation regulation [] + F(:,08) = 1 ./ (1+rE2*R) ;% E0 proliferation regulation [] + F(:,09) = 1 ./ (1+ a2*R) ;% activation regulation [] + F(:,10) = Km./ (Km +K) ;% K proliferation density dependence [] + F(:,11) = Em./ (Em +E) ;% E0 proliferation density dependence [] + F(:,12) = V ./ (V50+V) ;% activation antigen dependence [] + F(:,13) = F(:,4) .* F(:,7) .* F(:,10) ;% K proliferation modifier [] + F(:,14) = F(:,5) .* F(:,8) .* F(:,11) ;% E0 proliferation modifier [] + F(:,15) = F(:,6) .* F(:,9) .* F(:,12) ;% activation modifier [] + +end + +%% ------------------------------------------------------------------------ +% Prepare main outputs 'Y_OUT' ------------------------------------------- + +Y = abs(Y) ;% removing negatives resulting from numerical error + +% V,E,K [fold change] +Y_OUT = [ log10(Y(:,1)/Vi) ... V [log fold change] + , sum(Y(:,3:11),2)/Ei ... E [fold change] + , Y(:,2)/Ki ];% K [fold change] + +%% ------------------------------------------------------------------------ +% Add columns to 'Y_OUT' (if requested) ----------------------------------- + +if isempty(fullOut) ; return ; end + +Y_OUT = [ Y_OUT zeros(Pieces(end),20-3) ] ; + +Y_OUT(:,4) = Y(:,03)/E0i ;% E0 [fold change] +Y_OUT(:,5) = sum(Y(:,04:11),2)/EAi ;% EA [fold change] +Y_OUT(:,6) = Y(:,12) ;% X [pmol/kg] +Y_OUT(:,7) = Y(:,13) ;% C [pM] +Y_OUT(:,8) = Y(:,14) ;% R [] +Y_OUT(:,9) = (Y(:,13+N)+Y(:,14+N)) * (N~=0) ;% T [] + +F = terms(Y) ;% modifiers for 'rK','rE','a' + +Y_OUT(:,10) = F(:,01) ;% N803 effect [] +Y_OUT(:,11) = F(:,02) ;% tolerance effect [] +Y_OUT(:,12) = F(:,03) ;% N803 effect (with tolerance) [] +Y_OUT(:,13) = F(:,04) ;% K proliferation stimulation [] +Y_OUT(:,14) = F(:,05) ;% E0 proliferation stimulation [] +Y_OUT(:,15) = F(:,06) ;% activation stimulation [] +Y_OUT(:,16) = F(:,07)/F(1,07) ;% K proliferation regulation [fold change] +Y_OUT(:,17) = F(:,08)/F(1,08) ;% E0 proliferation regulation [fold change] +Y_OUT(:,18) = F(:,09)/F(1,09) ;% activation regulation [fold change] +Y_OUT(:,19) = F(:,10)/F(1,10) ;% K prolif density dependence [fold change] +Y_OUT(:,20) = F(:,11)/F(1,11) ;% E0 prolif density dependence [fold change] +Y_OUT(:,21) = F(:,12)/F(1,12) ;% activ antigen dependence [fold change] +Y_OUT(:,22) = F(:,13)/F(1,13) ;% K proliferation modifier [fold change] +Y_OUT(:,23) = F(:,14)/F(1,14) ;% E0 proliferation modifier [fold change] +Y_OUT(:,24) = F(:,15)/F(1,15) ;% activation modifier [fold change] + +% NK killing proportion +Y_OUT(:,25) = gK*Y(:,2) ./ ( gK*Y(:,2) + gE*sum(Y(:,04:11),2) ) ; + +end \ No newline at end of file