diff --git a/N803_shared.m b/N803_shared.m deleted file mode 100644 index 0510cb1..0000000 --- a/N803_shared.m +++ /dev/null @@ -1,255 +0,0 @@ -%% N803_shared.m - solves model for 2 cohorts with shared parameters -% -% /--------------------------------------------------------------\ -% | Date: 12/25/2022 | -% | Author: Jonathan Cody | -% | Affiliation: Purdue University | -% | Weldon School of Biomedical Engineering | -% | Pienaar Computational Systems Pharmacology Lab | -% \--------------------------------------------------------------/ -% -% Nomenclature: V = SIV virions [#/無] -% T8 = total CD8+ T cells [#/無] -% S0 = resting SIV-specific CD8+ T cells [#/無] -% Sa = active SIV-specific CD8+ T cells [#/無] -% N0 = resting non-SIV-specific CD8+ T cells [#/無] -% Na = active non-SIV-specific CD8+ T cells [#/無] -% X = N803 at absorption site [pmol/kg] -% C = N803 plasma concentration [pM] -% R = regulation [] (dimensionless quantity) -% -%% ======================================================================== -% INPUTS -% ======================================================================== -% -% SoluTimes = ascending vector of days at which to evaluate solution -% -% DoseTimes{c} = ascending vector of days at which to administer doses -% (elements of 'DoseTimes' must also be in 'SoluTimes') -% -% AllPars = vector of parameters (see list in function) -% -%% ======================================================================== -% OPTIONS -% ======================================================================== -% -% SkipTimes{c} = [min max] time point beyond which to skip model solving -% (outputs before 'min' will be made equal to output at 'min') -% (leave as [] to ignore and solve for all 'SoluTimes') -% -% oneCohort = scalar to run model for just one cohort ('1' or '2') -% (leave as [] to ignore and solve for both cohorts) -% -% All additional inputs will be passed as a cell vector to 'N803_model_2' -% and used to define options (see function for list) -% EX: N803_single(SoluTimes,DoseTimes,AllPars,'AbsTol',1e-2} -% will set ode solver absolute tolerance to 1e-2 -% -%% ======================================================================== -% OUTPUTS -% ======================================================================== -% -% Y_OUT(:,1) = V at points in 'SoluTimes' [log fold change] cohort 1 -% Y_OUT(:,2) = T8 at points in 'SoluTimes' [fold change] cohort 1 -% Y_OUT(:,3) = V at points in 'SoluTimes' [log fold change] cohort 2 -% Y_OUT(:,4) = T8 at points in 'SoluTimes' [fold change] cohort 2 -% -% PARS(1,:) = parameters for cohort 1 (see code) -% PARS(2,:) = parameters for cohort 2 (see code) -% -%% ======================================================================== -% FUNCTION -% ======================================================================== -function [Y_OUT,PARS] = ... - N803_shared(SoluTimes,DoseTimes,AllPars,SkipTimes,oneCohort,... - varargin) - -if isempty(oneCohort) ; RunCohort = [ 1 1 ] ; -elseif oneCohort == 1 ; RunCohort = [ 1 0 ] ; -else ; RunCohort = [ 0 1 ] ; -end - -Y_Cohort = cell(1,2) ;% cell for storing outputs -P_Cohort = cell(1,2) ;% cell for storing parmeters - -% Rename inputed parameters ----------------------------------------------- -Vi(1) = AllPars(01) ;% V initial value [log(#/mL)] (cohort 1) -Vi(2) = AllPars(02) ;% V initial value [log(#/mL)] (cohort 2) -SNi(1) = AllPars(03) ;% S+N initial value [#/無] (cohort 1) -SNi(2) = AllPars(04) ;% S+N initial value [#/無] (cohort 2) -fS(1) = AllPars(05) ;% initial frequency: S/(S+N) (cohort 1) -fSn = AllPars(06) ;% normalized S/(S+N) (co 2) -fS(2) = fS(1)+(0.3-fS(1))*fSn ;% initial frequency: S/(S+N) (cohort 2) - -q = AllPars(07) ;% V growth rate (if S+N were absent) [/d] -psi = AllPars(08) ;% ratio of base killing rates gN0/gS0 -V50S = AllPars(09) ;% 50% viral stimulation saturation for S [#/mL] -V50N = AllPars(10) ;% 50% viral stimulation saturation for N [#/mL] -mSn = AllPars(11) ;% normalized Sa reversion rate constant [] -mNn = AllPars(12) ;% normalized Na reversion rate constant [] - -SN50 = AllPars(13) ;% 50% S+N proliferation saturation [#/無] -p0n = AllPars(14) ;% nomalized S0/N0 proliferation rate constant [/d] -pS = AllPars(15) ;% Sa proliferation rate constant [/d] -pN = AllPars(16) ;% Na proliferation rate constant [/d] -d = AllPars(17) ;% S0/N0 death rate constant [/d] -dA = AllPars(18) ;% Sa/Na death rate constant [/d] - -Xi = AllPars(19) ;% X initial condition [pmol/kg] -ka = AllPars(20) ;% N803 absorption rate constant [/d] -ke = AllPars(21) ;% N803 elimination rate constant [/d] -vd = AllPars(22) ;% N803 'volume of distribution'/'bioavailability' [L/kg] -C50 = AllPars(23) ;% 50% N803 stimulation concentration [pM] (Cohort 1) -pm = AllPars(24) ;% S0/N0 maximum proliferation rate [] -aS1 = AllPars(25) ;% S activation stimulation factor [] -aN1 = AllPars(26) ;% N activation stimulation factor [] - -dR = AllPars(27) ;% R decay rate constant [/d] -sig = AllPars(28) ;% ratio of initial regulation due to N/S -p2 = AllPars(29) ;% S0/N0 proliferation regulation factor [] -gN2 = AllPars(30) ;% N killing regulation factor [] (cohort 1) - -%% ------------------------------------------------------------------------ -% Calculate some initial conditions & parameters -------------------------- - -Vi = 10.^(Vi - 3) ;% converting V initial value [#/無] -V50S = V50S/1000 ;% 50% viral stimulation saturation for S [#/無] -V50N = V50N/1000 ;% 50% viral stimulation saturation for N [#/無] -YS = Vi ./ (Vi+ V50S) ;% initial V/(V50S+V) [cohort 1,2] -YN = Vi ./ (Vi+ V50N) ;% initial V/(V50N+V) [cohort 1,2] - -% calculate initial R, and sS,sN -z = YS + sig*YN ; -Ri = [ 1 , (z(2)/z(1)) ] ;% initial regulation [cohort 1,2] -R = Ri(2) ;% initial regulation (cohort 2) -sS = dR / ( YS(1) + sig*YN(1) ) ;% R generation due to S0 activation [/d] -sN = sig*sS ;% R generation due to N0 activation [/d] - -% restrict mS and mN such that initial activation aS and aN are positive -US = 2*(2*pS/(pS+dA))^7 ; -UN = 2* 2*pN/(pN+dA) ; -mS = mSn*dA/(US-1) ;% Sa reversion rate constant [/d] -mN = mNn*dA/(UN-1) ;% Na reversion rate constant [/d] - -% restrict p0 to ensure same sign for S0/SA and for N0/NA -p0_max = min( d*(1+p2*Ri).*(SN50+SNi)/SN50 ) ;% maximum value for p0 -p0 = p0n * p0_max ;% S0/N0 prolif rate constant [/d] -p1 = pm/p0 ;% S0/N0 proliferation stimulation factor -pi = p0*SN50./(SN50+SNi)./(1+p2*Ri) ;% initial S0/N0 prolif rate [/d] - -% calculate aS0,aN0 and aS2,aN2 -aSi = (d-pi)/( mS*US/(dA+mS) - 1 ) ;% initial S0 activation rate [/d] -aNi = (d-pi)/( mN*UN/(dA+mN) - 1 ) ;% initial N0 activation rate [/d] -z = aSi./YS ; -aS2 = ( z(1)-z(2) ) / ( R*z(2)-z(1) ) ;% S0 activation rate constant [/d] -aS0 = aSi(1) / YS(1) * (1+aS2) ;% S activation regulation factor [] -z = aNi./YN ; -aN2 = ( z(1)-z(2) ) / ( R*z(2)-z(1) ) ;% N0 activation rate constant [/d] -aN0 = aNi(1) / YN(1) * (1+aN2) ;% N activation regulation factor [] - -% solve for initial ratios below (based on active steady-state) -ZS = US/(mS+dA) ; -for i = 1:7 - ZS = ZS + 2*(2*pS)^(i-1)/(pS+dA)^i ;% SAi/aSi/S0i -end -ZN = UN/(mN+dA) + 2/(pN+dA) ;% NAi/aNi/N0i - -% solve for initial values of S0,Sa,N0,Na each cohort -Si = SNi.*fS ;% initial S -Ni = SNi.*(1-fS) ;% initial N -S0 = Si./(1+ZS*aSi) ;% initial S0 -N0 = Ni./(1+ZN*aNi) ;% initial N0 -SA = S0.*(ZS*aSi) ;% initial SA -NA = N0.*(ZN*aNi) ;% initial NA - -% calculate gS0,gN0 and gS2 -beta = psi * NA ./ (1+gN2*Ri) ; z = beta(1) - beta(2) ; -a = z*R ; b = SA(1)*R - SA(2) + z*(1+R) ; c = SA(1) - SA(2) + z ; -gS2 = ( -b + sqrt(b^2 - 4*a*c) ) / (2*a) ;% S killing regulation factor [] -gS0 = q / ( SA(1)/(1+gS2) + beta(1) ) ;% Sa killing rate constant [無/#-d] -gN0 = psi * gS0 ;% Na killing rate constant [無/#-d] - -%% Do for each cohort (NOT indenting loop) ================================ -for c = 1:2 - -% solve for initial S1-8 and N1-2 -S = zeros(1,8) ;% initial S1-S8 -S(1) = 2*aSi(c)*S0(c) / (dA+pS) ;% S1 -for i = 2:7 - S(i) = 2*pS*S(i-1) / (dA+pS) ;% S2 to S7 -end -S(8) = 2*pS*S(7) / (dA+mS) ;% S8 - -N(1) = 2*aNi(c)*N0(c) / (dA+pN) ;% N1 -N(2) = 2*pN*N(1) / (dA+mN) ;% N2 - -%% ------------------------------------------------------------------------ -% Prepare parameter and initial value vectors and call 'N803_model_2' ----- - -Pars(01) = q ;% V growth rate (if S+N were absent) [/d] -Pars(02) = gS0 ;% Sa killing rate constant [無/#-d] -Pars(03) = gN0 ;% Na killing rate constant [無/#-d] - -Pars(04) = V50S ;% 50% viral stimulation saturation for S [#/無] -Pars(05) = V50N ;% 50% viral stimulation saturation for N [#/無] -Pars(06) = aS0 ;% S0 activation rate constant [/d] -Pars(07) = aN0 ;% N0 activation rate constant [/d] -Pars(08) = mS ;% Sa reversion rate constant [/d] -Pars(09) = mN ;% Na reversion rate constant [/d] - -Pars(10) = SN50 ;% 50% S+N proliferation saturation [#/無] -Pars(11) = p0 ;% S0/N0 proliferation rate constant [/d] -Pars(12) = pS ;% Sa proliferation rate constant [/d] -Pars(13) = pN ;% Na proliferation rate constant [/d] -Pars(14) = d ;% S0/N0 death rate constant [/d] -Pars(15) = dA ;% Sa/Na death rate constant [/d] - -Pars(16) = Xi ;% X initial condition [pmol/kg] -Pars(17) = ka ;% N803 absorption rate constant [/d] -Pars(18) = ke ;% N803 elimination rate constant [/d] -Pars(19) = vd ;% N803 'volume of distribution'/'bioavailability' [L/kg] -Pars(20) = C50 ;% 50% N803 stimulation concentration [pM] -Pars(21) = p1 ;% S0/N0 proliferation stimulation factor [] -Pars(22) = aS1 ;% S activation stimulation factor [] -Pars(23) = aN1 ;% N activation stimulation factor [] - -Pars(24) = sS ;% R generation due to S0 activation [/d] -Pars(25) = sN ;% R generation due to N0 activation [/d] -Pars(26) = dR ;% R decay rate constant [/d] -Pars(27) = gS2 ;% S killing regulation factor [] -Pars(28) = gN2 ;% N killing regulation factor [] -Pars(29) = p2 ;% S0/N0 proliferation regulation factor [] -Pars(30) = aS2 ;% S activation regulation factor [] -Pars(31) = aN2 ;% N activation regulation factor [] - -% V S0-8 N0-2 X C R initial values -Yic = [ Vi(c) S0(c) S N0(c) N 0 0 Ri(c) Ri(c) ] ; - -% if any( [ Pars Yic ] < 0 ) -% error('Negative parameters or initial values.') -% end - -% If 'SkipTimes' is empty, do not skip any times -if isempty(SkipTimes{c}) ; SkipTimes{c} = [-inf inf] ; end - -idLo = SoluTimes < SkipTimes{c}(1) ;% index of early times to skip soln -idHi = SoluTimes > SkipTimes{c}(2) ;% index of later times to skip soln -idSol = ~ ( idLo | idHi ) ;% index of times in 'SoluTimes' to solve - -if RunCohort(c) == 1 - Y_COH = N803_model_2(SoluTimes(idSol),DoseTimes{c},Pars,Yic,varargin) ; - Y_LO = ones(sum(idLo),1)*Y_COH(1 ,:) ;% constant Y for early times - Y_HI = ones(sum(idHi),1)*Y_COH(end,:) ;% constant Y for later times - Y_COH = [ Y_LO ; Y_COH ; Y_HI ] ;%#ok % total 'solution' matrix - Y_Cohort{c} = Y_COH ; - P_Cohort{c} = [ Pars Yic ] ; -end - -end - -%% Prepare main outputs 'Y_OUT' and 'PARS' ================================ - -Y_OUT = [ Y_Cohort{1} , Y_Cohort{2} ] ; -PARS = [ P_Cohort{1} ; P_Cohort{2} ] ; - -end \ No newline at end of file