ReplayBG is a digital twin-based methodology to assess new strategies for type 1 diabetes management.
Installation of ReplayBG can be easily performed by including the ReplayBG folder in the MATLAB® PATH. To do so:
-
Download the latest ReplayBG release (e.g.
X.Y.Z
) from Github --> here; -
Unzip the downloaded archive and move the
replay-bg-X.Y.Z
folder to the desired location, for example~/MATLAB/
; -
Open MATLAB®;
-
To use the ReplayBG functions in the current MATLAB® session (or within a specific script), add
replay-bg-X.Y.Z
source folder and subfolders to the current MATLAB® PATH by executing in the command window (or at the top of the script)
addpath(genpath('~/MATLAB/replay-bg-X.Y.Z/src'))
- Enjoy!
The first step consists of fitting the ReplayBG model agaist data to capture the underlying physiological dynamics and estimate its parameter distributions via Markov Chain Monte Carlo (MCMC). This is possible by providing to the tool a portion of data recorded from people with type 1 diabetes consisting of glucose recordings, carbohydrate intakes, and bolus + basal insulin infusions.
ReplayBG has different level of flexibility and can perform step 1 differently according to the user preferences. The easiest way to do this is calling the replayBG
function as:
replayBG('identification', data, BW, scenario, saveName);
where:
'identification'
specifies that the tool will be used to perform step 1, i.e., to identify the ReplayBG model on the given datadata
is a timetable which contains the data to be used by the tool. MUST contain a columnglucose
that contains the glucose measurements (in mg/dl), a columnbasal
that contains the basal insulin data (in U/min), a columnbolus
that contains the bolus insulin data (in U/min), a columnCHO
that contains the CHO intake data (in g/min).data
MUST be sampled on a homogeneous time grid and MUST, except for theglucose
column, not contain Nan values. In case ofscenario = 'multi-meal'
data
(see "Optional parameters" section below) MUST also contain a column of stringschoLabel
that contains for each non-zero value of theCHO
column, a character that specifies the type of CHO intake ('B'
for breakfast,'L'
for lunch,'D'
for dinner,'S'
for snack,'H'
for hypotreatment);BW
is the patient body weight (kg);scenario
: a vector of characters that specifies whether the given scenario refers to a single-meal scenario or a multi-meal scenario. Can be'single-meal'
or'multi-meal'
;saveName
is a vector of characters used to label, thus identify, each output file and result of ReplayBG.
ReplayBG accepts optional parameters that can be used to fit user's preferences. These include:
glucoseModel
: (optional, default:'IG'
) a vector of characters that specifies the glucose model to use. Can be'IG'
or'BG'
;cgmModel
: (optional, default:'IG'
) a vector of characters that specifies the glucose model to use as cgm measurement. Can be'CGM'
,'IG'
or'BG'
;pathology
: (optional, default:'t1d'
) a vector of characters that specifies the patient pathology. Can be't1d'
,'t2d'
,'pbh'
,'healthy'
. Note that't2d'
,'pbh'
, and'healthy'
functionalities will be matter of future work and do not work at the moment;scenario
: (optional, default:'single-meal'
) a vector of characters that specifies whether the given scenario refers to a single-meal scenario or a multi-meal scenario. Can be'single-meal'
or'multi-meal'
;sampleTime
: (optional, default:5
(min)) an integer that specifies thedata
sample time;seed
: (optional, default:randi([1 1048576])
) an integer that specifies the random seed. For reproducibility;maxETAPerMCMCRun
: (optional, default:inf
) a number that specifies the maximum time in hours allowed for each MCMC run;maxMCMCIterations
: (optional, default:inf
) an integer that specifies the maximum number of iterations for each MCMC run;maxMCMCRuns
: (optional, default:inf
) an integer that specifies the maximum number of MCMC runs;maxMCMCRunsWithMaxETA
: (optional, default:inf
) an integer that specifies the maximum number of MCMC runs having maximum ETA;adaptiveSCMH
: (optional, default:1
) a numerical flag that specifies whether to make the Single Components Metropolis Hastings algorithm adaptive or non-adaptive. Can be0
or1
.MCMCTheta0Policy
: (optional, default:'mean'
) a vector of characters defining the policy used by the MCMC procedure to set the initial MCMC chain values. Can be'mean'
or'last'
or'initial'
. Using'mean'
, the mean value of the MCMC chain obtained from the last MCMC run will be set as initial MCMC chain value to be used in the next MCMC run. Using'last'
, the last value of the MCMC chain obtained from the last MCMC run will be set as initial MCMC chain value to be used in the next MCMC run. Using'initial'
, the same initial value will be used for every run of MCMC;bayesianEstimator
: (optional, default:'mean'
) a vector of characters defining which Bayesian estimator to use to obtain a point estimate of model parameters. Can be'mean'
or'map'
. Using'mean'
the posterior mean estimater will be used. Using'map'
, the marginalized maximum-a-posteriori estimator will be used. Note that this point estimate is not used during step 2 (it is just an output of the identification procedure);preFilterData
: (optional, default:0
) a numerical flag that specifies whether to filter the glucose data before performing the model identification or not. Can be0
or1
. This might help the identification procedure. Filtering is performed using a non-causal fourth-order Butterworth filter having0.1*sampleTime
cut-off frequency;saveChains
: (optional, default:1
) a numerical flag that specifies whether to save the resulting mcmc chains in dedicated files (one for each MCMC run) for future analysis or not. Can be0
or1
;saveSuffix
: (optional, default:''
) a vector of char to be attached as suffix to the resulting output files' name;plotMode
: (optional, default:1
) a numerical flag that specifies whether to show the plot of the results or not. Can be0
or1
;enableLog
: (optional, default:1
) a numerical flag that specifies whether to log the output of ReplayBG not. Can be0
or1
;verbose
: (optional, default:1
) a numerical flag that specifies the verbosity of ReplayBG. Can be0
or1
.
Potentially, given the computational cost and the function of the parameter identification machinery (MCMC), this step can take forever.
These are some suggested usage of ReaplyBG depending on your needs.
Please remember that the rationale is to identify a dataset "once and for all" (that is the heavy part of ReplayBG) and then play with the simulations (see "Step 2: Use of ReplayBG for simulation").
replayBG( ...
'identification', data, BW, 'single-meal', saveName, ...
'maxETAPerMCMCRun', 6, ...
'maxMCMCRuns', 4, ...
'plotMode',1, ...
'verbose',1, ...
'seed', 1, ...
);
This will ensure a maximum running time of 24 hours (i.e., 4 x 6 hours). Consider increasing maxETAPerMCMCRun
if warnings on MCMC convergence are raised by ReplayBG.
If you prefer to run ReplayBG "silently" or in a server infrastructure, set plotMode
and verbose
to 0
.
replayBG( ...
'identification', data, BW, 'multi-meal', saveName, ...
'maxETAPerMCMCRun', 12, ...
'maxMCMCRuns', 4, ...
'plotMode',1, ...
'verbose',1, ...
'seed', 1, ...
);
This will ensure a maximum running time of 48 hours (i.e., 4 x 12 hours). Consider increasing maxETAPerMCMCRun
if warnings on MCMC convergence are raised by ReplayBG.
If you prefer to run ReplayBG "silently" or in a server infrastructure, set plotMode
and verbose
to 0
.
replayBG( ...
'identification', data, BW, scenario, saveName, ...
'maxETAPerMCMCRun', 5/60, ...
'maxMCMCRuns', 2, ...
'plotMode',1, ...
'verbose',1, ...
'seed', 1, ...
);
This will ensure a maximum running time of 10 minutes (i.e., 2 x 5 minutes). This is useful to have a grasp of ReplayBG functioning. Since 5 minutes MCMC runs can lead to very short MCMC chains, it can happen (almost for sure) that some warnings will be raised by ReplayBG alerting that the MCMC did not converge.
The second step of ReplayBG consists of using the model parameters identified during step 1 and run simulations on the same time window but with different, altered model inputs (meal and insulin) to predict "what would have happened", in terms of glycemic time course, if such inputs would have been used instead of the orginal recorded ones. This basically allows to test the potential impact on glycemia of new insulin/carbohydrate algorithms for the management of type 1 diabetes.
As in step 1, ReplayBG has different level of flexibility and can perform step 2 differently according to the user preferences. The easiest way to do this is calling the replayBG
function as:
replayBG('replay', data, BW, scenario, saveName);
where:
'replay'
specifies that the tool will be used to simulate the scenario specified by the givendata
;data
is a timetable which contains the data to be used by the tool. MUST contain a columnglucose
that contains the glucose measurements (in mg/dl), a columnbasal
that contains the basal insulin data (in U/min), a columnbolus
that contains the bolus insulin data (in U/min), a columnCHO
that contains the CHO intake data (in g/min).data
MUST be sampled on a homogeneous time grid and MUST, except for theglucose
column, not contain Nan values. In case ofscenario = 'multi-meal'
data
(see "Optional parameters" section below) MUST also contain a column of stringschoLabel
that contains for each non-zero value of theCHO
column, a character that specifies the type of CHO intake ('B'
for breakfast,'L'
for lunch,'D'
for dinner,'S'
for snack,'H'
for hypotreatment);BW
is the patient body weight (kg);scenario
: a vector of characters that specifies whether the given scenario refers to a single-meal scenario or a multi-meal scenario. Can be'single-meal'
or'multi-meal'
;saveName
is a vector of characters used to label, thus identify, each output file and result of ReplayBG. To use specific model parameters identified in step 1, you MUST be the same correspondingsaveName
used during step 1.
ReplayBG accepts optional parameters that can be used to fit user's preferences. These include:
glucoseModel
: (optional, default:'IG'
) a vector of characters that specifies the glucose model to use. Can be'IG'
or'BG'
;cgmModel
: (optional, default:'IG'
) a vector of characters that specifies the glucose model to use as cgm measurement. Can be'CGM'
,'IG'
or'BG'
;pathology
: (optional, default:'t1d'
) a vector of characters that specifies the patient pathology. Can be't1d'
,'t2d'
,'pbh'
,'healthy'
. Note that't2d'
,'pbh'
, and'healthy'
functionalities will be matter of future work and do not work at the moment;sampleTime
: (optional, default:5
(min)) an integer that specifies thedata
sample time;bolusSource
: (optional, default:'data'
) a vector of character defining whether to use, during replay, the insulin bolus data contained in thedata
timetable (ifdata
), or the boluses generated by the bolus calculator implemented via the providedbolusCalculatorHandler
function. Can be'data'
or'dss'
. It cannot be set ifmodality
is'identification'
;basalSource
: (optional, default:'data'
) a vector of character defining whether to use, during replay, the insulin basal data contained in thedata
timetable (ifdata
), or the basal generated by the controller implemented via the providedbasalControllerHandler
function (ifdss
), or fixed to the average basal rate used during identification (if'u2ss'
). Can be'data'
,'u2ss'
, or'dss'
. It cannot be set ifmodality
is'identification'
;choSource
: (optional, default:'data'
) a vector of character defining whether to use, during replay, the CHO data contained in thedata
timetable (ifdata
), or the CHO generated by the meal generator implemented via the providedmealGeneratorHandler
function. Can be'data'
or'generated'
. It cannot be set ifmodality
is'identification'
;CR
: (optional, default: 10) the carbohydrate-to-insulin ratio of the patient in g/U to be used by the integrated decision support system;CF
: (optional, default: 40) the correction factor of the patient in mg/dl/U to be used by the integrated decision support system;GT
: (optional, default: 120) the target glucose value in mg/dl to be used by the decsion support system modules;bolusCalculatorHandler
: (optional, default:'standardBolusCalculatorHandler'
) a vector of characters that specifies the name of the function handler that implements a bolus calculator to be used during the replay of a given scenario whenbolusSource
is'dss'
. The function must have 2 output, i.e., the computed insulin bolus (U/min) and a structure containing the dss hyperparameter. The function must have 7 inputs, i.e.,G
(mg/dl) a vector as long the simulation length containing all the simulated glucose concentrations up totimeIndex
(the other values are nan),mealAnnouncements
(g/min) a vector that contains the announced meal CHO intakes inputs for the whole replay simulation,bolus
(U/min) a vector that contains the bolus insulin input for the whole replay simulation,basal
(U/min) a vector that contains the basal insulin input for the whole replay simulation,time
(datetime) a vector that contains the time istants of current replay simulation,timeIndex
is a number that defines the current time istant in the replay simulation,dss
a structure containing the dss hyperparameters and the optionally providedbolusCalculatorHandlerParams
(dss
is also echoed in the output to enable memory-like features). Vectors contain one value for each integration step. The default bolus calculator implemented bystandardBolusCalculatorHandler
is the standard formula: B = CHO/CR + (GC-GT)/CF - IOB;bolusCalculatorHandlerParams
: (optional, default:[]
) a structure that contains the parameters to pass to the bolusCalculatorHandler function. It also serves as memory area for thebolusCalculatorHandler
function;basalHandler
: (optional, default:'defaultBasalHandler'
) a vector of characters that specifies the name of the function handler that implements a basal controller to be used during the replay of a given scenario whenbasalSource
is'dss'
. The function must have 2 output, i.e., the computed basal rate (U/min) and a structure containing the dss hyperparameter. The function must have 7 inputs, i.e.,G
(mg/dl) a vector as long the simulation length containing all the simulated glucose concentrations up totimeIndex
(the other values are nan),mealAnnouncements
(g/min) a vector that contains the announced meal CHO intakes inputs for the whole replay simulation,bolus
(U/min) a vector that contains the bolus insulin input for the whole replay simulation,basal
(U/min) a vector that contains the basal insulin input for the whole replay simulation,time
(datetime) a vector that contains the time istants of current replay simulation,timeIndex
is a number that defines the current time istant in the replay simulation,dss
a structure containing the dss hyperparameters and the optionally providedbasalHandlerParams
(dss
is also echoed in the output to enable memory-like features). Vectors contain one value for each integration step. The default basal controller implemented by'defaultBasalHandler'
is: if G < 70, basal = 0, otherwise basal = basal(1).basalHandlerParams
: (optional, default:[]
) a structure that contains the parameters to pass to the basalHandler function. It also serves as memory area for thebasalHandler
function;mealGeneratorHandler
: (optional, default:'defaultMealGeneratorHandler'
) a vector of characters that specifies the name of the function handler that implements a meal generator to be used during the replay of a given scenario whenchoSource
is'generated'
. The function must have 4 outputs, i.e., the generated meal (g/min) actually consumed by the virtual subject, the generated meal announcement (g/min) that is usually used to compute the corresponding insulin bolus, the type of the meal (must be'B'
,'L'
,'D'
, or'S'
), and a structure containing the dss hyperparameter. The function must have 8 inputs, i.e.,G
(mg/dl) a vector as long the simulation length containing all the simulated glucose concentrations up totimeIndex
(the other values are nan),meal
(g/min) a vector that contains the meal CHO intakes inputs for the whole replay simulation,mealAnnouncements
(g/min) a vector that contains the announced meal CHO intakes inputs for the whole replay simulation,bolus
(U/min) a vector that contains the bolus insulin input for the whole replay simulation,basal
(U/min) a vector that contains the basal insulin input for the whole replay simulation,time
(datetime) a vector that contains the time istants of current replay simulation,timeIndex
is a number that defines the current time istant in the replay simulation,dss
a structure containing the dss hyperparameters and the optionally providedmealGeneratorHandlerParams
(dss
is also echoed in the output to enable memory-like features). Vectors contain one value for each integration step. The default meal generator implemented bydefaultMealGeneratorHandler
is: put a snack meal of 50g of CHO in the first instant and announce only 40g.mealGeneratorHandlerParams
: (optional, default:[]
) a structure that contains the parameters to pass to themealGeneratorHandler
function. It also serves as memory area for themealGeneratorHandler
function;enableHypoTreatments
: (optional, default:0
) a numerical flag that specifies whether to enable hypotreatments during the replay of a given scenario. Can be0
or1
. Can be set only when ReplayBG is used in'replay'
mode;hypoTreatmentsHandler
: (optional, default:'adaHypoTreatmentsHandler'
) a vector of characters that specifies the name of the function handler that implements an hypotreatment strategy during the replay of a given scenario. The function must have 2 outputs, i.e., the hypotreatments carbohydrates intake (g/min) and thedss
structure (see section "Use case 2: Test an hypotreatment strategy with ReplayBG" below). The function must have 8 inputs, i.e.,G
(mg/dl) a glucose vector as long the simulation length containing all the simulated glucose concentrations up to timeIndex (the other values are nan),CHO
(g/min) a vector that contains the CHO intakes input for the whole replay simulation,hypotreatments
(g/min) a vector that contains the hypotreatments intakes input for the whole replay simulation (will be also populated by ReplayBG during the simulation if one or more hypotreatment is given by the function itself),bolus
(U/min) a vector that contains the bolus insulin input for the whole replay simulation,basal
(U/min) a vector that contains the basal insulin input for the whole replay simulation,time
(datetime) a vector that contains the time instants of current replay simulation,timeIndex
is a number that defines the current time instant in the replay simulation and thedss
structure. Vectors contain one value for each integration step. The default policy implemented by theadaHypoTreatmentsHandler
function is "take an hypotreatment of 10 g every 15 minutes while in hypoglycemia";enableCorrectionBoluses
: (optional, default:0
) a numerical flag that specifies whether to enable correction boluses during the replay of a given scenario. Can be0
or1
. Can be set only when ReplayBG is used in'replay'
mode;correctionBolusesHandler
: (optional, default:'correctsAbove250Handler'
) a vector of characters that specifies the name of the function handler that implements a corrective bolusing strategy during the replay of a given scenario. The function must have 2 output, i.e., the correction insulin bolus (U/min) and thedss
structure (see section "Use case 3: Test a corrective bolusing strategy with ReplayBG" below). The function must have 7 inputs, i.e.,G
(mg/dl) a glucose vector as long the simulation length containing all the simulated glucose concentrations up to timeIndex (the other values are nan),CHO
(g/min) a vector that contains the CHO intakes input for the whole replay simulation,bolus
(U/min) a vector that contains the bolus insulin input for the whole replay simulation,basal
(U/min) a vector that contains the basal insulin input for the whole replay simulation,time
(datetime) a vector that contains the time instants of current replay simulation,timeIndex
is a number that defines the current time instant in the replay simulation and thedss
structure. Vectors contain one value for each integration step. The default policy implemented by thecorrectsAbove250Handler
function is "take a corrective bolus of 1 U every 1 hour while above 250 mg/dl";hypoTreatmentsHandlerParams
: (optional, default:[]
) a structure that contains the parameters to pass to thehypoTreatmentsHandler
function (see section "Use case 2: Test an hypotreatment strategy with ReplayBG" below). It also serves as memory area for thehypoTreatmentsHandler
function;correctionBolusesHandlerParams
: (optional, default:[]
) a structure that contains the parameters to pass to thecorrectionBolusesHandler
function (see section "Use case 3: Test a corrective bolusing strategy with ReplayBG" below). It also serves as memory area for thecorrectionBolusesHandler
function;seed
: (optional, default:randi([1 1048576])
) an integer that specifies the random seed. For reproducibility;saveSuffix
: (optional, default:''
) a vector of char to be attached as suffix to the resulting output files' name;plotMode
: (optional, default:1
) a numerical flag that specifies whether to show the plot of the results or not. Can be0
or1
;enableLog
: (optional, default:1
) a numerical flag that specifies whether to log the output of ReplayBG not. Can be0
or1
;verbose
: (optional, default:1
) a numerical flag that specifies the verbosity of ReplayBG. Can be0
or1
.
Potentially, given the computational cost and the function of the parameter identification machinery (MCMC), this step can take forever.
These are some suggested usage of ReaplyBG depending on your needs.
Please remember that the rationale is to identify a dataset "once and for all" (that is the heavy part of ReplayBG) and then play with the simulations (see "Step 2: Use of ReplayBG for simulation").
The most simple use case of ReaplyBG consists of modifying the data
timetable used during step 1 setting new insulin/meal inputs and replaying the scenario. For example if you want to simulate "what would have happened if I had doubled the CHO intake?":
%Here, I suppose that data is the timetable used during identification
data.CHO = data.CHO*2;
replayBG( ...
'replay', data, BW, scenario, saveName, ...
'plotMode',1, ...
'verbose',1, ...
'seed', 1, ...
);
The same rationale can be applied to the insulin bolus/basal input.
Note that if you want to add a meal and scenario
is multi-meal
you also have to set the corresponding choLabel
, e.g.:
%Here, I suppose that data is the timetable used during identification
data.CHO(12) = 5; %Here I am adding a snack of 5 g/min at 12*`sampleTime` minutes after the starting of the simulation.
data.choLabel(12) = 'S';
replayBG( ...
'replay', data, BW, scenario, saveName, ...
'plotMode',1, ...
'verbose',1, ...
'seed', 1, ...
);
Let's say that you want to test an algorithm that advice to take an hypotreatment according to some strategy to avoid hypoglycemia.
You have to follow two steps:
- First define a function that implements your strategy, namely the "handler". This function, as stated above, must have 2 outputs, i.e., the hypotreatments carbohydrates intake (g/min) and the
dss
structure. The function must have 8 inputs, i.e.,G
(mg/dl) a glucose vector as long the simulation length containing all the simulated glucose concentrations up to timeIndex (the other values are nan),CHO
(g/min) a vector that contains the CHO intakes input for the whole replay simulation,hypotreatments
(g/min) a vector that contains the hypotreatments intakes input for the whole replay simulation (will be also populated by ReplayBG during the simulation if one or more hypotreatment is given by the function itself),bolus
(U/min) a vector that contains the bolus insulin input for the whole replay simulation,basal
(U/min) a vector that contains the basal insulin input for the whole replay simulation,time
(datetime) a vector that contains the time instants of current replay simulation,timeIndex
is a number that defines the current time instant in the replay simulation and thedss
structure. For example, here's the default handler implemented in ReplayBG that provides the user with an example to start playing with:
function [HT, dss] = adaHypoTreatmentsHandler(G,CHO,hypotreatments,bolus,basal,time,timeIndex,dss)
% function adaHypoTreatmentsHandler(G,CHO,hypotreatments,bolus,basal,time,timeIndex,dss)
% Implements the default hypotreatment strategy: "take an hypotreatment of
% 10 g every 15 minutes while in hypoglycemia".
%
% ---------------------------------------------------------------------
%
% Copyright (C) 2020 Giacomo Cappon
%
% This file is part of ReplayBG.
%
% ---------------------------------------------------------------------
HT = 0;
%If glucose is lower than 70...
if(G(timeIndex) < 70)
%...and if there are no CHO intakes in the last 15 minutes, then take an HT
if(timeIndex > 15 && ~any(hypotreatments((timeIndex - 15):timeIndex)))
HT = 15; % g/min
end
end
end
- Call ReplayBG by enabling the hypotreatments "module" and by providing the name of the handler just created:
%Here, I suppose that data is the timetable used during identification
replayBG( ...
'replay', data, BW, scenario, saveName, ...
'enableHypoTreatments', 1, ...
'hypoTreatmentsHandler', 'adaHypoTreatmentsHandler', ...
'plotMode',1, ...
'verbose',1, ...
'seed', 1, ...
);
Three notes:
- The handler function will be called for each integration time step (each "simulated" minute). This will slow down ReplayBG.
- If
scenario
is'single meal'
, thehypotreatments
input will contain only the hypotreatments generated by this function during the simulation. Ifscenario
is'multi-meal'
, thehypotreatments
input will ALSO contain the hypotreatments already present in the given data that labeled as such. CHO
does not contain hypotreatments.dss
is a structure that serves as memory area and contains decsion support related parameter passed to ReplayBG when you call thereplayBG
function (i.e., the parametersBW,
GT,
CR,
CF,
bolusCalculatorHandlerParams,
basalHandlerParams,
mealGeneratorHandlerParams,
hypoTreatmentsHandlerParams, and
correctionBolusesHandlerParams). As such you can use it to provide to your handler the parameter you need. Also, keep in mind that, serving as a memory area, it is possible to store values inside
dssand
hypoTreatmentsHandlerParams` and they will be available in the next call of the function (this is useful if you need to compute things like the insulin-on-board, or "when is the last time I gave an hypotreatment?").
Let's say that you want to test an algorithm that advice to take a corrective insulin bolus according to some strategy to avoid hyperglycemia.
You have to follow two steps:
- First define a function that implements your strategy, namely the "handler". This function, as stated above, must have 2 outputs, i.e., the hypotreatments carbohydrates intake (g/min) and the
dss
structure. The function must have 7 inputs, i.e.,G
(mg/dl) a glucose vector as long the simulation length containing all the simulated glucose concentrations up to timeIndex (the other values are nan),CHO
(g/min) a vector that contains the CHO intakes input for the whole replay simulation,bolus
(U/min) a vector that contains the bolus insulin input for the whole replay simulation,basal
(U/min) a vector that contains the basal insulin input for the whole replay simulation,time
(datetime) a vector that contains the time instants of current replay simulation,timeIndex
is a number that defines the current time instant in the replay simulation and thedss
structure. For example, here's the default handler implemented in ReplayBG that provides the user with an example to start playing with:
function [CB, dss] = correctsAbove250Handler(G,CHO,bolus,basal,time,timeIndex,dss)
% function correctsAbove250Handler(G,CHO,bolus,basal,time,timeIndex,dss)
% Implements the default correction bolus strategy: "take a correction
% bolus of 1 U every 1 hour while above 250 mg/dl".
%
% ---------------------------------------------------------------------
%
% Copyright (C) 2020 Giacomo Cappon
%
% This file is part of ReplayBG.
%
% ---------------------------------------------------------------------
CB = 0;
%If glucose is greater than 250...
if(G(timeIndex) > 250)
%...and if there are no boluses in the last 1 hour, then take a CB
if(timeIndex > 60 && ~any(bolus((timeIndex - 60):timeIndex)))
CB = 1; % U/min
end
end
end
- Call ReplayBG by enabling the correction insulin "module" and by providing the name of the handler just created:
%Here, I suppose that data is the timetable used during identification
replayBG( ...
'replay', data, BW, scenario, saveName, ...
'enableCorrectionBoluses', 1, ...
'correctionBolusesHandler', 'correctsAbove250Handler', ...
'plotMode',1, ...
'verbose',1, ...
'seed', 1, ...
);
Three notes:
- The handler function will be called for each integration time step (each "simulated" minute). This will slow down ReplayBG.
dss
is a structure that serves as memory area and contains decsion support related parameter passed to ReplayBG when you call thereplayBG
function (i.e., the parametersBW,
GT,
CR,
CF,
bolusCalculatorHandlerParams,
basalHandlerParams,
mealGeneratorHandlerParams,
hypoTreatmentsHandlerParams, and
correctionBolusesHandlerParams). As such you can use it to provide to your handler the parameter you need. Also, keep in mind that, serving as a memory area, it is possible to store values inside
dssand
correctionBolusesHandlerParams` and they will be available in the next call of the function (this is useful if you need to compute things like the insulin-on-board, or "when is the last time I gave a correction bolus?").
Let's say that you want to test an algorithm that computes the meal insulin bolus according to some strategy.
You have to follow two steps:
- First define a function that implements your strategy, namely the "handler". This function, as stated above, the function must have 2 output, i.e., the computed insulin bolus (U/min) and a structure containing the dss hyperparameter. The function must have 7 inputs, i.e.,
G
(mg/dl) a vector as long the simulation length containing all the simulated glucose concentrations up totimeIndex
(the other values are nan),mealAnnouncements
(g/min) a vector that contains the announced meal CHO intakes inputs for the whole replay simulation,bolus
(U/min) a vector that contains the bolus insulin input for the whole replay simulation,basal
(U/min) a vector that contains the basal insulin input for the whole replay simulation,time
(datetime) a vector that contains the time istants of current replay simulation,timeIndex
is a number that defines the current time istant in the replay simulation,dss
a structure containing the dss hyperparameters. For example, here's the default handler implemented in ReplayBG that provides the user with an example to start playing with:
function [B, dss] = standardBolusCalculatorHandler(G, mealAnnouncements,bolus,basal,time,timeIndex,dss)
% function standardBolusCalculatorHandler(G,mealAnnouncements,bolus,basal,time,timeIndex,dss)
% Implements the default insulin bolus calculator formula: B = CHO/CR + (GC -GT)/CF - IOB.
%
% ---------------------------------------------------------------------
%
% Copyright (C) 2022 Giacomo Cappon
%
% This file is part of ReplayBG.
%
% ---------------------------------------------------------------------
B = 0;
%If a meal is announced...
if(mealAnnouncements(timeIndex) > 0)
%...give a bolus
B = mealAnnouncements(timeIndex)/dss.CR + (G(timeIndex) - dss.GT) / dss.CF - iobCalculation(bolus(1:timeIndex),5);
end
end
function [IOB] = iobCalculation(insulin,Ts)
% define 6 hour curve
k1 = 0.0173;
k2 = 0.0116;
k3 = 6.75;
IOB_6h_curve = zeros(360,1);
for t = 1:360
IOB_6h_curve(t)= 1 - ...
0.75*((-k3/(k2*(k1-k2))*(exp(-k2*(t)/0.75)-1) + ...
k3/(k1*(k1-k2))*(exp(-k1*(t)/0.75)-1))/(2.4947e4));
end
IOB_6h_curve = IOB_6h_curve(Ts:Ts:end);
% IOB is the convolution of insulin data with IOB curve
IOB = conv(insulin, IOB_6h_curve);
IOB = IOB(length(insulin));
end
- Call ReplayBG by changing the
bolusSource
from'data'
to'dss'
and by providing the name of the handler just created:
%Here, I suppose that data is the timetable used during identification
replayBG( ...
'replay', data, BW, scenario, saveName, ...
'bolusSource','dss', ...
'bolusCalculatorHandler', 'standardBolusCalculator', ...
'plotMode',1, ...
'verbose',1, ...
'seed', 1, ...
);
Three notes:
- The handler function will be called for each integration time step (each "simulated" minute). This will slow down ReplayBG.
dss
is a structure that serves as memory area and contains decision support related parameter passed to ReplayBG when you call thereplayBG
function (i.e., the parametersBW,
GT,
CR,
CF,
bolusCalculatorHandlerParams,
basalHandlerParams,
mealGeneratorHandlerParams,
hypoTreatmentsHandlerParams, and
correctionBolusesHandlerParams). As such you can use it to provide to your handler the parameter you need. Also, keep in mind that, serving as a memory area, it is possible to store values inside
dssand
bolusCalculatorHandlerParams` and they will be available in the next call of the function (this is useful if you need to compute things like the insulin-on-board, or "when is the last time I gave an insulin bolus?").
Let's say that you want to test an algorithm that basal rate infusion according to some strategy.
You have to follow two steps:
- First define a function that implements your strategy, namely the "handler". This function, as stated above, the function must have 2 output, i.e., the computed basal rate (U/min) and a structure containing the dss hyperparameter. The function must have 7 inputs, i.e.,
G
(mg/dl) a vector as long the simulation length containing all the simulated glucose concentrations up totimeIndex
(the other values are nan),mealAnnouncements
(g/min) a vector that contains the announced meal CHO intakes inputs for the whole replay simulation,bolus
(U/min) a vector that contains the bolus insulin input for the whole replay simulation,basal
(U/min) a vector that contains the basal insulin input for the whole replay simulation,time
(datetime) a vector that contains the time istants of current replay simulation,timeIndex
is a number that defines the current time istant in the replay simulation,dss
a structure containing the dss hyperparameters. For example, here's the default handler implemented in ReplayBG that provides the user with an example to start playing with:
function [B, dss] = defaultBasalHandler(G, mealAnnouncements,bolus,basal,time,timeIndex,dss)
% function defaultBasalHandler(G,mealAnnouncements,bolus,basal,time,timeIndex,dss)
% Implements the default basal rate controller: if G < 70, basal = 0,
% otherwise basal = basal(1).
%
% ---------------------------------------------------------------------
%
% Copyright (C) 2022 Giacomo Cappon
%
% This file is part of ReplayBG.
%
% ---------------------------------------------------------------------
B = basal(1);
%If G < 70...
if(G(timeIndex) < 70)
%...set basal rate to 0.
B = 0;
end
end
- Call ReplayBG by changing the
basalSource
from'data'
to'dss'
and by providing the name of the handler just created:
%Here, I suppose that data is the timetable used during identification
replayBG( ...
'replay', data, BW, scenario, saveName, ...
'basalSource','dss', ...
'basalHandler', 'defaultBasalHandler', ...
'plotMode',1, ...
'verbose',1, ...
'seed', 1, ...
);
Three notes:
- The handler function will be called for each integration time step (each "simulated" minute). This will slow down ReplayBG.
dss
is a structure that serves as memory area and contains decision support related parameter passed to ReplayBG when you call thereplayBG
function (i.e., the parametersBW
,GT
,CR
,CF
,bolusCalculatorHandlerParams
,basalHandlerParams
,mealGeneratorHandlerParams
,hypoTreatmentsHandlerParams
, andcorrectionBolusesHandlerParams
). As such you can use it to provide to your handler the parameter you need. Also, keep in mind that, serving as a memory area, it is possible to store values insidedss
andbolusCalculatorHandlerParams
and they will be available in the next call of the function (this is useful if you need to compute things like the insulin-on-board, or "for how much did I suspend basal infusion?").
Let's say that you want to test an algorithm that generates meals according to some strategy.
You have to follow two steps:
- First define a function that implements your strategy, namely the "handler". This function, as stated above, the function must have 4 outputs, i.e., the generated meal (g/min) actually consumed by the virtual subject, the generated meal announcement (g/min) that is usually used to compute the corresponding insulin bolus, the type of the meal (must be
'B'
,'L'
,'D'
, or'S'
), and a structure containing the dss hyperparameter. The function must have 8 inputs, i.e.,G
(mg/dl) a vector as long the simulation length containing all the simulated glucose concentrations up totimeIndex
(the other values are nan),meal
(g/min) a vector that contains the meal CHO intakes inputs for the whole replay simulation,mealAnnouncements
(g/min) a vector that contains the announced meal CHO intakes inputs for the whole replay simulation,bolus
(U/min) a vector that contains the bolus insulin input for the whole replay simulation,basal
(U/min) a vector that contains the basal insulin input for the whole replay simulation,time
(datetime) a vector that contains the time istants of current replay simulation,timeIndex
is a number that defines the current time istant in the replay simulation,dss
a structure containing the dss hyperparameters. For example, here's the default handler implemented in ReplayBG that provides the user with an example to start playing with:
function [C, MA, type, dss] = defaultMealGeneratorHandler(G, meal, mealAnnouncements,bolus,basal,time,timeIndex,dss)
% function defaultMealGeneatorHandler(G,meal,mealAnnouncements,bolus,basal,time,timeIndex,dss)
% Implements the default meal generation policy: put a snack meal of 50g of CHO
% in the first instant and announce only 40g.
%
% ---------------------------------------------------------------------
%
% Copyright (C) 2022 Giacomo Cappon
%
% This file is part of ReplayBG.
%
% ---------------------------------------------------------------------
%Default output values
C = 0;
MA = 0;
type = '';
%If this is the first time instant...
if(timeIndex == 1)
%...generate a snack meal of 50g and announce just 40g.
C = 50;
MA = 40;
type = 'S';
end
end
- Call ReplayBG by changing the
mealSource
from'data'
to'generated'
and by providing the name of the handler just created:
%Here, I suppose that data is the timetable used during identification
replayBG( ...
'replay', data, BW, scenario, saveName, ...
'mealSource','generated', ...
'mealGeneratorHandler', 'defaultMealGeneratorHandler', ...
'plotMode',1, ...
'verbose',1, ...
'seed', 1, ...
);
Three notes:
- The handler function will be called for each integration time step (each "simulated" minute). This will slow down ReplayBG.
dss
is a structure that serves as memory area and contains decision support related parameter passed to ReplayBG when you call thereplayBG
function (i.e., the parametersBW,
GT,
CR,
CF,
bolusCalculatorHandlerParams,
basalHandlerParams,
mealGeneratorHandlerParams,
hypoTreatmentsHandlerParams, and
correctionBolusesHandlerParams). As such you can use it to provide to your handler the parameter you need. Also, keep in mind that, serving as a memory area, it is possible to store values inside
dssand
bolusCalculatorHandlerParams` and they will be available in the next call of the function (this is useful if you need to compute things like the insulin-on-board, or "when did I eat last time?").
Results are saved in the results/
folder of the replay-bg folder, specifically:
results/distributions/
: contains the identified ReplayBG model parameter distributions obtained via MCMC;results/logs/
: contains .txt files that log the command window output of ReplayBG. NB: .txt files will be empty if verbose = 0;results/mcmcChains/
: contains the MCMC chains, for each unknown parameter, obtained in each MCMC run;results/modelParameters/
: contains the model parameters identified using MCMC. Known model parameters are fixed to population values obtained from the literature. Unknown model parameters estimates are in draws..samples;results/workspaces/
: contains the core ReplayBG variables and data used in a specific ReplayBG call plus a structure calledanalysis
that contains useful metrics evaluated on the simulated glucose trace (this uses the AGATA toolbox). Name of the files start with'identification'
or'replay'
depending on how you called thereplayBG
function.
- Only the 'single-meal' mode has been extensively validated and evaluated in the referred paper. The 'multi-meal' model is working correctly but it is currently under development.
- The code is continuosly evolving and if you want to support its development please feel free to contact me at [email protected].