-
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
5 changed files
with
552 additions
and
13 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
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
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
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,254 @@ | ||
| %% N803_single_no_B.m - solves model for 2 cohorts without bystander T8 | ||
| % | ||
| % /--------------------------------------------------------------\ | ||
| % | Date: 07/01/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 [#/無] | ||
| % E0 = resting SIV-specific CD8+ T cells [#/無] | ||
| % Ea = active SIV-specific CD8+ T cells [#/無] | ||
| % B0 = resting bystander CD8+ T cells [#/無] | ||
| % Ba = active bystander 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 , initial conditions ] for cohort 1 (see code) | ||
| % PARS(2,:) = [ parameters , initial conditions ] for cohort 2 (see code) | ||
| % | ||
| %% ======================================================================== | ||
| % FUNCTION | ||
| % ======================================================================== | ||
| function [Y_OUT,PARS] = ... | ||
| N803_single_no_B(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 | ||
|
|
||
| nShared = 4 ;% number of shared parameters | ||
| nVaried = 19 ;% number of varied parameters | ||
|
|
||
| % Rename inputed parameters ----------------------------------------------- | ||
| Xi = AllPars(01) ;% X initial condition [pmol/kg] | ||
| ka = AllPars(02) ;% N803 absorption rate constant [/d] | ||
| ke = AllPars(03) ;% N803 elimination rate constant [/d] | ||
| vd = AllPars(04) ;% N803 'volume of distribution'/'bioavailability' [L/kg] | ||
|
|
||
| fE = 1 ;% initial frequency: E/(E+B) | ||
| psi = 0 ;% Ba/Ea killing rate ratio [gB0/gE0] | ||
| V50B = 0 ;% 50% viral stimulation saturation for B [#/mL] | ||
| mBn = 5 ;% normalized Ba reversion rate constant [] | ||
| pB = 1 ;% Ba proliferation rate constant [/d] | ||
| gB2 = 0 ;% initial B killing regulation [] | ||
| aB2 = 0 ;% initial B activation regulation [] | ||
|
|
||
| %% Do for each cohort (NOT indenting loop) ================================ | ||
| for c = 1:2 | ||
|
|
||
| n = nShared + nVaried*(c-1) ; | ||
| % Rename inputed parameters ----------------------------------------------- | ||
| Vi = AllPars(01+n) ;% V initial value [log(#/mL)] | ||
| EBi = AllPars(02+n) ;% E+B initial value [#/無] | ||
| fEA = AllPars(03+n) ;% initial frequency: Ea/E | ||
|
|
||
| q = AllPars(04+n) ;% V growth rate (if E+B were absent) [/d] | ||
| V50E = AllPars(05+n) ;% 50% viral stimulation saturation for E [#/mL] | ||
| mEn = AllPars(06+n) ;% normalized Ea reversion rate constant [] | ||
|
|
||
| EB50 = AllPars(07+n) ;% 50% E+B proliferation saturation [#/無] | ||
| pE = AllPars(08+n) ;% Ea proliferation rate constant [/d] | ||
| d = AllPars(09+n) ;% E0/B0 death rate constant [/d] | ||
| dA = AllPars(10+n) ;% Ea/Ba death rate constant [/d] | ||
|
|
||
| C50 = AllPars(11+n) ;% 50% N803 stimulation concentration [pM] (Cohort 1) | ||
| pm = AllPars(12+n) ;% E0/B0 maximum proliferation rate [] | ||
| aE1 = AllPars(13+n) ;% E activation stimulation factor [] | ||
| aB1 = AllPars(14+n) ;% B activation stimulation factor [] | ||
|
|
||
| sig = AllPars(15+n) ;% sB/sE regulation generation rate ratio | ||
| dR = AllPars(16+n) ;% R decay rate constant [/d] | ||
| gE2 = AllPars(17+n) ;% initial E killing regulation [] | ||
| p2 = AllPars(18+n) ;% initial E0/B0 proliferation regulation [] | ||
| aE2 = AllPars(19+n) ;% initial E activation regulation [] | ||
|
|
||
| %% ------------------------------------------------------------------------ | ||
| % Calculate some initial conditions & parameters -------------------------- | ||
|
|
||
| Vi = 10^(Vi-3) ;% V initial value [#/無] | ||
| V50E = V50E/1000 ;% 50% viral stimulation saturation for E [#/無] | ||
| V50B = V50B/1000 ;% 50% viral stimulation saturation for B [#/無] | ||
|
|
||
| % restrict mE and mB such that initial activation aE and aB are positive | ||
| UE = (2*pE/(pE+dA))^7 ; | ||
| UB = 2*pB/(pB+dA) ; | ||
| mE = mEn*dA/(UE-1) ;% Ea reversion rate constant [/d] | ||
| mB = mBn*dA/(UB-1) ;% Ba reversion rate constant [/d] | ||
|
|
||
| % solve for initial ratios below (based on active steady-state) | ||
| ZE = UE/(mE+dA) ; | ||
| for i = 1:7 | ||
| ZE = ZE + (2*pE)^(i-1)/(pE+dA)^i ;% EAi/aEi/E0i | ||
| end | ||
| ZB = 1/(pB+dA) + UB/(mB+dA) ;% BAi/aBi/B0i | ||
|
|
||
| WE = 1 ; | ||
| for i = 1:7 | ||
| WE = WE + (mE+dA)*(pE+dA)^(i-1)/(2*pE)^i ;% EAi/E8i | ||
| end | ||
| WB = 1 + (mB+dA)/(2*pB) ;% BAi/B2i | ||
|
|
||
| QE = mE/WE - 1/ZE ;% collection | ||
| QB = mB/WB - 1/ZB ;% collection | ||
|
|
||
| fBA = 1/( 1 + QB/QE*(1-fEA)/fEA ) ;% initial frequency: Ba/B | ||
|
|
||
| % solve for E and B initial conditions | ||
| Ei = EBi*fE ;% initial E | ||
| Bi = EBi*(1-fE) ;% initial B | ||
| EA = Ei*fEA ;% initial Ea | ||
| E0 = Ei*(1-fEA) ;% initial E0 | ||
| BA = Bi*fBA ;% initial Ea | ||
| B0 = Bi*(1-fBA) ;% initial E0 | ||
|
|
||
| E = zeros(1,8) ;% initial E1-E8 | ||
| E(8) = EA/WE ;% E8 | ||
| E(7) = E(8) * (mE+dA)/(2*pE) ;% E7 | ||
| for i = 6:-1:1 | ||
| E(i) = E(i+1) * (pE+dA) / (2*pE) ;% E6 to E1 | ||
| end | ||
| B = BA/WB * [ (mB+dA) / (2*pB) , 1 ] ;% initial B1-B2 | ||
|
|
||
| % solve for rate constants | ||
| aEi = E(1) / E0 * (pE+dA) ;% initial activation rate for E0 | ||
| aBi = aEi * (UE*mE/(mE+dA)-1) / (UB*mB/(mB+dA)-1) ;% for B0 | ||
| aE0 = aEi * (V50E+Vi)/Vi * (1+aE2) ;% E0 activation rate constant [/d] | ||
| aB0 = aBi * (V50B+Vi)/Vi * (1+aB2) ;% B0 activation rate constant [/d] | ||
|
|
||
| pi = d - QE*EA/E0 ;% initial proliferation rate for E0/B0 | ||
| p0 = pi * (EB50+EBi)/EB50 * (1+p2) ;% E0/B0 proliferation rate con [/d] | ||
| p1 = pm/p0 ;% E0/B0 proliferation stimulation factor | ||
|
|
||
| gE0 = q / ( EA/(1+gE2) + psi*BA/(1+gB2) ) ;% Ea killing rate con [無/#-d] | ||
| gB0 = psi*gE0 ;% Ba killing rate constant [無/#-d] | ||
|
|
||
| sE = dR / ( Vi/(V50E+Vi) + sig*Vi/(V50B+Vi) ) ;% reg gen by E act | ||
| sB = sig*sE ;% regulation generation by B activation | ||
| Ri = 1 ;% initial regulation | ||
|
|
||
| %% ------------------------------------------------------------------------ | ||
| % Prepare parameter and initial value vectors and call 'N803_model_2' ----- | ||
|
|
||
| Pars(01) = q ;% V growth rate (if E+B were absent) [/d] | ||
| Pars(02) = gE0 ;% Ea killing rate constant [無/#-d] | ||
| Pars(03) = gB0 ;% Ba killing rate constant [無/#-d] | ||
|
|
||
| Pars(04) = V50E ;% 50% viral stimulation saturation for E [#/無] | ||
| Pars(05) = V50B ;% 50% viral stimulation saturation for B [#/無] | ||
| Pars(06) = aE0 ;% E0 activation rate constant [/d] | ||
| Pars(07) = aB0 ;% B0 activation rate constant [/d] | ||
| Pars(08) = mE ;% Ea reversion rate constant [/d] | ||
| Pars(09) = mB ;% Ba reversion rate constant [/d] | ||
|
|
||
| Pars(10) = EB50 ;% 50% E+B proliferation saturation [#/無] | ||
| Pars(11) = p0 ;% E0/B0 proliferation rate constant [/d] | ||
| Pars(12) = pE ;% Ea proliferation rate constant [/d] | ||
| Pars(13) = pB ;% Ba proliferation rate constant [/d] | ||
| Pars(14) = d ;% E0/B0 death rate constant [/d] | ||
| Pars(15) = dA ;% Ea/Ba 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 ;% E0/B0 proliferation stimulation factor [] | ||
| Pars(22) = aE1 ;% E activation stimulation factor [] | ||
| Pars(23) = aB1 ;% B activation stimulation factor [] | ||
|
|
||
| Pars(24) = sE ;% R generation due to E0 activation [/d] | ||
| Pars(25) = sB ;% R generation due to B0 activation [/d] | ||
| Pars(26) = dR ;% R decay rate constant [/d] | ||
| Pars(27) = gE2 ;% E killing regulation factor [] | ||
| Pars(28) = gB2 ;% B killing regulation factor [] | ||
| Pars(29) = p2 ;% E0/B0 proliferation regulation factor [] | ||
| Pars(30) = aE2 ;% E activation regulation factor [] | ||
| Pars(31) = aB2 ;% B activation regulation factor [] | ||
|
|
||
| % V E0-8 B0-2 X C R initial values | ||
| Yic = [ Vi E0 E B0 B 0 0 Ri Ri ] ; | ||
|
|
||
| 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<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 |
Oops, something went wrong.