diff --git a/main_shared_calibration_2022_08_18.m b/main_shared_calibration_2022_08_18.m deleted file mode 100644 index 526c7ee..0000000 --- a/main_shared_calibration_2022_08_18.m +++ /dev/null @@ -1,160 +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 parameters -% shared between cohorts, fitting to both cohorts simultaneously. -% Parameter and IC calculation is done in 'N803_shared.m'. -% ODEs are solved in 'N803_model_2.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}{:})) ;% S+N iv [#/uL] (co 1) -GuessBOUNDS(:,04) = [1 1]*mean(cat(1,IVs{4}{:})) ;% S+N iv [#/uL] (co 2) -GuessBOUNDS(:,05) = [.03 .3] ;% initial frequency: S/(S+N) (co 1) -GuessBOUNDS(:,06) = [.01 1] ;% normalized S/(S+N) (co 2) - -GuessBOUNDS(:,07) = [.01 1] ;% V growth rate (if S+N were absent) [/d] -GuessBOUNDS(:,08) = [.01 1] ;% Na/Sa killing rate ratio [gN0/gS0] -GuessBOUNDS(:,09) = [1e4 1e6] ;% 50% viral stim saturation for S [#/mL] -GuessBOUNDS(:,10) = [1e4 1e6] ;% 50% viral stim saturation for N [#/mL] -GuessBOUNDS(:,11) = [1 10] ;% normalized Sa reversion rate constant [] -GuessBOUNDS(:,12) = [1 10] ;% normalized Na reversion rate constant [] - -GuessBOUNDS(:,13) = [300 3e3] ;% 50% S+N proliferation saturation [#/uL] -GuessBOUNDS(:,14) = [0.1 1] ;% S0/N0 normalized prolif rate constant [/d] -GuessBOUNDS(:,15) = [1 4] ;% Sa proliferation rate constant [/d] -GuessBOUNDS(:,16) = [0.5 2] ;% Na proliferation rate constant [/d] -GuessBOUNDS(:,17) = [.01 .05] ;% S0/N0 death rate constant [/d] -GuessBOUNDS(:,18) = [.1 .5] ;% Sa/Na death rate constant [/d] - -GuessBOUNDS(:,19) = [1 1]*880 ;% X initial condition [pmol/kg] -GuessBOUNDS(:,20) = [1 1]*0.8 ;% N803 absorption rate constant [/d] -GuessBOUNDS(:,21) = [1 1]*2.1 ;% N803 elimination rate constant [/d] -GuessBOUNDS(:,22) = [1 1]*1.3 ;% N803 'vol of dist'/'bioavail' [L/kg] -GuessBOUNDS(:,23) = [1 1e3] ;% 50% N803 stimulation concentration [pM] -GuessBOUNDS(:,24) = [.1 2] ;% S0/N0 maximum prolif rate [] -GuessBOUNDS(:,25) = [0.1 1e3] ;% S activation stimulation factor [] -GuessBOUNDS(:,26) = [0.01 100] ;% N activation stimulation factor [] - -GuessBOUNDS(:,27) = [0.05 1] ;% R decay rate constant [/d] -GuessBOUNDS(:,28) = [0.1 100] ;% sN/sS regulation generation ratio -GuessBOUNDS(:,29) = [0.01 100] ;% S0/N0 prolif reg factor [] -GuessBOUNDS(:,30) = [0.01 100] ;% N 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(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' ; - -% 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 = 0 ; - -%% ------------------------------------------------------------------------ -% 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) ;