Clinically Relevant Mathematical Model for the BCG-based Treatment Of Type 1 Diabetes

This work introduces a model of immunotherapy treatment, namely the Bacillus Calmette-Guerin (BCG) vaccine, of type 1 diabetes (T1D). The model takes into consideration a clinically relevant interaction network between multiple immune cells and compartments. A set of ordinary differential equations (ODEs) is introduced to capture the connectivity between these variables and clinical presentation of the disease. Four subsets of the T1D patients and healthy controls that exhibit normal and high-level glucose consumption are evaluated. The results that obtained for mice, suggest that BCG treatment of the T1D patients that follow healthy eating habits normalizes glucose to levels observed in non-diabetic controls. Furthermore, glucose consumption profoundly influences disease progression. The stable equilibrium state with constant glucose levels is not attainable without repeated BCG treatment. This outcome suggests that immunotherapy may modulate molecular and cellular manifestations of the disease but it does not eliminate T1D. Of note, our data indicate that the BCG immunotherapy treatment may benefit healthy controls on a high-glucose diet. One may speculate the preventive BCG treatment to provide long-term health benefits in this specific cohort. Author summary We proposed a clinically relevant mathematical model of Bacillus Calmette-Guerin (BCG) based immunotherapy for type 1 diabetes (T1D) treatment. The model shows that BCG treatment is able to delay the T1D effects and to provide long-term health benefits while it may modulate molecular and cellular manifestations of the disease but it does not eliminate T1D. The main advantage of the proposed model is the ability to personalize the treatment protocol according to the patient’s metabolism and diet.


Bacillus Calmette-Guérin (BCG) vaccine is one of the oldest immunotherapies in the
clinical practice [14]. BCG treatment of the T1D patients was introduced relatively 25 recently [15]. BCG was originally developed to treat tuberculosis (TB) and early 26 (non-muscle invasive) stages of bladder cancer [16,17]. A detailed analysis of the 27 literature describes the utility of BCG to a broad spectrum of autoimmune, allergic, and 28 induced adaptive immune conditions [18,19,20]. Of particular significance to this work, 29 several instances of BCG vaccines administration in the childhood yielded a dramatic 30 improvement of the T1D progression [21]. 31 Namely, a randomized eight-year-long prospective examination of T1D patients with 32 the long-term disease who received two doses of BCG vaccine was reported by 33 Kühtreiber et al. [15]. The authors showed that BCG furnished a robust and lasting 34 control of the blood sugar levels [15]. Furthermore, the authors discovered that BCG 35 induced a systemic shift in glucose metabolism from oxidative phosphorylation to 36 aerobic glycolysis [22,23]. 37 Multiple molecular and cellular variables affect the onset of T1D. We reasoned that 38 a clinically relevant mathematical model that recaptures these parameters, as well as 39 the reported BCG-induced clinical observations, may offer a practical insight into 40 immunotherapy of the disease [24,25]. Earlier, Maree and Kublik [24] introduced an 41 algorithm that describes the interaction between β-cells, T-cells, resting macrophages, 42 active macrophages, and immune system-related cells. The authors analyzed the β-cells 43 free equilibria which affect T1D. However, the model did not take into consideration the 44 decreased glucose intake, one of the main factor associated with T1D [3]. In our view, outcomes for a mouse model. Our results and presentation are organized as follows: 58 1. Section 2 provides biological background of T1D and BCG-related treatment used 59 for the mathematical model; 60 2. Section 3 offers the analysis of the proposed model to yield a set of clinically 61 relevant equilibria states; 62 3. Section 4 deals with numerical analysis of the model using healthy subjects, T1D 63 patients and corresponding data from the NOD mouse model of T1D. It further defines 64 a baseline behavior followed by the sensitivity analysis of metabolism-and 65 treatment-related parameters including optimal treatment regimen. The suggested model is based on the clinically validated pulsed administration of BCG 71 to increase the amount of β-cells in the body to prevent T1D [27]. The model describes 72 diverse interactions between BCG and the immune system to arrive at the expected 73 outcome of the treatment on the T-cell and β-cell counts.

75
In T1D or insulin-dependent diabetes, the pancreas produces little or no insulin, a 76 hormone that facilitates glucose to enter cells in order to produce energy. The 77 autoimmune destruction of pancreatic beta cells is one of the key hallmarks of T1D.

78
Glucose metabolism is mediated by multiple types of cells that interact and are 79 controlled by intrinsic and extrinsic factors. Specific cell types include β-cells, T-cells, 80 resting and activated macrophages and dendritic cells [28]. The depleted pool of beta 81 cells triggers insulin deficiency and additional pathophysiology including life-threatening 82 hypoglycemia, ketoacidosis, cardiovascular disease, retinopathy, diabetic renal disease, 83 and neuropathy. In a healthy organism, the basic population of β-cells is supplemented 84 with additional cells as a response to the increased glucose levels after meals [29]. 85 Once the resulting pool of β-cells metabolizes available glucose, they get inactivated 86 during insulin generation [28] and their interaction with Th-cells [30]. This interaction 87 induces production of dendritic cells [28,29] that in turn prompt resting macrophages to 88 process intracellular antigenic peptides and to eliminate glucose from the circulation [29]. 89 Notably, active macrophages produce IL-1 and TNF to recruit resting macrophages.

90
Once this step is completed, active macrophages get deactivated and return to their 91 resting state [28,29]. Importantly, β-cell, T-cell, dendritic cell populations are tightly 92 controlled to exhibit exponential decay [28]. It is speculated that as BCG is distributed 93 and eliminated post-injection, it engages the dendritic cells and activated macrophages 94 [31,32] to result in the elevated β-cells concentration and overall balancing of glucose 95 metabolism. A summary of this cellular interaction network is shown on Fig. 1.

96
Considering the evidence accumulated thus far, significant attention was dedicated 97 to the design and development of immunotherapies to address T1D. Of these agents, the 98 BCG vaccine attracted particular attention primarily due to an array of immune and 99 induced adaptive immune responses exemplified by its clinical performance in multiple 100 sclerosis [33] and non-muscle invasive bladder cancer [34]. Notably, the BCG vaccine 101 was reported to activate immune cells and accelerate glycolysis [32]. Bolus injection of 102 the BCG vaccine-induced significant remission in T1D patients vs controls [31]. Furthermore, our approach accounts for cell count dynamics resulting from cell-cell 112 interactions as well as pharmacokinetics and pharmacodynamics processes. These 113 processes are captured using the following system of non-linear, ordinary differential 114 equations (ODEs). represented using a Dirac delta function [35]). ii) BCG is depleted as a result of natural 119 exponential decay with half life of µ −1 B . iii) BCG levels decrease proportionally p > 0 to 120 the interactions with effector cells including macrophages and dendritic cells [36]. Eq.
(2), the dM R (t) dt describes the amount of resting macrophages as a function of 122 time. This process is affected by several positive and negative factors. Positive factors 123 are: i) the endogenous supply or recruitment rate of macrophages a, ii) recruitment rate 124 γ of macrophages due to IL-1 and TNF-α produced by activated macrophages and iii) 125 rate of addition of macrophages by deactivation back to resting macrophages awaiting 126 activation. Negative factors include i) rate of decrease in the macrophage population c 127 due to the natural efflux of macrophages [24] and ii) rate e of macrophages activation 128 due to their interaction with dendritic cells; this factor is proportional to the antigenic 129 peptide and macrophage's population size.
reflects the amount of activated macrophages as a function of time. 131 This variable is affected by the interaction between macrophages and dendritic cells at 132 the rate e to increase the number of activated macrophages. Processing of the 133 intracellular antigenic peptides by activated macrophages affords the opposite effect on 134 the macrophage population to lead to resting macrophages at the rate k.
Eq. (4), dD(t) dt describes levels of dendritic cells as a function of time. The variable is 136 affected by i) β cell antigenic peptides released by β cells as a result of cell-cell 137 interactions with the autolytic Th-lymphocytes [37]; ii) antigenic peptides produced by 138 the BCG infected cells at the rate φ; iii) dendritic cells that are cleared at the rate µ D . 139 Eq. (5), dG(t) dt introduces the amount of glucose as a function of time. The variable 140 is affected by i) dietary glucose in the amount of g > 0 that is introduced every τ G time 141 units (represented using a Dirac delta function [35]); ii) active macrophages that eliminate glucose at the rate w; iii) glucose consumption due to the natural insulin 143 generated by the β cell and released by degrading β cells at the rate υ [38].
Eq. (6), dT (t) dt defines levels of autolytic T cells as a function of time. This variable is 145 affected by i) the natural supply or the replacement rate s T of autolytic T cells; ii) 146 proliferation rate s due to the immune cells-induced release of cytokines and chemokines 147 iii) the T cell population decrease due to lack of stimulation at the rate of µ T .
Eq. (7), dβ(t) dt describes levels of β cells as a function of time. This variable is 149 influenced by i) the natural supply or replacement rate s B of β of beta cells; ii) 150 glucose-dependent natural production of β cells at the rate ζ ; iii) β cell count decrease 151 due to β cell-Th-cell interaction [39]; iv) β cells count decrease due to the natural 152 attrition at the rate µ B ; v) β cells count decrease in a rate ξ due to generation of 153 insulin in order to process the glucose.
Considering the above, we arrive at the following system of non-linear, first-order, 155 fully-depended ODEs: Average parameter values for the model are summarized in Table. 1, for the case of 157 a mouse.

158
Equilibria and stability analysis 159 Equilibrium state presumes that the dynamic system does not change without external 160 stimuli. Stable equilibria are not affected by minor factors which can be taken out of 161 consideration when modeling complex biological systems. In our proposed model (see 162 Section ), the equilibria of interest presume the dendritic cell population is equal or 163 close to zero. This is owing to either weak immunologic response caused by the disease 164 burden or due to the weakened response by the immune system to the injected BCG. 165 Therefore, we assume D(t) = 0 which leads to an equilibrium: whereg := Σ ∞ n=0 gδ(t − nτ G ) such that qs T = µ β µ T , υµ β = qs T , and 167 ( ξg υ − s β ) 2 ≥ 4(µ β − qs T υ )(g υ ). Indeed, we obtain B * = 0 as expected. In addition, 168 M * R = a/c which is the normal rate of M A cells due to the recurrent and natural decay 169 of the M R cells population indicative of the non-compromised immune response due to 170 the treatment. Furthermore, G * =g/(υ * 1,2 ) ) which implies that the glucose levels are 171 periodic and exhibit a first-order influence on the amount of β-cells.

172
For the system of Eq. (8), we evaluated the Jacobian J: where Υ := wqs T ck and (Υaγ − µ T γ − sck) 2 ≥ 4(Υγ 2 )(Υa 2 − µ T a). This equilibrium is Therefore, As w, q, s t , γ, a, µ T > 0 by definition, the equilibrium takes the form: and due to the fact that β * = 0 the glucose levels raise to lead to the unstable 188 equilibrium as well as to unrealistic clinical scenario. As a result, the equilibrium in Eq. 189 (11) is unstable.

191
We solve the model (Eq. (8)) numerically using the ode15s function in Matlab (version 192 2020b) [43,44], where the model's initial conditions are as follows [15]:  Table 1, does not affect healthy individuals 202 (Figs. 2a-2b. Importantly, the protocol delays T1D markers (e.g., G > 110 for a day on 203 average) in T1D affected cohort on a healthy diet (Fig. 2c. The baseline treatment 204 protocol is not effective in T1D patients on high glucose diet.
Our data suggest that there is a relationship between the patient's dietary habits, 206 such as daily glucose intake (g) and genetic factors, such as metabolism and β cell 207 production due to glucose stimulus (s β ). Sensitivity analysis g ∈ [0, 50], s β ∈ [800, 1200] 208 shows the effect of both parameters on glucose levels three months after the last BCG 209 treatment (Fig. 3).
A linear approximation is chosen to balance between the accuracy vs the sampled to afford a coefficient of determination R 2 = 0.749. Therefore, the influence of g is 217 estimated to be 10 2 times larger than the influence of s β .

218
Because the eating habits and physiology are unique to each patient, the optimal 219 treatment is expected to differ between patients. In order to find the most relevant     3,4,5,6]). Green (circle) pixels represent successful treatment and red (star) pixels represent unsuccessful treatment. The parameters taken from Table 1. set of border pixels defining a surface can be approximated using the least mean square 232 (LMS) method [45]. The family function for the surface approximation is selected as 233 follows: to produce a coefficient of determination R 2 = 0.874. model. In the context of the BCG therapy, our data propose that glucose consumption 256 exert more influence (est. 10 2 difference) on the disease progression than metabolism, 257 especially when evaluated over a clinically relevant treatment time (Fig. 3).

258
Figs. 4c and 4d illustrates the profound effect of glucose on the outcome of the T1D. 259 Namely, T1D patients featured in Fig. 4d are expected to perform much worse 260 compared to patients shown on Fig. 4c even following the insulin treatment.

261
Intriguingly, data analysis summarized on Fig. 4d indicates that the BCG treatment 262 may still benefit healthy controls on a high-glucose diet. One may speculate the 263 preventive BCG treatment to provide long-term health benefits in this specific cohort. 264 Eqs. (17-20) represent the border function between the successful and unsuccessful 265 treatment protocols, differing in the amount of BCG injected, the rate and number of 266 the injections. 267 We conclude that it is possible to find the optimal treatment protocol by setting a 268 specific desired value for properties comprising the protocol to recommend a 269 personalized treatment that corresponds to the patient's glucose intake.
270 Moreover, the equilibrium described by Eq. (9) suggests a long-term system stability. 271 This outcome is clinically relevant as it indicates that the BCG treatment may have a 272 defined longitudinal effect and should be repeated as needed. Eq. (11) corroborates this 273 conclusion to imply that a biologically stable state of the system with stable glucose 274 levels is not attainable without repeating the BCG treatment. In other words, the BCG 275 treatment seems to modulate molecular and cellular manifestations of the disease but 276 does not eliminate T1D.

277
In our view, the proposed system of parameters and equations introduces both 278 accurate and practical clinical models for the immunotherapy-based treatment of T1D. 279 Our model further proposes that additional interventional modalities, for example 280 mechanistically different immunotherapy including vaccines or engineered cells may 281 enhance both targeted and symptomatic treatment outcomes.