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 Jan 8, 2023
1 parent 98022fd commit 7033560
Show file tree
Hide file tree
Showing 6 changed files with 214 additions and 101 deletions.
109 changes: 54 additions & 55 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
%
% /--------------------------------------------------------------\
% | Date: 12/25/2022 |
% | Date: 01/08/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 [#/無]
% Nomenclature: V = SIV virions [#/uL]
% T8 = total CD8+ T cells [#/uL]
% S0 = resting SIV-specific CD8+ T cells [#/uL]
% Sa = active SIV-specific CD8+ T cells [#/uL]
% N0 = resting non-SIV-specific CD8+ T cells [#/uL]
% Na = active non-SIV-specific CD8+ T cells [#/uL]
% X = N803 at absorption site [pmol/kg]
% C = N803 plasma concentration [pM]
% R = regulation [] (dimensionless quantity)
Expand Down Expand Up @@ -94,17 +94,17 @@
% Rename inputed parameters -----------------------------------------------

q = Pars(01) ;% V growth rate (if S+N were absent) [/d]
gS0 = Pars(02) ;% Sa killing rate constant [無/#-d]
gN0 = Pars(03) ;% Na killing rate constant [無/#-d]
gS0 = Pars(02) ;% Sa killing rate constant [uL/#-d]
gN0 = Pars(03) ;% Na killing rate constant [uL/#-d]

V50S = Pars(04) ;% 50% viral stimulation saturation for S [#/無]
V50N = Pars(05) ;% 50% viral stimulation saturation for N [#/無]
V50S = Pars(04) ;% 50% viral stimulation saturation for S [#/uL]
V50N = Pars(05) ;% 50% viral stimulation saturation for N [#/uL]
aS0 = Pars(06) ;% S0 activation rate constant [/d]
aN0 = Pars(07) ;% N0 activation rate constant [/d]
mS = Pars(08) ;% Sa reversion rate constant [/d]
mN = Pars(09) ;% Na reversion rate constant [/d]

SN50 = Pars(10) ;% 50% S+N proliferation saturation [#/無]
SN50 = Pars(10) ;% 50% S+N proliferation saturation [#/uL]
p0 = Pars(11) ;% S0/N0 proliferation rate constant [/d]
pS = Pars(12) ;% Sa proliferation rate constant [/d]
pN = Pars(13) ;% Na proliferation rate constant [/d]
Expand Down Expand Up @@ -184,11 +184,11 @@

Y_OUT = [ Y_OUT, Y, zeros(Pieces(end),46) ] ;

V = Y(:,01) ;% SIV virions [#/無]
S0 = Y(:,02) ;% resting specific CD8+ T cells [#/無]
N0 = Y(:,11) ;% resting non-SIV-specific CD8+ T cells [#/無]
Sa = sum( Y(:,03:10) ,2) ;% active specific CD8+ T cells [#/無]
Na = sum( Y(:,12:13) ,2) ;% active non-SIV-specific T cells [#/無]
V = Y(:,01) ;% SIV virions [#/uL]
S0 = Y(:,02) ;% resting specific CD8+ T cells [#/uL]
N0 = Y(:,11) ;% resting non-SIV-specific CD8+ T cells [#/uL]
Sa = sum( Y(:,03:10) ,2) ;% active specific CD8+ T cells [#/uL]
Na = sum( Y(:,12:13) ,2) ;% active non-SIV-specific T cells [#/uL]
Saf = Y(:,10) ;% last generation of Sa
Naf = Y(:,13) ;% last generation of Na
Sap = Sa - Saf ;% Sa that is proliferating
Expand All @@ -198,12 +198,11 @@
p = F(:,13) ;% S0/N0 proliferation rate [/d]
aS = F(:,14) ;% S0 activation rate [/d]
aN = F(:,15) ;% N0 activation rate [/d]
gS = F(:,16) ;% Sa killing rate [#/無-d]
gN = F(:,17) ;% Na killing rate [#/無-d]
gS = F(:,16) ;% Sa killing rate [#/uL-d]
gN = F(:,17) ;% Na killing rate [#/uL-d]

Y_OUT(:,20) = F(:,01) ;% N803 effect []

Y_OUT(:,21) = (S0+N0+Sa+Na) ;% CD8+ T cells [#/無]
Y_OUT(:,21) = (S0+N0+Sa+Na) ;% CD8+ T cells [#/uL]
Y_OUT(:,22) = Sa ;% Sa
Y_OUT(:,23) = Na ;% Na
Y_OUT(:,24) = (S0+Sa) ;% S
Expand All @@ -214,35 +213,35 @@
Y_OUT(:,29) = log10( Sa ) ;% Sa (log)
Y_OUT(:,30) = log10( Na ) ;% Na (log)

Y_OUT(:,31) = p.*S0 ;% S0 proliferation term [#/無/d]
Y_OUT(:,32) = p.*N0 ;% N0 proliferation term [#/無/d]
Y_OUT(:,33) = d.*S0 ;% S0 death term [#/無/d]
Y_OUT(:,34) = d.*N0 ;% N0 death term [#/無/d]
Y_OUT(:,31) = p.*S0 ;% S0 proliferation term [#/uL/d]
Y_OUT(:,32) = p.*N0 ;% N0 proliferation term [#/uL/d]
Y_OUT(:,33) = d.*S0 ;% S0 death term [#/uL/d]
Y_OUT(:,34) = d.*N0 ;% N0 death term [#/uL/d]
Y_OUT(:,35) = p ;% S0/N0 proliferation rate [/d]
Y_OUT(:,36) = F(:,02) ;% C multiplier for p* (log)
Y_OUT(:,37) = F(:,05) ;% R multiplier for p* (log)
Y_OUT(:,38) = F(:,10) ;% density-dep multiplier for p* (log)

Y_OUT(:,39) = aS.*S0 ;% S0 activation term [#/無/d]
Y_OUT(:,40) = aN.*N0 ;% N0 activation term [#/無/d]
Y_OUT(:,41) = pS*Sap ;% Sa proliferation term [#/無/d]
Y_OUT(:,42) = pN*Nap ;% Na proliferation term [#/無/d]
Y_OUT(:,43) = dA*Sa ;% Sa death term [#/無/d]
Y_OUT(:,44) = dA*Na ;% Na death term [#/無/d]
Y_OUT(:,45) = mS*Saf ;% Sa reversion term [#/無/d]
Y_OUT(:,46) = mN*Naf ;% Na reversion term [#/無/d]
Y_OUT(:,39) = aS.*S0 ;% S0 activation term [#/uL/d]
Y_OUT(:,40) = aN.*N0 ;% N0 activation term [#/uL/d]
Y_OUT(:,41) = pS*Sap ;% Sa proliferation term [#/uL/d]
Y_OUT(:,42) = pN*Nap ;% Na proliferation term [#/uL/d]
Y_OUT(:,43) = dA*Sa ;% Sa death term [#/uL/d]
Y_OUT(:,44) = dA*Na ;% Na death term [#/uL/d]
Y_OUT(:,45) = mS*Saf ;% Sa reversion term [#/uL/d]
Y_OUT(:,46) = mN*Naf ;% Na reversion term [#/uL/d]
Y_OUT(:,47) = aS ;% S0 activation rate [/d]
Y_OUT(:,48) = aN ;% N0 activation rate [/d]
Y_OUT(:,49) = F(:,03) ;% C multiplier for aS*
Y_OUT(:,50) = ( F(:,04) + F(:,12) )./F(:,12) ;% C multiplier for aS*
Y_OUT(:,50) = ( F(:,04) + F(:,12) )./F(:,12) ;% C multiplier for aN*
Y_OUT(:,51) = F(:,06) ;% R multiplier for aS*
Y_OUT(:,52) = F(:,07) ;% R multiplier for aN*
Y_OUT(:,53) = F(:,11) ;% V multiplier for aS*
Y_OUT(:,54) = F(:,12) ;% V multiplier for aN*

Y_OUT(:,55) = q*V ;% V growth term [#/無/d]
Y_OUT(:,56) = gS.*Sa.*V ;% S killing term [#/無/d]
Y_OUT(:,57) = gN.*Na.*V ;% N killing term [#/無/d]
Y_OUT(:,55) = q*V ;% V growth term [#/uL/d]
Y_OUT(:,56) = gS.*Sa.*V ;% S killing term [#/uL/d]
Y_OUT(:,57) = gN.*Na.*V ;% N killing term [#/uL/d]
Y_OUT(:,58) = gS.*Sa ;% S killing rate [/d]
Y_OUT(:,59) = gN.*Na ;% N killing rate [/d]
Y_OUT(:,60) = F(:,08) ;% R multiplier for gS*
Expand Down Expand Up @@ -270,11 +269,11 @@
error('ODE solver stalled.')
end

V = Y(01) ;% SIV virions [#/無]
S0 = Y(02) ;% resting specific CD8+ T cells [#/無]
Sa = Y(03:10) ;% active specific CD8+ T cells [#/無] (S1 to S8)
N0 = Y(11) ;% resting non-SIV-specific CD8+ T cells [#/無]
Na = Y(12:13) ;% active non-SIV-specific T cells [#/無] (N1 to N2)
V = Y(01) ;% SIV virions [#/uL]
S0 = Y(02) ;% resting specific CD8+ T cells [#/uL]
Sa = Y(03:10) ;% active specific CD8+ T cells [#/uL] (S1 to S8)
N0 = Y(11) ;% resting non-SIV-specific CD8+ T cells [#/uL]
Na = Y(12:13) ;% active non-SIV-specific T cells [#/uL] (N1 to N2)
X = Y(14) ;% N803 at absorption site [pmol/kg]
C = Y(15) ;% N803 plasma concentration [pM]
R = Y(16:17) ;% immune regulation []
Expand All @@ -283,8 +282,8 @@
p = F(13) ;% S0/N0 proliferation rate [/d]
aS = F(14) ;% S0 activation rate [/d]
aN = F(15) ;% N0 activation rate [/d]
gS = F(16) ;% Sa killing rate [#/無-d]
gN = F(17) ;% Na killing rate [#/無-d]
gS = F(16) ;% Sa killing rate [#/uL-d]
gN = F(17) ;% Na killing rate [#/uL-d]
s = F(18)+F(19) ;% regulation generation rate [/d]

dY = Y ;% dY/dt
Expand All @@ -307,19 +306,19 @@

function J = jacob(~,Y)

V = Y(01) ;% SIV virions [#/無]
S0 = Y(02) ;% resting specific CD8+ T cells [#/無]
Sa = sum(Y(03:10)) ;% active specific CD8+ T cells [#/無]
N0 = Y(11) ;% resting non-SIV-specific CD8+ T cells [#/無]
Na = Y(12)+Y(13) ;% active non-SIV-specific T cells [#/無]
V = Y(01) ;% SIV virions [#/uL]
S0 = Y(02) ;% resting specific CD8+ T cells [#/uL]
Sa = sum(Y(03:10)) ;% active specific CD8+ T cells [#/uL]
N0 = Y(11) ;% resting non-SIV-specific CD8+ T cells [#/uL]
Na = Y(12)+Y(13) ;% active non-SIV-specific T cells [#/uL]
C = Y(15) ;% N803 plasma concentration [pM]
F = terms(Y') ;% modifiers for 'gS0','gN0','p0','aS0','aN0'

p = F(13) ;% S0/N0 proliferation rate [/d]
aS = F(14) ;% S0 activation rate [/d]
aN = F(15) ;% N0 activation rate [/d]
gS = F(16) ;% Sa killing rate [#/無-d]
gN = F(17) ;% Na killing rate [#/無-d]
gS = F(16) ;% Sa killing rate [#/uL-d]
gN = F(17) ;% Na killing rate [#/uL-d]

J = zeros(17) ;

Expand Down Expand Up @@ -400,8 +399,8 @@

function F = terms(Y)

V = Y(:,01) ;% SIV virions [#/無]
SN = sum(Y(:,02:13),2) ;% CD8+ T cells [#/無]
V = Y(:,01) ;% SIV virions [#/uL]
SN = sum(Y(:,02:13),2) ;% CD8+ T cells [#/uL]
C = Y(:,15) ;% N803 plasma concentration [pM]
R = Y(:,17) ;% immune regulation []

Expand All @@ -427,8 +426,8 @@
F(:,13) = p0 * F(:,10).* F(:,02).* F(:,05) ;% S0/N0 proliferation [/d]
F(:,14) = aS0* F(:,11).* F(:,03).* F(:,06) ;% S0 activation rate [/d]
F(:,15) = aN0* (F(:,12)+F(:,04)).* F(:,07) ;% N0 activation rate [/d]
F(:,16) = gS0* F(:,08) ;% Sa killing rate [#/無-d]
F(:,17) = gN0* F(:,09) ;% Na killing rate [#/無-d]
F(:,16) = gS0* F(:,08) ;% Sa killing rate [#/uL-d]
F(:,17) = gN0* F(:,09) ;% Na killing rate [#/uL-d]

F(:,18) = sS* F(:,11).* F(:,03) ;% regulation generation by aS [/d]
F(:,19) = sN* (F(:,12) + F(:,04)) ;% regulation generation by aN [/d]
Expand Down
114 changes: 114 additions & 0 deletions N803_plotter_2_parameters.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,114 @@
% N803_parameters.m - makes violin plots of parameter distributions

close all
clear all %#ok<CLALL>
addpath Analyzer
addpath Data

load('N803_MCMC','C')

%% Model Parameters =======================================================

% Retrieving parameter boundaries -----------------------------------------
warning off ; load('N803_shared_Results','Results') ; warning on
% q vS vN h pS pN d dA C50 aS1 aN1 dR p2 gN2
id1 = [ 7 9 10 13 15 16 17 18 23 25 26 27 29 30 ] ;
id2 = [ 1 4 5 10 12 13 14 15 16 18 19 22 25 24 ] ;
Bounds = NaN*ones(2,27) ;
Bounds(:,id2) = Results(1).setup.GuessBOUNDS(:,id1) ;
Bounds(:,4:5) = Bounds(:,4:5)/1000 ;% converting vS,vN to [#/uL]

% Plotting shared parameters ----------------------------------------------
PARS = C(1).ParSample(:,[1:15,20:31]) ;% parameter sample
Size = [6.5 3] ;% figure size
figName = 'Fig S1 Parameters' ;

f = figure('Name', figName, ...
'Color', [1 1 1], ...
'PaperUnits', 'Inches', ...
'PaperSize', Size, ...
'Units', 'Inches', ...
'Position', [0 0 Size]) ;

Names = {'q','g_S','g_N','V_{50,S}','V_{50,N}','a_S','a_N','m_S','m_N',...
'H','p','p_S','p_N','d','d_A','C_{50}','\rho','\alpha_S','\alpha_N',...
's_S','s_N','d_R','\lambda_S','\lambda_N','\phi','\zeta_S','\zeta_N'} ;

parplotter(PARS, ...
'AxesLims', [0 28 -3 4],...
'AxesPos', [.075 .2 .9 .75],...
'Bounds', Bounds,...
'Factors', [1 1e3 1e3 ones(1,24)],...converting gS0,gN0 to [nL/#-d]
'figHand', f,...
'logScale', 1,...
'Names', Names,...
'ViolColor', [.5 .5 .5]) ;

set(gca,'FontSize',10,'FontName','Arial','LineWidth',1,'XTickLabelRotation',315)
set(gcf,'InvertHardcopy','off')
print(figName,'-dtiff','-r300')

%% Sampled Parameter Convolutions =========================================
PARS = C(1).ParSample(:,1:31) ;% parameter sample

% Calculating sampled parameter convolutions ------------------------------
YIC1 = C(1).ParSample(:,32:end) ;% initial conditions
YIC2 = C(2).ParSample(:,32:end) ;% initial conditions
SN1 = sum( YIC1(:,2:13),2 ) ;% initial: S+N (co 1)
SN2 = sum( YIC2(:,2:13),2 ) ;% initial: S+N (co 2)
R1 = YIC1(:,end) ;% initial R (co 1)
R2 = YIC2(:,end) ;% initial R (co 2)
mS = PARS(:,08) ;% Sa reversion rate constant [/d]
mN = PARS(:,09) ;% Na reversion rate constant [/d]
h = PARS(:,10) ;% 50% S+N proliferation saturation [#/uL]
p0 = PARS(:,11) ;% S0/N0 proliferation rate constant [/d]
pS = PARS(:,12) ;% Sa proliferation rate constant [/d]
pN = PARS(:,13) ;% Na proliferation rate constant [/d]
d = PARS(:,14) ;% S0/N0 death rate constant [/d]
dA = PARS(:,15) ;% Sa/Na death rate constant [/d]
p1 = PARS(:,21) ;% S0/N0 proliferation stimulation factor []
p2 = PARS(:,29) ;% S0/N0 proliferation regulation factor []

psi = PARS(:,03)./PARS(:,02) ;% gN0/gS0
sig = PARS(:,25)./PARS(:,24) ;% sN/sS
fS1 = sum( YIC1(:,2:10),2 ) ./ SN1 ;% initial: S/(S+N) (co 1)
fS2 = sum( YIC2(:,2:10),2 ) ./ SN2 ;% initial: S/(S+N) (co 2)
fSn = (fS2-fS1)./ (.3-fS1) ;% normalized initial S/(S+N) (co 2)
US = 2*(2*pS./(pS+dA)).^7 ;
UN = 2* 2*pN./(pN+dA) ;
mSn = mS.*(US-1)./dA ;% normalized Sa reversion rate constant []
mNn = mN.*(UN-1)./dA ;% normalized Na reversion rate constant []
p0_max1 = d.*(1+p2.*R1).*(h+SN1)./h ;
p0_max2 = d.*(1+p2.*R2).*(h+SN2)./h ;
p0_max = min( [p0_max1 p0_max2],[],2 ) ;
p0n = p0./ p0_max ;% nomalized S0/N0 proliferation rate constant [/d]
pm = p1.*p0 ;% S0/N0 maximum proliferation rate []

% Plotting sampled parameter convolutions ---------------------------------
PARCONS = [ fS1 fSn psi mSn mNn p0n pm sig ] ;
id1 = [ 05 06 08 11 12 14 24 28 ] ;
Bounds = Results(1).setup.GuessBOUNDS(:,id1) ;
Size = [3 3*5/7] ;% figure size
figName = 'Fig S2 Convolutions' ;% figure name

Names = {'f_S','\xi_S','g_N/g_S','\mu_S','\mu_N','\psi','p\cdot\rho','s_N/s_S'} ;

f = figure('Name', figName, ...
'Color', [1 1 1], ...
'PaperUnits', 'Inches', ...
'PaperSize', Size, ...
'Units', 'Inches', ...
'Position', [0 0 Size]) ;

parplotter(PARCONS,...
'AxesLims', [0 9 -2 2],...
'AxesPos', [.2 .25 .75 .7],...
'Bounds', Bounds,...
'figHand', f,...
'logScale', 1,...
'Names', Names,...
'ViolColor', [.5 .5 .5]) ;

set(gca,'FontSize',10,'FontName','Arial','LineWidth',1,'XTickLabelRotation',315)
set(gcf,'InvertHardcopy','off')
print(figName,'-dtiff','-r300')
Loading

0 comments on commit 7033560

Please sign in to comment.