-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
1 changed file
with
259 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,259 @@ | ||
| %% N803_shared_noN.m - solves model for 2 cohorts with shared parameters | ||
| % | ||
| % Edits of 'N803_shared.m' to zero out non-SIV-specific CD8+ T cells noted !! | ||
| % | ||
| %#ok<*NASGU> suppressing 'variable is not used' | ||
| % | ||
| % /--------------------------------------------------------------\ | ||
| % | 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_noN(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 = -c/b ;% ( -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([1:6 8 10:30 ]) 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<AGROW> % 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 |