From 9b04ce8c47ff1045ff498f4b0768fbb9be97906a Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Sun, 1 Jan 2023 13:16:37 -0600 Subject: [PATCH] Add files via upload --- N803_shared_noN.m | 24 ++-- main_shared_noN_calibration_2023_01_01.m | 154 +++++++++++++++++++++++ 2 files changed, 166 insertions(+), 12 deletions(-) create mode 100644 main_shared_noN_calibration_2023_01_01.m diff --git a/N803_shared_noN.m b/N803_shared_noN.m index d30f462..0bd23bd 100644 --- a/N803_shared_noN.m +++ b/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,26 +101,25 @@ 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 ; @@ -128,7 +127,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.') diff --git a/main_shared_noN_calibration_2023_01_01.m b/main_shared_noN_calibration_2023_01_01.m new file mode 100644 index 0000000..901d2ce --- /dev/null +++ b/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 [#/無] (co 1) +GuessBOUNDS(:,04) = [1 1]*mean(cat(1,IVs{4}{:})) ;% E iv [#/無] (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 [#/無] +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) ;