**3. Results**

### *3.1. OGTT Modeling*

In this approach, we first propose a model, which is adapted from the model MINMOD [21], to reproduce the OGTT results obtained in [17,18] from the protocol described in Figure 1. We refer to [23] for a review of glucose regulation models, and, together with [24–27], for a modeling of the OGTT in a more complex and exhaustive manner. Finally, we highlight the work [28] which contains a very detailed model of glucose response after a meal. The oral minimal model defined in [29] assembles three "minimal" models in order to build a quantitative minimal model of both glucose and insulin evolution in the context of OGTT. However, this model is still too complex for our simple purpose. The trade off of not using a model as detailed as [29] is that our conclusion will only be qualitative: the fitted models and parameters cannot be used for quantitative predictive purposes unlike those in [21,29].

The MINMOD model from [21] is not designed for OGTT, but for intravenous glucose tolerance tests (IvGTT). Therefore, we cannot consider that the plasma glucose is already at its maximal concentration at t = 0, as it is done for IvGTT studies.

Complex OGTT models such as [28] use compartmental modeling to represent the multiple stages of glucose distribution in the body, and to obtain the glucose rate of appearance in plasma after the meal. As a first approximation, we propose a simpler modeling using direct experimental results measuring the glucose rate of appearance after ingestion in rat, as follows. Following Wielinga et al. [30], we set the maximum of the rate of glucose appearance in the rat circulation ∼30 min after the meal (we consider that the time food spends in stomach is close to zero for the glucose solution). Similarly, the initial value of the rate of appearance was ∼70% of its maximum (see Figure 2). With this approximation, we eluded the steep increase of the glucose rate of appearance in the first few minutes of the experiments. This is justified as we focus on the simulation of the plasma glucose concentration with the first data point 10 min after glucose feeding. However, this implies imprecision with respect to the numerical variation of insulin between the fasting period and the peak of the insulin production. The curve of the glucose rate of appearance as a function of time was modeled by the continuous function GRA(*t*):

$$\mathbf{G}\_{\rm RA}(t) = K \frac{1}{\sigma \sqrt{2\pi}} \mathbf{e}^{\frac{-(t-\mu)^2}{2\sigma^2}},\tag{2}$$

where 1 *σ* √<sup>2</sup>*π e* −(*<sup>t</sup>*−*μ*)<sup>2</sup> <sup>2</sup>*σ*<sup>2</sup> is a Gaussian function centred on *μ*. Fitting to the experimental curve in ([30], Figure 4), we chose *μ* = 30 min that is the time of maximal appearance rate, and *σ* = 35 to obtain ∼70% of the maximum rate at *t* = 0. The value of the parameter *K* is determined by the actual quantity of glucose fed to the rats, such that all the glucose has to be absorbed over the [0, + ∞[ time interval. Let mGlc be the mass of glucose fed to the rats (see Table 2 for the chosen rat body mass in the simulations at PND21, -26, and -60 with respect to the measurement in [17,18]) and *VBlood* the rat blood volume (see Table 2) as given in [31]. Then, given an administrated concentration of glucose mGlc/*VBlood*, the value of *K* (Table 2) is the solution of the following equation:

$$\int\_0^{+\infty} \mathbf{K} \frac{1}{\sigma \sqrt{2\pi}} \mathbf{e}^{\frac{-(t-\mu)^2}{2\sigma^2}} \mathbf{d}t = \frac{\mathbf{m}\_{\text{Glc}}}{V\_{\text{Blood}}} . \tag{3}$$

**Figure 2.** Approximation of the plasma glucose rate of appearance in 21-day-old rats (PND21).


Upon secretagogue (here glucose) stimulation, insulin release by *β*-cells occurs in two overlapping phases [32]. The first phase, hereafter called 'fast', corresponds to the mobilization of granules belonging to the 'readily releasable pool' (RRP) by fusion to the plasma membrane and quasi-immediate excretion: it lasts for a few minutes. The second ('slow') phase lasts longer because it involves more complex phenomena such as trafficking to the surface of more deeply stored granules, and production and maturation of new insulin granules up to increased insulin transcription and synthesis to replenish the RRP. The duration of this second phase of the order of 1–2 h allows the organism to start responding to the insulin increase by the uptake of glucose in processing organs such as liver, muscles, and adipose tissues, and parallel insulin clearance by the liver. Measuring glucose clearance in the circulation over time thus integrates all these phenomena and probes the efficiency of insulin secretion by islets of Langerhans in response to glucose and the insulin sensitivity of glucose processing tissues. Of note, the first and second phases of insulin production cannot be distinguished in OGTT as delays in glucose levels in the circulation depend on intestinal absorption. The MINMOD model correctly simulates the slow phase of insulin production, but not the fast one outside of the IVGTT experimental context. To address this problem, we added a state variable representing the insulin already present and ready to be released in the blood circulation. Finally, the adapted MINMOD model is given by Equation (4).

$$\begin{aligned} \dot{\mathbf{G}} &= -p\_1(\mathbf{G}(t) - \mathbf{G}\_b) - r\_{\mathbf{C}d} \mathbf{X}(t) \mathbf{G}(t) + \mathbf{G}\_{\text{RA}}(t) \\ \dot{\mathbf{X}} &= -p\_2 \mathbf{X}(t) + p\_3(\mathbf{I}(t) - \mathbf{I}\_b) \\ \dot{\mathbf{I}} &= -n \mathbf{I}(t) + \gamma \left( \mathbf{G}(t) - h \right) t + p\_4 \mathbf{I}\_b(t) \\ \dot{\mathbf{I}}\_s &= -p\_4 \mathbf{I}\_s(t) \end{aligned} \tag{4}$$

We recall that, in the ODE system in Equation (4), G is the glucose concentration in circulating blood, X is the rate of glucose withdrawal by muscles and adipocytes due to insulin, I is the insulin concentration in circulating blood, and I*s* is the insulin concentration stored in the *β*-cells and ready to

be released. Here, we did not model the C-peptide concentration nor the liver production or absorption of glucose unlike in the more detailed model of [29].

The parameters *p*1, *p*2, *p*3, *n*, *γ*, G*b*, I*b*, and *h* are original parameters of the MINMOD model. However, they were refitted in our experimental context. We refer to Section 2.2 for a more detailed introduction to the original MINMOD model. In addition to the parameters of the original MINMOD model, we added, for our purpose of modeling OGTT and cadmium impact, three additional parameters. The flux associated with glucose intake is modeled by GRA(*t*) in Equation (2). The efficiency of the rate of glucose withdrawal X(*t*) is modeled by *rCd*. This parameter *rCd* was always taken equal to 1.0 in the simulations corresponding to the data of the control group. Finally, the parameter *p*4 models the secretion rate of the readily releasable pool of insulin modeled by the variable <sup>I</sup>*s*(*t*).

The datasets were obtained from three groups of pups [18], namely control, Cd1, and Cd2 born from female rats with background, medium, and relatively high cadmium burden, respectively, still all corresponding to low exposure doses [18]. OGTT were performed on pups at weaning (21 days after birth, that is PND21), a few days later after shifting on a regular non contaminated chow (PND26), and pcorfive weeks later (PND60). To estimate the goodness of fit of a given simulation compared to the experimental data, we used the root of the weighted least squares error:

$$\varepsilon(\mathbf{k}) = \sqrt{\sum\_{i} \mathcal{W}\_{i} (\mathbf{x}\_{\text{exp},i} - \mathbf{x}\_{\text{simu}}(t\_{i}, \mathbf{k}))^{2}} \,\,\,\tag{5}$$

where **k** is a parameter set, **<sup>x</sup>***exp*,*<sup>i</sup>* are the mean values associated to an experimental dataset (see Table 1 for its corresponding population), and **<sup>x</sup>***simu*(*<sup>t</sup>*, **k**) its associated simulation of the OGTT. The weight *Wi* is determined by the equation:

$$\mathcal{W}\_{\dot{i}} = \frac{1}{\nu\_{\dot{i}}^2 \left(\sum\_{i} \mathbf{x}\_{exp,i}^2\right)} \sim$$

where *ν*2 *i*are the variance to the mean associated to the *i*th data point.

It follows that, for a given parameter set **k**, the lower the fitting error *<sup>ε</sup>*(**k**), the better the fitting of the mean of the experimental data. When fitting the experimental results, this implies a bias in favor of the mean values of the data points as the error *ε*(**k**) is 0 if the simulation goes through all mean points. Let us remark that the fitting error will decrease when the variances *ν*2 *i* increase: an experimental point is easier to fit when its experimental uncertainty increases. This adjustment allows giving more importance to the fitting of points which are experimentally in close positions as they will be the ones which really matter in the decrease of *<sup>ε</sup>*(**k**). Here, the uncertainty on the experimental results not only comes from the precision of the measurements, but also, and mainly, from the individual variations within groups.

It is important to note that the error is *specific* to a dataset, and that errors associated to two different datasets cannot be compared. Only comparing the effect of two parameter sets to mean data points (and associated variances) corresponding to a given group (control, Cd1, or Cd2) at a given period (PND21, -26, or -60) makes sense. In Figure 3, we provide a schematic representation of our method to test the possible target of cadmium at PND21, PND26, and PND60, with respect to both our OGTT modeling and the experimental data.

**Figure 3.** Modeling method to test hypotheses on the possible targets of cadmium. In the diagram, the index *i* spans hypotheses 1.1–3.3.
