Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@cody2
cody2 authored Feb 17, 2022
1 parent 69a7663 commit b4022cb
Showing 1 changed file with 81 additions and 72 deletions.
153 changes: 81 additions & 72 deletions N803_model_2.m
View file
Original file line number Diff line number Diff line change
@@ -1,19 +1,19 @@
%% N803_model_2.m - solves ODE model for N-803 treatment of SIV
%
% Extra outputs are absolute instead of fold-change (otherwise identical)
%
% /--------------------------------------------------------------\
% | Date: 11/05/2021 |
% | Date: 02/14/2022 |
% | Author: Jonathan Cody |
% | Affiliation: Purdue University |
% | Weldon School of Biomedical Engineering |
% | Pienaar Computational Systems Pharmacology Lab |
% \--------------------------------------------------------------/
%
% Version 2D+ (also allowing regulation calibration)
%
% 12/03/2021 - added some extra terms to the full output matrix
% 11/27/2021 - replaced mu with fA -> initial frequency: (Ea+Ba)/(E+B)
% Version 2G with expanded outputs
%
% 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 [#/無]
Expand Down Expand Up @@ -56,16 +56,18 @@
% ========================================================================
%
% Y_OUT(:,1) = V at points in 'SoluTimes' [log fold change]
% Y_OUT(:,2) = E at points in 'SoluTimes' [fold change]
% Y_OUT(:,2) = T8 at points in 'SoluTimes' [fold change]
% Y_OUT(:,3) = R at points in 'SoluTimes' [fold change]
%
% If 'fullOut' is non-empty, more columns are included (see code)
%
% CalcPars = vector of calculated parameters
%
%% ========================================================================
% FUNCTION
% ========================================================================
function Y_OUT = N803_model_2(SoluTimes,DoseTimes,Pars,...
EquatedPars,fullOut,timeOut,warnOff)
function [Y_OUT,CalcPars] = N803_model_2(SoluTimes,DoseTimes,Pars,...
EquatedPars,fullOut,timeOut,warnOff)

% convert inputs to column vectors
if size(SoluTimes ,1)==1 ; SoluTimes = SoluTimes' ; end
Expand Down Expand Up @@ -196,76 +198,83 @@

Y_OUT = [ log10(Y(:,1)/Vi) ... V [log fold change]
, sum(Y(:,2:13),2)/EBi ... E [fold change]
, Y(:,15+nR)/Y(1,15+nR) ];% R [fold change]
, Y(:,15+nR) ];% R

% Calculated parameters
CalcPars = [mE mB gE0 gB0 p0 aE0 aB0 p1 p2 aE2 aB2 gE2 gB2 Ri] ;

%% ------------------------------------------------------------------------
% Add columns to 'Y_OUT' (if requested) -----------------------------------
% NOTE: ERC1 = effective first-order growth rate constant

if isempty(fullOut) ; return ; end

Y_OUT = [ Y_OUT zeros(Pieces(end),32) ] ;

E0 = Y(:,02) ;% resting specific CD8+ T cells [#/無]
B0 = Y(:,11) ;% resting bystander CD8+ T cells [#/無]
Ea = sum( Y(:,03:10) ,2) ;% active specific CD8+ T cells [#/無]
Ba = sum( Y(:,12:13) ,2) ;% active bystander T cells [#/無]
R = Y(:,15+nR) ;% immune regulation []
[F,T] = terms(Y) ;% modifiers for 'gE0','gB0','p0','aE0','aB0'

Y_OUT(:,04) = (E0+Ea) / (E0(1)+Ea(1)) ;% E [fold change]
Y_OUT(:,05) = (B0+Ba) / (B0(1)+Ba(1)) ;% B [fold change]
Y_OUT(:,06) = E0 / E0(1) ;% E0 [fold change]
Y_OUT(:,07) = B0 / B0(1) ;% B0 [fold change]
Y_OUT(:,08) = Ea / Ea(1) ;% Ea [fold change]
Y_OUT(:,09) = Ba / Ba(1) ;% Ba [fold change]
Y_OUT(:,10) = (E0+Ea)./(E0+B0+Ea+Ba) ;% E frequency in (E+B)
Y_OUT(:,11) = (Ea+Ba)./(E0+B0+Ea+Ba) ;% activated frequency in (E+B)
Y_OUT(:,12) = (Ea)./(E0+Ea) ;% activated frequency in E
Y_OUT(:,13) = (Ba)./(B0+Ba) ;% activated frequency in B
Y_OUT(:,14) = Y(:,14) ;% X [pmol/kg]
Y_OUT(:,15) = Y(:,15) ;% C [pM]
Y_OUT(:,16) = T ;% T []

Y_OUT(:,17) = F(:,01) ;% N803 effect []
Y_OUT(:,18) = F(:,02) ;% tolerance effect []
Y_OUT(:,19) = F(:,03) ;% N803 effect (with tolerance) []
Y_OUT(:,20) = log10( F(:,04) ) ;% C-lfc in p
Y_OUT(:,21) = log10( F(:,05) ) ;% C-lfc in aE
Y_OUT(:,22) = log10( ( F(:,06) + F(:,14) )./F(:,14) ) ;% C-lfc in aB
Y_OUT(:,23) = log10( F(:,07)/F(1,07) ) ;% R-lfc in p
Y_OUT(:,24) = log10( F(:,08)/F(1,08) ) ;% R-lfc in aE
Y_OUT(:,25) = log10( F(:,09)/F(1,09) ) ;% R-lfc in aB
Y_OUT(:,26) = log10( F(:,10)/F(1,10) ) ;% R-lfc in gE
Y_OUT(:,27) = log10( F(:,11)/F(1,11) ) ;% R-lfc in gB
Y_OUT(:,28) = log10( F(:,12)/F(1,12) ) ;% EB-lfc in p
Y_OUT(:,29) = log10( F(:,13)/F(1,13) ) ;% V-lfc in aE
Y_OUT(:,30) = log10( F(:,14)/F(1,14) ) ;% V-lfc in aB
Y_OUT(:,31) = log10( F(:,15)/F(1,15) ) ;% lfc in p
Y_OUT(:,32) = log10( F(:,16)/F(1,16) ) ;% lfc in aE
Y_OUT(:,33) = log10( F(:,17)/F(1,17) ) ;% lfc in aB
Y_OUT(:,34) = log10( F(:,18)/F(1,18) ) ;% lfc in regulation generation

% Ba regulation stimulation proportion
Y_OUT(:,35) = sB*aB0*(F(:,14)+F(:,06)) ./ F(:,18) ;

% lfc in killing
Kill = gB0*Ba./(1+gB2*R) + gE0*Ea./(1+gE2*R) ;
Y_OUT(:,36) = log10( Kill / Kill(1) ) ;

% Ba killing proportion
Y_OUT(:,37) = gB0*Ba./(1+gB2*R) ./ Kill ;

% Average E proliferation rate
Y_OUT(:,38) = ( p0*F(:,15).*E0 + pE*sum( Y(:,03:9) ,2) ) ./ (E0+Ea) ;

% Average B proliferation rate
Y_OUT(:,39) = ( p0*F(:,15).*B0 + pB*Y(:,12) ) ./ (B0+Ba) ;

% Average E killing rate
Y_OUT(:,40) = gE0*Ea./(1+gE2*R) ./ (E0+Ea) ;

% Average B killing rate
Y_OUT(:,41) = gB0*Ba./(1+gB2*R) ./ (B0+Ba) ;
E0 = Y(:,02) ;% resting specific CD8+ T cells [#/無]
B0 = Y(:,11) ;% resting bystander CD8+ T cells [#/無]
Ea = sum( Y(:,03:10) ,2) ;% active specific CD8+ T cells [#/無]
Ba = sum( Y(:,12:13) ,2) ;% active bystander T cells [#/無]
Eap = sum( Y(:,03:09) ,2) ;% Ea that is proliferating
Bap = Y(:,12) ;% Ba that is proliferating
Eaf = Y(:,10) ;% last generation of Ea
Baf = Y(:,13) ;% last generation of Ba
R = Y(:,15+nR) ;% immune regulation []
[F,T] = terms(Y) ;% see 'terms' function
gE = gE0*F(:,10) ;% specific killing rate [#/無-d]
gB = gB0*F(:,11) ;% bystander killing rate [#/無-d]
p = p0*F(:,15) ;% E0/B0 proliferation rate [/d]
aE = aE0*F(:,16) ;% specific activation rate [/d]
aB = aB0*F(:,17) ;% bystander activation rate [/d]
RgenE = sE*aE0* F(:,13).* F(:,05) ;% R gen by aE
RgenB = sB*aB0*( F(:,14) + F(:,06) ) ;% R gen by aB

Y_OUT(:,04) = Y(:,14) ;% X [pmol/kg]
Y_OUT(:,05) = Y(:,15) ;% C [pM]
Y_OUT(:,06) = T ;% T []
Y_OUT(:,07) = F(:,01) ;% N803 effect []
Y_OUT(:,08) = F(:,02) ;% tolerance effect []
Y_OUT(:,09) = F(:,03) ;% N803 effect (with tolerance) []

Y_OUT(:,10) = (E0+Ea) ;% E
Y_OUT(:,11) = (B0+Ba) ;% B
Y_OUT(:,12) = (E0+Ea)./(E0+B0+Ea+Ba) ;% E frequency in (E+B)
Y_OUT(:,13) = (Ea+Ba)./(E0+B0+Ea+Ba) ;% activated frequency in (E+B)
Y_OUT(:,14) = (Ea)./(E0+Ea) ;% activated frequency in E
Y_OUT(:,15) = (Ba)./(B0+Ba) ;% activated frequency in B

Y_OUT(:,16) = E0 ;% E0
Y_OUT(:,17) = B0 ;% B0
Y_OUT(:,18) = p - d ;% E0/B0 expansion ERC1
Y_OUT(:,19) = mE*Eaf./E0 - aE ;% E0 activation ERC1
Y_OUT(:,20) = mB*Baf./B0 - aB ;% B0 activation ERC1
Y_OUT(:,21) = log10( F(:,04) ) ;% C multiplier for p* (log)
Y_OUT(:,22) = log10( F(:,07) ) ;% R multiplier for p* (log)
Y_OUT(:,23) = log10( F(:,12) ) ;% density-dep multiplier for p* (log)

Y_OUT(:,24) = Ea ;% Ea
Y_OUT(:,25) = Ba ;% Ba
Y_OUT(:,26) = (pE*Eap - dA*Ea) ./ Ea ;% EA expansion ERC1
Y_OUT(:,27) = (pB*Bap - dA*Ba) ./ Ba ;% BA expansion ERC1
Y_OUT(:,28) = (aE.*E0 - mE*Eaf) ./ Ea ;% EA activation ERC1
Y_OUT(:,29) = (aB.*B0 - mB*Baf) ./ Ba ;% BA activation ERC1
Y_OUT(:,30) = log10( F(:,05) ) ;% C multiplier for aE* (log)
Y_OUT(:,31) = log10( ( F(:,06) + F(:,14) )./F(:,14) ) ;% same for aB* (log)
Y_OUT(:,32) = log10( F(:,08) ) ;% R multiplier for aE* (log)
Y_OUT(:,33) = log10( F(:,09) ) ;% R multiplier for aB* (log)
Y_OUT(:,34) = log10( F(:,13) ) ;% V multiplier for aE* (log)
Y_OUT(:,35) = log10( F(:,14) ) ;% V multiplier for aB* (log)

Y_OUT(:,36) = q - gE.*Ea - gB.*Ba ;% V ERC1
Y_OUT(:,37) = log10( gE.*Ea ) ;% E killing ERC1 (log)
Y_OUT(:,38) = log10( gB.*Ba ) ;% B killing ERC1 (log)
Y_OUT(:,39) = log10( F(:,10) ) ;% R multiplier for gE* (log)
Y_OUT(:,40) = log10( F(:,11) ) ;% R multiplier for gB* (log)
Y_OUT(:,41) = gB.*Ba ./ (gE.*Ea + gB.*Ba) ;% B killing fraction

Y_OUT(:,42) = log10( RgenE ) ;% R gen by aE (log)
Y_OUT(:,43) = log10( RgenB ) ;% R gen by aB (log)
Y_OUT(:,44) = RgenB ./ (RgenE + RgenB) ;% R gen by aB fraction

%% ========================================================================
% NESTED FUNCTIONS
Expand Down Expand Up @@ -305,8 +314,8 @@
C = Y(15) ;% N803 plasma concentration [pM]
F = terms(Y') ;% modifiers for 'gE0','gB0','p0','aE0','aB0'

gE = gE0*F(10) ;% specific killing rate [/d]
gB = gB0*F(11) ;% bystander killing rate [/d]
gE = gE0*F(10) ;% specific killing rate [#/無-d]
gB = gB0*F(11) ;% bystander killing rate [#/無-d]
p = p0*F(15) ;% E0/B0 proliferation rate [/d]
aE = aE0*F(16) ;% specific activation rate [/d]
aB = aB0*F(17) ;% bystander activation rate [/d]
Expand Down

0 comments on commit b4022cb

Please sign in to comment.