Add files via upload

N803_shared_noN.m

@@ -1,7 +1,7 @@
 %% N803_shared_noN.m - solves model 2S for 2 cohorts with shared parameters
 %
 %     /--------------------------------------------------------------\
-%     | Date:         12/27/2022                                     |
+%     | Date:         01/01/2023                                     |
 %     | Author:       Jonathan Cody                                  |
 %     | Affiliation:  Purdue University                              |
 %     |               Weldon School of Biomedical Engineering        |
@@ -76,12 +76,12 @@
 % 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)
-Ei(1) = AllPars(03) ;% E+B initial value [#/無] (cohort 1)
-Ei(2) = AllPars(04) ;% E+B initial value [#/無] (cohort 2)
+Ei(1) = AllPars(03) ;% E initial value [#/無] (cohort 1)
+Ei(2) = AllPars(04) ;% E initial value [#/無] (cohort 2)
 
 q    = AllPars(05)  ;% V growth rate (if E were absent) [/d]
 V50E = AllPars(06)  ;% 50% viral stimulation saturation for E [#/mL]
-V50R = AllPars(07)  ;% 50% viral stimulation saturation for B [#/mL]
+V50R = AllPars(07)  ;% 50% viral stimulation saturation for R [#/mL]
 m    = AllPars(08)  ;% EA reversion rate constant [/d]
 
 E50  = AllPars(09)  ;% 50% E proliferation saturation [#/無]
@@ -101,34 +101,33 @@
 
 dR   = AllPars(22)  ;% R decay rate constant [/d]
 p2   = AllPars(23)  ;% E0 proliferation regulation factor []
-g2   = AllPars(24)  ;% EA killing regulation factor []
 
 %% ------------------------------------------------------------------------
 % Calculate some initial conditions & parameters --------------------------
 
 Vi = 10.^(Vi - 3) ;% converting V initial value [#/無]
 V50E = V50E/1000  ;% 50% viral stimulation saturation for E [#/無]
-V50R = V50R/1000  ;% 50% viral stimulation saturation for B [#/無]
+V50R = V50R/1000  ;% 50% viral stimulation saturation for R [#/無]
 YE = Vi ./ (Vi+ V50E) ;% initial V/(V50E+V) [cohort 1,2]
 YR = Vi ./ (Vi+ V50R) ;% initial V/(V50B+V) [cohort 1,2]
 
 % calculate initial R (normalized to Cohort 1) and s
 s0 = dR / YR(1)  ;% R generation rate [/d]
-Ri = [ 1 , s0*YR(2)/dR ]  ;% initial regulation [cohort 1,2]
+R = [ 1 , s0*YR(2)/dR ]  ;% initial regulation [cohort 1,2]
 
 % restrict p0 to ensure same sign for E0/EA and for B0/BA
-p0_max = min( d*(1+p2*Ri).*(E50+Ei)/E50 )  ;% maximum value for p0
+p0_max = min( d*(1+p2*R).*(E50+Ei)/E50 )  ;% maximum value for p0
 p0 = p0n * p0_max  ;% E0/B0 prolif rate constant [/d]
 p1 = pm/p0         ;% E0/B0 proliferation stimulation factor 
-pi = p0*E50./(E50+Ei)./(1+p2*Ri)  ;% initial E0/B0 prolif rate [/d]
+pi = p0*E50./(E50+Ei)./(1+p2*R)  ;% initial E0/B0 prolif rate [/d]
 
 % calculate initial ratio: E8/E0/a
 U = 2/(m+dA) * (2*pA/(pA+dA))^7 ;
 
 % calculate a0, a2
 ai = (d-pi)/( m*U - 1 )  ;% initial E0 activation rate [/d]
 F = ai./YE  ;
-a2 = ( F(1)-F(2) ) / ( Ri(2)*F(2)-F(1) );% E0 activation rate constant [/d]
+a2 = ( F(1)-F(2) ) / ( R(2)*F(2)-F(1) );% E0 activation rate constant [/d]
 a0 = ai(1) / YE(1) * (1+a2)           ;% E0 activation regulation factor []
 
 % calculate initial ratio of EA/E8
@@ -141,7 +140,8 @@
 E0 = Ei.*(m/W)./ ( (m/W) - pi + d + ai )  ;% initial E0
 EA = Ei - E0  ;% initial EA
 
-% calculate g0
+% calculate g0 and g2
+g2 = (EA(1)-EA(2)) / (EA(2)-EA(1)*R(2))  ;% EA killing regulation factor []
 g0 = q * (1+g2) / EA(1)  ;% EA killing rate constant [無/#-d]
 
 %% Do for each cohort (NOT indenting loop) ================================
@@ -187,7 +187,7 @@
 Pars(24) = a2    ;% E0 activation regulation factor []
 
 %       V     E0-8    X C R       initial values
-Yic = [ Vi(c) E0(c) E 0 0 Ri(c) Ri(c) ] ;
+Yic = [ Vi(c) E0(c) E 0 0 R(c) R(c) ] ;
 
 % if any( [ Pars Yic ] < 0 )
 %     error('Negative parameters or initial values.')

main_shared_noN_calibration_2023_01_01.m

@@ -0,0 +1,154 @@
+%% main.m = calls model analysis functions stored in 'Analyzer' folder
+
+addpath Analyzer
+addpath Data
+
+%% ========================================================================
+%	NOTES
+%  ========================================================================
+
+% Script performs MSLS parameter estimation and Bayesian MCMC uncertainty
+% quantification for an ODE model of N-803 treatment of SIV.
+
+% Training data is from two SIV cohorts:
+% Cohort 1 (low VL) = Ellis et. al. J Virol. 2018. (N803_DataSet_1.mat)
+% Cohort 2 (high VL) = Ellis et. al. J Virol. 2022. (N803_DataSet_2.mat)
+
+% This is a repeat of 'main_shared_calibration_2022_08_18.m' without
+% non-SIV-specific CD8+ T cells
+
+% Model is run for both cohorts (from two fixed points), with parameters
+% shared between cohorts, fitting to both cohorts simultaneously.
+% Parameter and IC calculation is done in 'N803_shared_noN.m'.
+% ODEs are solved in 'N803_model_2S.m'.
+% Files within 'Analyzer' folder perform MSLS and MCMC.
+
+%% ========================================================================
+%	INPUTS for 'calibrator.m'
+%  ========================================================================
+
+% X and Y = cell arrays for data x-values and y-values --------------------
+load('N803_DataSet_1') ; DS1 = DataSet ;
+load('N803_DataSet_2') ; DS2 = DataSet ;
+X = {DS1(1).X_data,DS1(3).X_data,DS2(1).X_data,DS2(3).X_data} ;
+Y = {DS1(1).Y_data,DS1(3).Y_data,DS2(1).Y_data,DS2(3).Y_data} ;
+% Apply censoring and normalization to data (see 'prepper' function)
+[Xp,Yp,Sid,Cid,IVs] = ...
+        prepper(X,Y, {2,[],2,[]}, {2 ,[],2 ,[]}, {0,0,0,0}, [1 0 1 0]) ;
+[X,Y] = prepper(X,Y, {2,[],2,[]}, {[],[],[],[]}, {0,0,0,0}, [1 0 1 0]) ;
+
+% GuessBOUNDS = array for boundaries of initial parameter guesses -----
+GuessBOUNDS(:,01) = [1 1]*mean(cat(1,IVs{1}{:})) ;% V iv [log(#/mL)] (co 1)
+GuessBOUNDS(:,02) = [1 1]*mean(cat(1,IVs{3}{:})) ;% V iv [log(#/mL)] (co 2)
+GuessBOUNDS(:,03) = [1 1]*mean(cat(1,IVs{2}{:})) ;% E iv [#/µL] (co 1)
+GuessBOUNDS(:,04) = [1 1]*mean(cat(1,IVs{4}{:})) ;% E iv [#/µL] (co 2)
+
+GuessBOUNDS(:,05) = [.01 1]    ;% V growth rate (if E were absent) [/d]
+GuessBOUNDS(:,06) = [1e4 1e6]  ;% 50% viral stim saturation for E [#/mL]
+GuessBOUNDS(:,07) = [1e4 1e6]  ;% 50% viral stim saturation for R [#/mL]
+GuessBOUNDS(:,08) = [.005 .05] ;% EA reversion rate constant []
+
+GuessBOUNDS(:,09) = [300 3e3]  ;% 50% E proliferation saturation [#/µL]
+GuessBOUNDS(:,10) = [0.1 1]    ;% E0 normalized prolif rate constant [/d]
+GuessBOUNDS(:,11) = [1 4]      ;% EA proliferation rate constant [/d]
+GuessBOUNDS(:,12) = [.01 .05]  ;% E0 death rate constant [/d]
+GuessBOUNDS(:,13) = [.1 .5]    ;% EA death rate constant [/d]
+
+GuessBOUNDS(:,14) = [1 1]*880  ;% X initial condition [pmol/kg]
+GuessBOUNDS(:,15) = [1 1]*0.8  ;% N803 absorption rate constant [/d]
+GuessBOUNDS(:,16) = [1 1]*2.1  ;% N803 elimination rate constant [/d]
+GuessBOUNDS(:,17) = [1 1]*1.3  ;% N803 'vol of dist'/'bioavail' [L/kg]
+GuessBOUNDS(:,18) = [1 1e3]    ;% 50% N803 stimulation concentration [pM]
+GuessBOUNDS(:,19) = [.1 2]     ;% E0 maximum prolif rate []
+GuessBOUNDS(:,20) = [0.1 1e3]  ;% E activation stimulation factor []
+GuessBOUNDS(:,21) = [0.01 100] ;% R generation stimulation factor []
+
+GuessBOUNDS(:,22) = [0.05 1]   ;% R decay rate constant [/d]
+GuessBOUNDS(:,23) = [0.01 100] ;% E0 proliferation regulation factor []
+
+% Inputs for model function and cost function -----------------------------
+Xeval = min(cat(1,X{:})) : 1 : max(cat(1,X{:}))  ;% x-values (days)
+Doses = { 7*[0:3 5:8 37:40] ... x-values for doses (Cohort 1)
+        , 7*(0:2:4)         };% x-values for doses (Cohort 2)
+Skips = { [-7 max(cat(1,X{1:2}))] ... [min max] time to solve (Cohort 1)
+        , [-7 max(cat(1,X{3:4}))] };% [min max] time to solve (Cohort 2)
+
+% modelFunction = handle for calling model function -------------------
+modelFunction = @(t,P) N803_shared_noN(t,Doses,P,Skips,[],...
+'warnOff',1,'AbsTol',1e-3,'RelTol',1e-2,'odeSolver','23s','timeOut',1)  ;
+
+% FitBOUNDS = array for boundaries of parameter fitting -------------------
+FitBOUNDS = GuessBOUNDS ;
+
+% costFunction = handle for calling cost function -------------------------
+costFunction = @(P) costFun_NLL(P,modelFunction,X,Y,'Xevaluate',Xeval) ;
+
+% costFunction(GuessBOUNDS(1,:)) ;% uncomment to debug your costFunction
+
+% filePrefix = string to prefix to file names -----------------------------
+filePrefix = 'N803_shared_noN' ;
+
+% timeLimit = number of minutes allowed for a single optimization run -----
+timeLimit = 30 ;
+
+% numberGuesses = number of initial parameter guesses ---------------------
+numberGuesses = 1000 ;
+
+% numberWorkers = number of parallel processes (enter [] or 0 for serial)
+numberWorkers = 12 ;
+
+%% ------------------------------------------------------------------------
+%   Optional inputs for results summary figures ('calibrator_plot.m')
+%  ------------------------------------------------------------------------
+
+plotSetup.Xdata = Xp ;
+plotSetup.Ydata = Yp ;
+plotSetup.CensorID  = Cid ;
+plotSetup.Fits2Plot = 1:82 ;
+plotSetup.modelFunction = modelFunction ;
+plotSetup.nFits2Sum = 82 ;
+plotSetup.SubjectID = Sid ;
+plotSetup.XLabels   = {{'Days Post-Treatment'} ...
+                      ,{'Days Post-Treatment'} ...
+                      ,{'Days Post-Treatment'} ...
+                      ,{'Days Post-Treatment'}} ;
+plotSetup.YLabels   = {{'vRNA copy','log fold \Delta'} ...
+                      ,{'CD8^{+} T cell','fold \Delta'}...
+                      ,{'vRNA copy','log fold \Delta'} ...
+                      ,{'CD8^{+} T cell','fold \Delta'}} ;
+plotSetup.XRANGE = 7*[ -4 41 ; -4 41 ; -4 6 ; -4 6 ] ;
+plotSetup.YRANGE = [ -3 1 ; 0 6 ; -2 2 ; 0 6 ] ;
+plotSetup.PlotGrid = [2 2] ;
+plotSetup.PlotOuts = [1 3 2 4] ;
+plotSetup.plotPars = 0 ;
+plotSetup.xPrecise = 1 ;
+
+%% ========================================================================
+%	PERFORM MSLS
+%  ========================================================================
+
+% Perform calibration using 'calibrator.m'
+setup = struct('GuessBOUNDS',GuessBOUNDS,'FitBOUNDS',FitBOUNDS,...
+               'costFunction',costFunction,'timeLimit',timeLimit) ;
+
+Results = calibrator(setup,filePrefix,numberGuesses,numberWorkers) ;
+
+% Plot results using 'calibrator_plot.m'
+calibrator_plot(Results,plotSetup,filePrefix)
+
+%% ========================================================================
+%	PERFORM MCMC
+%  ========================================================================
+
+% InitialRanks = ranks (by NLL) of parameter sets to initialize MCMC ------
+InitialRanks = 1:25:101 ;
+
+% MCMC sample size --------------------------------------------------------
+numSamples = 1e6 ;
+
+% save every 'saveEach' iterations ----------------------------------------
+saveEach = 1000 ;
+
+% Perform uncertainty quantification using 'MCMCsampler.m'
+ResultsFile = [filePrefix '_Results'] ;
+MCMCsampler(ResultsFile,InitialRanks,numSamples,saveEach) ;