(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

_{max}= λ

_{min}+ η.

The definitions and units of all state variables and parameters are reported respectively in Tables 1 and 2.

Table 1

Definition of state variables

State variables: symbol and meaning |
Units |
---|---|

t, time |
[mo] |

Β, β-cell mass as Millions of active β-cells |
[Mc] |

I, fasting serum insulin concentration (monthly average) |
[pM] |

G, fasting plasma glucose concentration (monthly average) |
[mM] |

h, glucose effect on pancreas |
[#] |

x, replaces Glucose variations with variations of a scale-free pure number |
[#] |

λ, net rate constant for β-cell growth (or decay), resulting from the difference between production (replication) rate and mortality (apoptosis) rate |
[mo ^{−1}] |

η, possible lambda excursion, represents pancreatic replication ‘reserve’ |
[mo ^{−1}] |

K _{xgI}, second-order insulin-dependent glucose tissue uptake rate per pM insulin, represents insulin sensitivity |
[min ^{−1}/pM] |

γ, resistance to insulin as the number of mM of glucose at which glycemia stabilizes for a single pM of insulinemia |
[mM * pM] |

λ _{max}, maximum (positive) value of λ, i.e. the maximum possible difference between replication rate and apoptosis rate |
[mo ^{−1}] |

I _{maxB}, insulin attainable levels, expressed as the maximal contribution of a million β-cells to fasting insulin plasma concentration |
[pM/Mc] |

A, glycosylated Haemoglobin, HbA1c (monthly average) |
[%] |

Table 2

Definition of model parameters

Model parameters: symbol and meaning |
Units |
---|---|

B _{0} = B(t_{0}): initial condition on β-cell mass |
[Mc] |

B _{max}: maximal size of β-cell mass |
[Mc] |

I _{0}: insulinema at age t_{0} |
[pM] |

G _{0}: glycemia at age t_{0} |
[mM] |

G _{λ}: the glycemia of maximal sensitivity of the regulation of β-cell population dynamics |
[mM] |

T _{gl}: liver glucose output. This parameter can be considered as a constant basal, insulin-independent glucose output, while insulin-dependent glucose output (decreasing with increasing levels of insulin) is taken into account through the action of KxgI |
[mM/min] |

K _{xg} (FAST): first-order insulin-independent glucose tissue uptake rate |
[min ^{−1}] |

T _{igB} (FAST): maximal pancreatic insulin secretion per million β-cells |
[pM/min/Mc] |

K _{xi} (FAST): apparent first-order elimination rate constant for insulin |
[min ^{−1}] |

K _{xiStart}: apparent first-order elimination rate constant for insulin at baseline (e.g. at age 18 years) |
[min ^{−1}] |

K _{xiEnd}: Apparent first-order elimination rate constant for insulin at the end of a normal life (e.g. at age 90 years) |
[min ^{−1}] |

K _{xgI0}: value of K_{xgI} at t_{0} |
[min ^{−1}/pM] |

t _{I}: time of midpoint Κ_{xgI} decrease (time at 0.5 * K_{xgI0}) |
[mo] |

ν _{I} : exponent for the decrease in K_{xgI} |
[#] |

ρ: ratio of 1st order to 2nd order (Insulin Dependent) rate constants for tissue glucose uptake from plasma |
[min ^{−1}/(min^{−1} /pM)] = pM |

I _{maxB0}: value of I_{maxB} at t_{0} |
[pM/Mc] |

λ _{min}: the minimum (negative) value of λ, i.e. the maximum net apoptosis rate |
[mo ^{−1}] |

η _{0}: value of η at t_{0} (determined) |
[mo ^{−1}] |

T _{η}: the spontaneous recovery rate of the pancreas |
[mo ^{−2}] |

K _{ηG}: Rate constant expressing pancreatic glucose toxicity |
[mo ^{−1}mM^{−1}] |

G _{h}: centering glucose concentration for Hill-shaped glycemia effect on pancreatic insulin release, set to 9 mM |
[mM] |

ν _{h}: power coefficient for Hill-shaped glycemia effect on pancreatic insulin release, set to 4 in order to achieve linearity between 5 and 10 mM, essentially zero before 3.5 mM, essentially saturated above 25 mM |
[#] |

h _{0}: The value of the h(G) function at t_{0} (determined) |
[#] |

ξ: the scaling factor for the X support variable, computed so as to center the X sigmoid curve on the value λ = 0 |
[#] |

A _{0} = A(t_{0}): initial condition on glycosylated Haemoglobin (HbA1c). |
[%] |

K _{ag} : the rate constant of production of glycosylated haemoglobin from circulating glucose |
[%/mo/mM] |

K _{xa}: the spontaneous elimination rate constant of (glycosylated) Haemoglobin |
[mo ^{−1}] |

t _{0}: starting age, with system at equilibrium, in months |
[mo] |

t _{η10}, Half-life of Pancreatic reserve at 10 mM Glucose |
[mo] |

Equations 3 and 4 simply state that prevailing glycemias at any historical moment are inversely related to the prevailing insulinemias, while prevailing insulinemias are directly related to (an increasing function of) prevailing glycemias as well as being directly related to available β-cell mass and secretory capacity. The “prevailing” glycemia and insulinemia levels described in the model do not strictly refer to either fasting or post-prandial values, but are in fact representative of daily steady-state blood concentrations, hourly oscillations being averaged out when considering the evolution of the system over a time span of months or years. Prevailing glycemia and insulinemia are derived as equilibrium values, given the existing level of insulin sensitivity and the current functional β-cell mass available for insulin secretion: they substantially reflect real-life post-absorptive values.

The model hypothesizes, in a straightforward population dynamics fashion, that the variation in β-cell mass B (Eq. 1) is proportional to current mass in logistic growth fashion, with carrying capacity B

_{max}, multiplied by a net-growth rate λ.In turn, λ depends (Eq. 6) on prevailing glucose concentrations, in the sense that it varies from a minimum negative value (λ

_{min}) to a maximum value that is dependent both on the prevailing pancreatic reserve described below (i.e. λ_{max}= λ_{min}+ η) and on glucose level according to a sigmoidal 3rd-degree Hill function, with λ = λ_{min}when G = 0, and λ tending to λ_{max}as G tends to infinity.The Hill function is normalized (it is a function of x rather than G) so as to have maximum slope at the glycemia value G

_{λ}. If G_{λ}is approximately equal to normal fasting glycemia (e.g. around 5 mM), this is the same as assuming that the regulation mechanisms of the pancreas work optimally in a neighborhood of the target fasting glycemia.The term η represents the ‘pancreatic reserve’, i.e. the current ability of the pancreas to increase its β-cell proliferation rate if sufficiently stimulated by the ambient glucose concentration. Pancreatic reserve, however, is not a fixed quantity: it varies with time (Eq. 2), in such a way that it always tends to an equilibrium level. This equilibrium is the result of two competing forces: glucose-toxicity-induced shrinkage of pancreatic reserve and spontaneous recovery of the pancreas. Pancreatic reserve (at equilibrium) is thus inversely proportional to prevailing glycemias. In other words, Eq. 6 prescribes that hyperglycemia will always acutely stimulate pancreatic β-cell replication (to the extent that this is allowed by the current pancreatic reserve level), while Eq. 2 prescribes that, chronically, sustained hyperglycemia will make pancreatic reserve decrease. Notice that, if pancreatic reserve is sufficiently damaged, i.e. if η < −λ

_{min}, then λ_{min}+ η < 0 and β-cell population would decline, in the absence of therapy, no matter what the prevailing glucose concentration. It should be kept in mind that this effect could well stem both from an anti-proliferative and from a pro-apoptotic original lesion.Equations 3 and 4 derive, as hinted to above, from the equilibrium condition for glycemia and insulinemia respectively, given that the model focuses on changes over months or years. At each point in slow time, therefore, the fast dynamics of glucose and insulin are assumed to be at equilibrium. This equilibrium may be conceived as being the average glycemia produced by the average insulinemias and by the current levels of insulin sensitivity. Conversely, average insulinemia is determined by the prevailing glycemias, by the β-cell mass available at the time, and by a coefficient of insulin production at maximal stimulation per β-cell mass unit. It should be pointed out that using fasting or average glucose and insulin concentrations amounts conceptually to the same thing, as we only need indicator quantities for the prevailing state of the system over a period of a day or a week.

The function h (Eq. 5) captures the hypothesis that glycemia stimulates insulin production by the pancreas in a sigmoidal fashion, starting at zero, sloping up approximately linearly between 5 and 10 mM glycemia, and approaching asymptotically a maximum as glycemia increases.

We found it necessary to stipulate that, as observed in the literature [30], the apparent first order rate constant of elimination of insulin from plasma decreases with age (Eq. 7): this allows the model to replicate the observed decline of insulin secretory function, associated with preserved or even increased levels of serum insulin, with advancing age [31].

The variable γ (inversely related to the index of insulin sensitivity K

_{xgI}) expresses resistance to insulin, as the concentration at which glucose stabilizes for each pM of insulin concentration.For the purposes of the present simulations, it has been assumed that the insulin-independent glucose tissue uptake rate constant K

_{xg}is small (within the range of considered glycemias) compared to insulin-dependent glucose uptake [19], and that consequently the ratio ρ is also approximately zero.The variable I

_{maxB}is the maximal contribution of a million β-cells to fasting plasma insulin concentration. Its value is determined by the maximal insulin secretion rate per million β-cells T_{igB}(determined by equilibrium conditions at t_{0}), and by the actual first-order apparent rate of elimination of insulin from plasma, K_{xi}, which is considered here to decrease with age (linearly, in first approximation), as mentioned above. It should be kept in mind that T_{igB}in our simulations is made to decrease when supposing primary defects of insulin secretory function with advancing age, and conversely it is made to increase in response to secretagogues.In order to link observable glycated hemoglobin with the glucose dynamics, a simple linear model of the kinetics of HbA1c (indicated as A in Eq. 8) is part of DPM, with increase determined by prevailing glycemias and by the concentration of native HbA0, and decrease linearly determined by the continuous destruction of red blood cells.

### 2.1 Parameter Assessment

Parameters for the model have been carefully assessed from literature data and from the results of directed experiments, with the goal of maintaining consistency among the several sources. Only human data have been used. Although limited information was available for some parameters, most were based on empirical observations derived using well-established methodologies. The majority of parameters were obtained from in vivo studies but in a few cases only in vitro or post-mortem data were available. Preference was generally given to studies with larger numbers of subjects, but study quality and design were also considered (e.g. age range of study population, level of validation of key measurements). For each parameter, we sought to identify a typical value, reflecting a non-diseased population, keeping in mind that the range of reasonable values for each parameter would include those seen in diabetic individuals. Typical values and ranges at time = t

_{0}were generally taken from data obtained in young, healthy adults (age 18–30).Parameter values used in this simulation were as reported in [29], with the following specifications. Moderate secretion and sensitivity defects were assumed to have a minimum reachable value of the corresponding coefficients at 15 % normal, whereas severe defects were supposed to cause a decrease of the coefficient down to 2.5 % normal. In both cases the age at 50 % defect was assumed to be 25 years. The threshold of glycemia for the diagnosis of diabetes was 7 mM. Proportional improvement by classical sensitizers or secretagogues was assumed to be 50 % with respect to the current “natural” value as determined by the progression of disease. Definitive improvement of insulin sensitivity produced by bariatric surgery was assumed to 7.5 % of normal (i.e. a small fraction of supposedly perfect insulin sensitivity as seen in young age) above the current “natural” value due to age and disease. Definitive decrease of glucose toxicity produced by β-cell protectors was assumed to be 20 % of normal.

The model’s numerical integration, starting with given parameter values, has been implemented in a mixed Matlab (© 1994–2007 The MathWorks, Inc.) and C/C++ (GCC, © Free Software Foundation, Inc.) environment, and is freely available as a service to academic users through the CNR IASI BioMatLab website [32].

## 3 Results

### 3.1 Qualitative Analysis of the Behavior of Solutions

The original published model [24] was proven to admit three equilibrium points: a physiological steady-state, which is locally asymptotically stable, consisting of the (normal) basal (i.e. ‘early-life’) values of glycemia and insulinemia G

_{b}and I_{b}; a pathological ‘severe diabetes’ steady-state, which is also locally asymptotically stable, consisting of zero levels of β-cell mass and insulinemia, to which there corresponds a hyperglycemic state (G_{h}> G_{b}); and an unstable saddle point, with intermediate glycemia (G_{b}< G_{s}< G_{h}).