From aeec59d1640088f4262acd703b9cd253f2aabbf0 Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Sun, 1 Jan 2023 13:16:23 -0600 Subject: [PATCH] Delete main_shared_noN_calibration_2022_11_02.m --- main_shared_noN_calibration_2022_11_02.m | 157 ----------------------- 1 file changed, 157 deletions(-) delete mode 100644 main_shared_noN_calibration_2022_11_02.m diff --git a/main_shared_noN_calibration_2022_11_02.m b/main_shared_noN_calibration_2022_11_02.m deleted file mode 100644 index 293cee6..0000000 --- a/main_shared_noN_calibration_2022_11_02.m +++ /dev/null @@ -1,157 +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) - -% 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 [] -GuessBOUNDS(:,24) = [0.01 100] ;% EA killing 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') ; - -% 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) ; - -% tic -% costFunction(GuessBOUNDS(1,:)) ;% uncomment to debug your costFunction -% disp(toc) - -% 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 = 2000 ; - -% 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:9:82 ; - -% 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) ;