From 82675f04f509580cc47900b8287f7b99759a0665 Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Mon, 15 May 2023 10:23:14 -0500 Subject: [PATCH] Delete main_singles_calibration_2022_09_25.m --- main_singles_calibration_2022_09_25.m | 165 -------------------------- 1 file changed, 165 deletions(-) delete mode 100644 main_singles_calibration_2022_09_25.m diff --git a/main_singles_calibration_2022_09_25.m b/main_singles_calibration_2022_09_25.m deleted file mode 100644 index e74bb6e..0000000 --- a/main_singles_calibration_2022_09_25.m +++ /dev/null @@ -1,165 +0,0 @@ -%% 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) - -% Model is run for both cohorts (from two fixed points), with NO parameters -% shared between cohorts, fitting to both cohorts seperately. -% Parameter and IC calculation is done in 'N803_single.m'. -% ODEs are solved in 'N803_model_2.m'. -% Files within 'Analyzer' folder perform MSLS and MCMC. - -%% ======================================================================== -% INPUTS for 'calibrator.m' -% ======================================================================== - -% GuessBOUNDS = array for boundaries of initial parameter guesses --------- -GuessBOUNDS(:,03) = [.03 .3] ;% initial frequency: S/(S+N) -GuessBOUNDS(:,04) = [.009 .9] ;% initial frequency: Sa/S -GuessBOUNDS(:,05) = [.01 1] ;% V growth rate (if S+N were absent) [/d] -GuessBOUNDS(:,06) = [.01 1] ;% Na/Sa killing rate ratio [gN0/gS0] -GuessBOUNDS(:,07) = [1e4 1e6] ;% 50% viral stim saturation for S [#/mL] -GuessBOUNDS(:,08) = [1e4 1e6] ;% 50% viral stim saturation for N [#/mL] -GuessBOUNDS(:,09) = [1 10] ;% normalized Sa reversion rate constant [] -GuessBOUNDS(:,10) = [1 10] ;% normalized Na reversion rate constant [] -GuessBOUNDS(:,11) = [300 3e3] ;% 50% S+N proliferation saturation [#/uL] -GuessBOUNDS(:,12) = [1 4] ;% Sa proliferation rate constant [/d] -GuessBOUNDS(:,13) = [0.5 2] ;% Na proliferation rate constant [/d] -GuessBOUNDS(:,14) = [.01 .05] ;% S0/N0 death rate constant [/d] -GuessBOUNDS(:,15) = [.1 .5] ;% Sa/Na death rate constant [/d] -GuessBOUNDS(:,16) = [1 1]*880 ;% X initial condition [pmol/kg] -GuessBOUNDS(:,17) = [1 1]*0.8 ;% N803 absorption rate constant [/d] -GuessBOUNDS(:,18) = [1 1]*2.1 ;% N803 elimination rate constant [/d] -GuessBOUNDS(:,19) = [1 1]*1.3 ;% N803 'vol of dist'/'bioavail' [L/kg] -GuessBOUNDS(:,20) = [1 1e3] ;% 50% N803 stimulation concentration [pM] -GuessBOUNDS(:,21) = [.1 2] ;% S0/N0 maximum proliferation rate [] -GuessBOUNDS(:,22) = [0.1 1e3] ;% S activation stimulation factor [] -GuessBOUNDS(:,23) = [0.01 100] ;% N activation stimulation factor [] -GuessBOUNDS(:,24) = [0.1 100] ;% sN/sS regulation generation rate ratio -GuessBOUNDS(:,25) = [0.05 1] ;% R decay rate constant [/d] -GuessBOUNDS(:,26) = [0.01 100] ;% initial S killing regulation [] -GuessBOUNDS(:,27) = [0.01 100] ;% initial N killing regulation [] -GuessBOUNDS(:,28) = [0.01 100] ;% initial S0/N0 proliferation regulation [] -GuessBOUNDS(:,29) = [0.01 100] ;% initial S activation regulation [] -GuessBOUNDS(:,30) = [0.01 100] ;% initial N activation regulation [] - -% timeLimit = number of minutes allowed for a single optimization run ----- -timeLimit = 30 ; - -% numberGuesses = number of initial parameter guesses --------------------- -numberGuesses = 500 ; - -% numberWorkers = number of parallel processes (enter [] or 0 for serial) -numberWorkers = 0 ; - -%% ------------------------------------------------------------------------ -% Do for each cohort... -for c = 1:2 - -% X and Y = cell arrays for data x-values and y-values -------------------- -load(['N803_DataSet_' num2str(c)]) -X = {DataSet(1).X_data,DataSet(3).X_data} ; -Y = {DataSet(1).Y_data,DataSet(3).Y_data} ; -% Apply censoring and normalization to data (see 'prepper' function) -[Xp,Yp,Sid,Cid,IVs] = prepper(X,Y,{2,[]},{ 2,[]},{0,0},[1 0]) ; -[X,Y] = prepper(X,Y,{2,[]},{[],[]},{0,0},[1 0]) ; - -% GuessBOUNDS = array for boundaries of initial parameter guesses --------- -GuessBOUNDS(:,01) = [1 1]*mean(cat(1,IVs{1}{:})) ;% V iv [log(#/mL)] -GuessBOUNDS(:,02) = [1 1]*mean(cat(1,IVs{2}{:})) ;% S0+Sa iv [#/uL] - -% FitBOUNDS = array for boundaries of parameter fitting ------------------- -FitBOUNDS = GuessBOUNDS ; - -% 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 inf] ;% [min max] time point beyond which to skip solving - -% modelFunction = handle for calling model function ----------------------- -modelFunction = @(t,P) N803_single(t,Doses{c},P,Skips,... - 'warnOff',1,'AbsTol',1e-3,'RelTol',1e-2,'odeSolver','23s') ; - -% costFunction = handle for calling cost function ------------------------- -costFunction = @(P) costFun_NLL(P,modelFunction,X,Y,'Xevaluate',Xeval) ; - -% tic -% costFunction(GuessBOUNDS(1,:)) ;% uncomment to debug your costFunction -% disp(toc) - -% filePrefix = string to prefix to file names ----------------------------- -filePrefix = ['N803_cohort_' num2str(c)]; - -%% ------------------------------------------------------------------------ -% Optional inputs for results summary figures ('calibrator_plot.m') -% ------------------------------------------------------------------------ - -plotSetup.Xdata = Xp ; -plotSetup.Ydata = Yp ; -plotSetup.modelFunction = modelFunction ; -plotSetup.CensorID = Cid ; -plotSetup.Fits2Plot = 1:101 ; -plotSetup.nFits2Sum = 101 ; -plotSetup.SubjectID = Sid ; -plotSetup.XLabels = {{'Days Post-Treatment'} ... - ,{'Days Post-Treatment'}} ; -plotSetup.YLabels = {{'vRNA copy','log fold \Delta'} ... - ,{'CD8^{+} T cell','fold \Delta'}} ; -switch c - case 1 - plotSetup.XRANGE = 7*[ -4 41 ; -4 41 ] ; - plotSetup.YRANGE = [ -5 3 ; 0 6 ] ; - case 2 - plotSetup.XRANGE = 7*[ -4 6 ; -4 6] ; - plotSetup.YRANGE = [ -4 4 ; 0 6 ] ; -end -plotSetup.PlotGrid = [2 1] ; -plotSetup.PlotOuts = 1:2 ; -plotSetup.plotPars = 0 ; - -%% ======================================================================== -% FUNCTION CALLS -% ======================================================================== - -% 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) - -end - -%% ======================================================================== -% 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 ; - -parfor c = 1:2 - - % Perform uncertainty quantification using 'MCMCsampler.m' - ResultsFile = ['N803_cohort_' num2str(c) '_Results'] ; - MCMCsampler(ResultsFile,InitialRanks,numSamples,saveEach) ; -end \ No newline at end of file