## Abstract

Few attempts have been made to model mathematically the progression of type 2 diabetes. A realistic representation of the long-term physiological adaptation to developing insulin resistance is necessary for effectively designing clinical trials and evaluating diabetes prevention or disease modification therapies. Writing a good model for diabetes progression is difficult because the long time span of the disease makes experimental verification of modeling hypotheses extremely awkward. In this context, it is of primary importance that the assumptions underlying the model equations properly reflect established physiology and that the mathematical formulation of the model give rise only to physically plausible behavior of the solutions. In the present work, a model of the pancreatic islet compensation is formulated, its physiological assumptions are presented, some fundamental qualitative characteristics of its solutions are established, the numerical values assigned to its parameters are extensively discussed (also with reference to available cross-sectional epidemiologic data), and its performance over the span of a lifetime is simulated under various conditions, including worsening insulin resistance and primary replication defects. The differences with respect to two previously proposed models of diabetes progression are highlighted, and therefore, the model is proposed as a realistic, robust description of the evolution of the compensation of the glucose-insulin system in healthy and diabetic individuals. Model simulations can be run from the authors' web page.

- glucose
- insulin resistance
- β-cell mass
- mathematical models
- type 2 diabetes mellitus

glycemia and insulinemia are regulated in the short term (minutes to hours) through a negative feedback loop involving β-cell glucose-stimulated release of insulin from β-cells as well as insulin-mediated increased tissue glucose uptake and decreased liver gluconeogenesis and glycogenolysis (21). Also, there is evidence from rodents that chronic hyperglycemia may also induce negative feedback by increasing the net rate of β-cell replication (β-cell replication or differentiation from progenitor ductal cells − β-cell apoptosis) (9, 23, 32, 33, 59), although it has been observed that, after prolonged severe hyperglycemia, β-cell replication rates may be depressed, possibly due to glucose toxicity (66).

The interplay of insulin sensitivity, pancreatic β-cell responsiveness, and β-cell population dynamics underlies the development of type 2 diabetes mellitus (T2DM). Whatever the genetic or environmental determinants, a progressive increase of glycemia is likely determined by the progressively increasing liver and peripheral tissue insulin resistance and a corresponding failure of the pancreatic ability to compensate insulin resistance by elevating circulating insulin concentrations. In populations, this phenomenon appears to underlie the progression from the normal state to the prediabetic states of impaired glucose tolerance or impaired fasting glucose to frank T2DM. As an example, a significantly increased β-cell mass has been observed in normoglycemic subjects with obesity-induced insulin resistance (37), whereas individuals with T2DM have decreased β-cell mass relative to lean or obese nondiabetic controls (16).

Therefore, it appears that most cases of T2DM result from a combination of insulin resistance and a relative decrease in β-cell mass and/or function. In fact, isolated insulin resistance, as can be seen in severely obese individuals and in pregnant women, and reduced β-cell mass, occurring for instance in preclinical type 1 diabetes mellitus (T1DM) (24) and in pancreatectomized animals after islet transplantation (60), do not usually by themselves lead to appreciable hyperglycemia.

In the effort to understand quantitatively the interrelationships of the several players in the control of glycemia, mathematical models (5, 8, 10, 40, 44, 50, 63) have been used, typically aimed at representing the short-term evolution of glycemia and insulinemia over the time course of perturbation experiments like the intravenous glucose tolerance test, the oral glucose tolerance test, or the euglycemic hyperinsulinemic clamp. In these cases, statistical estimation of the model parameters from experimentally observed concentrations allows the snapshot determination of the degree of insulin resistance or β-cell function at a particular time in the life of the subject.

Mathematical models have also been proposed for the representation of the long-term interaction of β-cell population, fasting plasma glucose (FPG), and fasting serum insulin (FSI) in the development of T2DM (19, 61) or as discrete-state dynamical systems with covariate-dependent transition probabilities among states (22, 57, 68). In this case, however, the possibility of actually estimating the relevant parameters from observations is greatly reduced by the time scale over which the development of T2DM takes place and by the practical impossibility of serially observing the size of the β-cell population.

Much of our understanding of the progression of disease is derived from cross-sectional studies and from rational extrapolations of limited experimental results. In fact, we do not have prolonged longitudinal studies, and in this context mathematical modeling provides a precise way to coherently assemble diverse pieces of information to approximate the overall structure of the phenomenon.

The slow time scale of the process, as described by, say, FPG, FSI, or glycosylated hemoglobin (Hb A_{1c}), coupled with the lack of reliable early (surrogate) biomarkers for the development of T2DM, makes it imperative that the mathematical structure of the proposed model be thoroughly understood and verified. In this case, in fact, simple approximation to observed curves is essentially impossible, and errors in the assumptions translate into unrealistic projections beyond the relatively short observation span that is granted to studies. Grossly exemplifying, if a very naive model were to assume a linear increase of glycemia at rates observed early in the prediabetic stage, it would probably fail to predict the occurrence of T2DM in most patients who actually do develop it. In this sense, it is paramount to verify that the assumptions made when writing the model are consistent with its desired explanatory power, that the nonestimable parameters (derived from the literature or established by consensus of experts) are as precisely assessed as possible and that confidence regions are also assessed to provide a range of credible predictions, and, finally, that the mathematically derived qualitative properties of the model under study agree with commonly accepted physiology.

The goal of the present work is to offer a carefully constructed model of the evolution of glycemia, insulinemia, and β-cell mass over the several decades during which the development of T2DM typically takes place. This model will be accompanied by parameter values, deduced from the available literature, whose validity for human studies will be discussed in detail.

Although the reported model was originally independently built from first principles (its derivation is detailed in appendix 1), its structure turned out in the end to be rather similar to that of the model previously published by Topp et al. (61). The differences between the two will be discussed regarding the concept of “pancreatic reserve” introduced in the present model to explain the phenomenon of glucose toxicity, the postulated evolution of insulin clearance, and the choice of parameter values that make the model approximate commonly observed human data. Conversely, the philosophy underlying the present models and the models used by Topp et al. differs markedly from that underlying the model used by de Winter and colleagues (19, 20). Such difference and its implications will also be discussed.

Finally, the ability of the present model to represent the time course of the glucose-insulin system of a single individual over many years, in the absence or presence of significantly worsening insulin resistance or declining pancreatic reserve, will be graphically illustrated.

## METHODS

### Model Structure

The considered model's equations are (see appendix 1 for the derivation of the model and for the meaning of the suffix “slow”) with and

The definitions and units of all state variables and parameters are reported in Table 1.

As a starting point in the consideration of the model, we may notice how *Eqs. 3* and *4* simply state that prevailing glycemias at any historical moment are inversely related to the prevailing insulinemias, whereas prevailing insulinemias are directly related to (an increasing function of) prevailing glycemias as well as available β-cell mass and secretory capacity.

Clearly, at any given historical time, the levels of glycemia and insulinemia will represent an equilibrium point between these opposing tendencies. We will now consider in turn the several equations of the model, discussing the assumptions underlying the form of each.

The model hypothesizes, in a straightforward population dynamics fashion, that the variation in β-cell mass B (*Eq. 1*) equals the current mass multiplied by a net growth rate constant λ, where a negative value implies net population death.

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 on both the prevailing pancreatic reserve described below (i.e., λ_{max} = λ_{min} + η) and glucose level according to a sigmoidal third-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., ∼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 toward an equilibrium level. This equilibrium represents the steady state of two competing forces, glucose/toxicity-induced shrinkage of pancreatic reserve and spontaneous recovery of the pancreas. Pancreatic reserve 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), whereas *Eq. 2* prescribes that, chronically, sustained hyperglycemia will make pancreatic reserve decrease toward zero. 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.

*Eqs. 3*_{slow} and *4*_{slow} derive, as hinted to above, from the equilibrium condition for glycemia and insulinemia, respectively, given that the “slow” model focuses on changes over months or years. Therefore, at each point in slow time, 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 need only indicator quantities for the prevailing state of the system over a period of 1 day or 1 wk.

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 (35), the apparent first-order rate constant of elimination of insulin from plasma decreases with age (*Eq. 7*_{slow}); this allows the model to replicate the observed decline of insulin secretory function with advancing age (29).

The variable γ is the converse of glucose effectiveness [insulin sensitivity (K_{xgI})] and expresses resistance to insulin as the concentration at which glucose stabilizes for each picomolar 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 with insulin-dependent glucose uptake (50) and that, consequently, the ratio ρ is also approximately zero.

The variable I_{maxB} is the maximal contribution of 1 million β-cells to fasting plasma insulin concentration. Its value is determined by the maximal insulin secretion rate per million β-cells T_{igB} (supposed constant for the purpose of the simulations reported in the present work and 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.

To link observable glycated hemoglobin dynamics with the glucose dynamics, a simple linear model of the kinetics of Hb A_{1c} (indicated as A in *Eq. 8*_{slow}) has been hypothesized, with increase determined by prevailing glycemias and by the concentration of native Hb A_{0} and decrease linearly determined by the continuous destruction of red blood cells.

#### Parameter assessment.

An extensive (if directed) review of the literature was undertaken to identify physiological values for the model parameters defined above. 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 postmortem 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). In vivo studies utilized adult subjects, and pediatric or developmental data were not considered.

Most parameter values were taken from a single source reference (i.e., not averaged between studies), but supportive references have also been provided. For each parameter, we sought to identify a typical value, reflecting a nondiseased population, keeping in mind that the range of reasonable values for each parameter would include those seen in diabetic individuals. The parameter values are listed below with justification for the choice of values, particularly when extrapolations have been made from the referenced work or when more than one reference has been used to derive a parameter estimate. Typical values and ranges at time = t_{0} were generally taken from data obtained in young, healthy adults (age 18–30 yr); t_{0} has been defined as age 18 yr for the present work.

#### β-cell mass (B).

The largest human autopsy series reporting parameters related to β-cell mass or number is from Butler et al. (16). However, these studies reported only fractional β-cell volume and did not take into account the variable size of the pancreata. Also, the studied sample included few lean nondiabetic subjects, and these were elderly (mean 78 yr). A second autopsy series by Sakuraba et al. (56) reported a mean β-cell mass of 1.14 g. These results are consistent with other reports (52, 67). Assuming water density (1 g/ml) and a cuboidal shape with sides of 10 μm (1), the mass of a single β-cell can be estimated to be ∼1 ng. Therefore, a typical human pancreas is estimated to contain 1.14 × 10^{9} β-cells (1.14 g ÷ 10^{−9} g/cell). Other approaches yield similar estimates. For example, the human pancreas is estimated to contain 300,000–1,500,000 islets (55). Each islet contains ∼2,000 cells (58), based on 822 × 10^{6} cells in 344 × 10^{3} islet equivalents, of which ∼50% are β-cells (14). These estimates would predict ∼1 × 10^{9} β-cells in a pancreas containing 1,000,000 islets. On the other hand, the mass or number of β-cells is known to differ between individuals. β-Cells were 1.6% of pancreatic volume in lean nondiabetics (16) but ranged from 0.2 to 10% in another study (53). β-Cell mass is also known to be decreased ≥90% in long-standing type 1 diabetes mellitus. Assuming from this that β-cell number can vary as much as 50-fold, we considered reasonable a range of β-cell number from 0.1 × 10^{9} to 8 × 10^{9} cells/pancreas.

#### FSI (I) and FPG (G).

Typical values for these parameters were derived from the National Health and Nutrition Examination Survey data set (pooled years 1999–2004) (17a).

#### Net β-cell replication rate (λ).

This parameter represents the net rate constant for β-cell growth (or decay) resulting from the difference between production rate and mortality (apoptosis) rate. The production rate, in turn, reflects the sum of the rates of both replication of existing β-cells and neogenesis from non-β-cell precursors. Although current methods do not allow precise estimates of these parameters in humans, the available literature was used to define a range of reasonable limits.

Typical β-cell replication rate was estimated as 0.4% of β-cells/day on the basis of bromodeoxyuridine labeling of β-cells in human islets maintained in culture or transplanted into nude mice (62). This value was derived from the labeling index of islets from young adult donors after adjustment for the duration of replication in the same study (0.575 days). This replication rate appears consistent with data using Ki-67 immunostaining of human autopsy specimens if the same correction is applied (11, 16). We have adjusted for this duration of replication for all estimates of replication and neogenesis (since the methods used typically provide only a cross-sectional view of the percentage of cells undergoing cell division and a rate can be estimated only after adjustment for duration of the process). The range of reasonable replication rates was taken as ∼0–0.7%/day (16, 62).

Neogenesis of β-cells is believed to arise from transdifferentiation of pancreatic duct cells (12). These ductal cells are reported to comprise up to one-third of pancreatic mass and replicate at a rate of 0.05%/day (12). When extrapolating reported data to age 20 yr, ∼2% of ductal cells are insulin positive (16), so we have assumed that this is the proportion of the total duct cells that have a β-cell fate. On the basis of this, the rate of β-cell production from neogenesis can be estimated as one-third of pancreas mass represented by duct cells multiplied by 122 g as the average mass of the human pancreas (56) multiplied by 0.05% of duct cells replicating/day multiplied by 2% of duct cells differentiating into β-cells, or ∼0.37 mg β-cells/day from neogenesis. Therefore, using the typical β-cell mass of 1.14 g, neogenesis would contribute 0.03%/day to total β-cell mass, and the sum of replication and neogenesis would contribute 0.43%/day new cells to β-cell mass, with the majority deriving from replication of existing β-cells. The rate of neogenesis also appears to decline with age (16). A range of duct cell replication of 0.02–0.09% has been observed (12). This would translate to a range of neogenesis of 0.15–0.72 mg/day, which represents a relatively minor contribution to the whole β-cell production rate, ∼10%.

The frequency of apoptosis has been estimated as 2.8% of β-cells, according to data obtained from autopsy specimens in nondiabetics without pancreatic disease and using caspase-3 as a marker of apoptosis (47). Others have reported similar rates (45). Assuming a duration of apoptosis of ∼24 h, as seen in other cell types (28), these data suggest that ∼3% of β-cells undergo apoptosis each day and rates as high as 10% are estimated from data in other groups (16, 47). We believe that these rates of cell death are implausible and that the use of postmortem specimens may grossly overestimate in vivo apoptotic rates, since (if combined with the above estimates for β-cell replication) they would predict rapid β-cell failure as a near inevitability. On the other hand, the rates of β-cell replication estimated above are more plausible and are supported by independent sources of data. Therefore, we have confided in the replication rates and have assessed net β-cell replication rate (λ), the difference between rates of new β-cell growth and the rate of β-cell death, with reference to published replication rates only. Specifically, λ_{max} could conceivably be 0.7%/day using the maximum replication rate noted above and assuming no apoptosis. This would predict that β-cell mass could double in a little more than 3 mo under these hypothetical extreme conditions. A minimum net replication rate (λ_{min}) of −0.7% was arbitrarily chosen, which would symmetrically imply halving of β-cell mass over the same time period. A rate λ that would produce doubling of the β-cell mass in 1 yr would be approximately +0.2%/day, or ∼6%/mo, and we have therefore (conservatively) set λ_{max} = +6%/mo, and symmetrically, λ_{min} = −6%/mo. We expect that these crude assumptions will be tested and refined as our model is verified against observed data.

We have finally set λ_{0} equal to zero, assuming an equilibrium situation with no net change in β-cell mass in the healthy individual at age 20 yr.

β-cell replication rates are actually reported to decline gradually with age (62). We could have allowed the “typical” value of λ to decline with age with the assumption that the apoptosis rate does not change (42). However, the primary nature of this decline in β-cell net replication rate is not clear, because it may actually be secondary to accumulated glucose toxicity or be an apparent effect of a decline in β-cell mass. Both of these phenomena are accounted for by terms already included in the present model. For this reason, and following Ockham's razor criterion [“Numquam ponenda est pluralitas sine necessitate.” (Never assume many causes without need) (65)], we have preferred not to add to the current model an additional arbitrary function representing the supposed primary decrement of λ over time.

Pancreatic reserve (η) is a measure of the maximal replication rate of pancreatic β-cells, depending on the current state of pancreatic health. It has been modeled as moving from some starting value (η_{0}) and then changing [potentially increasing to a maximum (η_{max}) or decreasing toward zero] depending on the prevailing glycemia levels. Hyperglycemia is supposed to be toxic to β-cell replication (66), and hence, sustained hyperglycemia will lead to a decrease of η.

Therefore, glucose toxicity is represented as a different and separate mechanism with respect to the “instantaneous” (in slow time) proliferation capability of the β-cells. In other words, at any given level of pancreatic health (i.e., at any given value of η), the proliferation response of β-cells to increases in glycemia is always positive.

Since we had set maximal net replication rate at +6%/mo, and from this we had set maximal net apoptosis rate at −6%/mo, it follows that the maximal range η has been set at 12%/mo.

#### Glucose toxicity (K_{nG}).

This is the rate constant expressing the effect of (hyper)glycemia on the pancreatic reserve η. This parameter embodies an operational definition of glucose toxicity as the relationship between (prevailing) glucose concentration and λ. In vitro studies with human islets by Maedler et al. (43) demonstrated that β-cell replication progressively decreases and apoptosis increases as ambient glucose is increased from 5.5 to 11.1 to 33.3 mM. Similar effects of glucose on β-cell proliferation and/or apoptosis have been reported in other studies (25, 42), although on the other hand, Tyrberg et al. (62) did not observe a significant change in β-cell proliferation under hyperglycemic conditions.

It is somewhat difficult to reconcile in vitro observations like those by Maedler et al. (43) with a common-sense understanding of the long-term physiology in normal subjects. Specifically, we find that the quantitative relationship between glycemia and β-cell turnover observed over the course of several days in these studies leads to untenable predictions regarding β-cell population dynamics when extrapolated over longer time intervals. We believe it is very reasonable to assume that, at a young age in a normal subject, the overall rates of apoptosis and replication/proliferation should balance so as to produce neither a progressive increase nor a progressive loss of β-cell mass. Therefore, if we assume, following Maedler et al. (43), that equilibrium occurs at 5.5 mM glycemia (with both apoptosis and replication occurring at 100% of the rate required to maintain this balance), we find that, should the subject's glycemia increase (due to altered lifestyle or to beginning insulin resistance), even to as little as 6 or 6.5 mM (below the threshold for the clinical identification of diabetes mellitus), apoptosis would be predicted to increase and replication would fall. As a result, β-cell mass would be predicted to decrease progressively, insulin production would be affected, and the conditions for a further increase of glycemia would occur. The exponential continuation of this positive feedback loop would determine the eventual development of diabetes mellitus in all subjects who accidentally increase their glycemia for a few days. For the sake of exemplification, we could compute, using the relationships derived from the Maedler et al. (43) data, that a subject with a sustained glycemia of ≥6.5 mM would experience a 10% decrease in β-cell replication and a 30% increase in apoptosis (relative to the previous equilibrium at 5.5 mM). This would result in a net loss of β-cell mass on the order of 40% of 0.43%/day, or ∼0.2%/day or 5%/mo, with reduction of the same β-cell mass to 10% of normal in ∼44 mo. In other words, every subject who is found to maintain a glycemia of ≥6.5 mM should become an insulin-dependent diabetic in <4 yr. Thus, the direct extrapolation of these glucose toxicity in vitro data to the parameters governing the long-term behavior of a model of diabetes progression does not appear practicable.

For the purpose of the present simulations, glucose toxicity (K_{ηG}) has therefore been arbitrarily set as 0.02·mo^{−1}·mM^{−1} glucose concentrations so that (hypothesizing a situation of equilibrium if glycemia remains at 5 mM) at a moderate but permanent hyperglycemia of 10 mM it takes between 2 and 3 yr to essentially abolish the ability of the pancreas to increase net replication rate above zero. Within the same period of time, if a permanent hyperglycemia of 15 mM were maintained, the maximum net replication rate would become significantly negative (of the order of −2%/mo), and continuous β-cell population loss would then result.

#### Insulin sensitivity (K_{xgI}).

This coefficient represents the second-order insulin-dependent glucose uptake rate constant per unit of insulin and therefore reflects insulin sensitivity (K_{xgI}) of peripheral tissues. We have used the sensitivity index (S_{I}) derived from studies using frequently sampled intravenous glucose tolerance tests for this purpose (5, 49, 50). A value for this parameter (1.08 × 10^{−4} min/pM) has also been given for the young (4). This value is higher than the larger values in the Insulin Resistance Atherosclerosis Study dataset (29), but mean age was higher in the latter study (∼57 yr). Minimum and maximum reasonable values were taken from visual inspection of data in Kahn (36), and the ∼20-fold range of S_{I} in this study is consistent with other reports (46). The minimum possible value was taken as zero since a significant number of individuals in the Insulin Resistance Atherosclerosis Study data set did have an S_{I} of 0 (29), although the “zero-S_{I}” problem is a well-recognized methodological limitation of the so-called “minimal model” approach and has been exhaustively discussed (50). The maximum possible value could be a subjective selection of 10 times the typical value. In the simulations shown, however, K_{xgI} is made to decrease, as an external input function into the model, from an initial high value of 1.0 × 10^{−4} min/pM at t_{0} (18 yr of age), with a sigmoidal shape whose steepness and whose time at 50% decrease is arbitrarily made to vary widely depending on the simulated pathological condition (see Fig. 1).

#### Maximal insulin secretion rate (T_{igB}) and maximal insulin secretory capacity (I_{maxB}).

For the purpose of the present model, maximal insulin secretory capacity (I_{maxB}) has been defined as the ratio of maximal insulin secretion per million β-cells per liter of distribution space [maximal insulin secretion rate (T_{igB})] and K_{xi}. In other words, we make a distinction between secretion rate (T_{igB}) and the effect, in terms of concentration, that this maximal secretion rate attains (I_{maxB}). The distinction is relevant when (like for the present simulations) the apparent first-order elimination rate for insulin is made to decline with age, as experimentally observed.

An estimate of maximal insulin secretory rate (ISR) under stimulated conditions was taken from Breda et al. (13). In this work, a maximal mean ISR of 800 pmol/min was observed in healthy, normoglycemic subjects when plasma glucose was artificially raised to 15 mM. This is similar to other reports (26, 51). Assuming 1 × 10^{9} β-cells per pancreas, this corresponds to a typical maximum insulin secretory capacity of 0.8 pmol·10^{6} β-cells^{−1}·min^{−1}. Ferrannini et al. (26) reported maximal ISRs ranging from ∼250 to 1,450 pmol/min in individuals ranging from normal glucose tolerant to T2DM (assuming a BSA of 1.7 m^{2} for all groups). This corresponds to a range of 0.25–1.45 pmol·10^{6} β-cells^{−1}·min^{−1}, assuming equal contribution from 10^{9} β-cells in all groups. Estimates of stimulated insulin secretion rates have also been derived from in vitro studies of human islets. Ling and Pipeleers (41) reported a maximum insulin release rate of 5.6 ± 1.5 pmol·10^{6} β-cells^{−1}·min^{−1} when islets were incubated in 20 mM glucose.

Although insulin secretory response has been shown to decline with age (29), it is not clear whether this represents a decrease in β-cell function (i.e., T_{igB} and I_{maxB}) or a decrease in β-cell mass. In fact, in a study of freshly prepared cadaveric human islets, insulin release in response to high glucose did not differ between younger and older donors, whereas insulin secretion in response to low glucose actually increased with age (34). In a second study, human islets from older donors secreted less insulin than islets from young donors in response to hyperglycemia, yet basal insulin secretion appeared to be preserved (42).

These values could be used to define the minimum and maximum reasonable values for T_{igB} and I_{maxB}. However, given the functional form of the present model, and given the assumption that the system is at equilibrium at t_{0}, the values of T_{igB} and I_{maxB} results are determined once the values of other model parameters (I_{0}, K_{xiStart}, etc.) have been fixed. When setting, as for the present simulations, I_{0} = 50 pM, G_{0} = 5 mM, G_{h} = 9 mM, ν_{h} = 4, V_{i} = 0.25 l/kg body wt, body wt = 70 kg, K_{xiStart} = 0.05 min^{−1}, and B_{0} = 1,000 Mc, we obtain T_{igB} = 0.00287 pmol·min^{−1}·l^{−1}·Mc^{−1}, which translates into a whole body maximal secretion rate of 503 pmol/min. With an initial glycemia of 4.5 mM, we would obtain T_{igB} = 0.00425 pmol·min^{−1}·l^{−1}·Mc^{−1} and a whole body maximal secretion rate of 744 pmol/min so that the range of variation of the insulin secretion rate seems to be perfectly in agreement with the observations cited.

*K _{xi}* is the first-order elimination rate constant for insulin and reflects clearance from the plasma. A typical value and ranges were taken from Walton et al. (63). Minimum and maximum possible values are defined as approximately ±3 SD from the mean in the same study. Insulin clearance rates have been shown to decline with age (35). On the basis of this report, insulin clearance was assumed to follow a linear decline of ∼0.002 min/decade or 2 × 10

^{−4}min/yr in the current simulations (it goes from ∼0.05 at age 18 yr to 0.035 at age 90 yr).

#### Glucose effectiveness (K_{xg}).

This parameter is the first-order insulin-independent glucose uptake rate constant. Insulin-independent first-order glucose uptake into peripheral tissues has been named “glucose effectiveness” in the minimal model analysis (5) of intravenous glucose tolerance test data (where it is generally known as S_{g}). Typical values for K_{xg} may be taken from Clausen et al. (18), and the mean and variance in this study are consistent with another report by Walton et al. (63). Both studies involve relatively large samples (100–200 young, healthy volunteers), and in both studies the estimates of S_{G} or K_{xg} depended on using the minimal model as the mathematical representation of the glucose-insulin system.

Some considerations on this parameter are, however, in order. First of all, other authors (27) have already pointed out that neither K_{xg} nor S_{g} is truly insulin independent since they both reflect basal insulin conditions. Furthermore, common diabetological understanding of glucose elimination from plasma considers two phenomena to be of primary importance: *1*) a second-order insulin-dependent glucose tissue uptake by insulin-sensitive tissues (mainly muscle and adipose tissue), represented in the present work by the coefficient K_{xgI}, and *2*) an essentially zero-order constant (at rest) glucose uptake by the brain and, to a lesser extent, other tissues (in the present work, T_{gl} represents the net difference of zero-order liver production and zero-order brain uptake of glucose). First-order insulin-independent glucose elimination from plasma occurs when glycemia exceeds the renal threshold, in which case glucose is eliminated with urine in a linear, plasma concentration-dependent fashion. When glycemia is below the renal threshold, as occurs in most normal or treated subjects, there is difficulty identifying a physiological process to which a linear term for glucose elimination may correspond. The linear term S_{g} has been included in some modeling representations (5) because in its absence the fit of the model to data was not good. With other models (50) this is not the case, and K_{xg} may well be zero. For these reasons, although provision has been made for possibly including K_{xg} in the current model, its numerical value (and consequently the numerical value of the constant ρ) has been set to zero.

*K _{ag}* is the rate constant of production of Hb A

_{1c}from circulating glucose. Mortenson and Volund (48) estimated the rate constant for the irreversible formation of Hb A

_{1c}in vivo as 7.75 × 10

^{−6}mM/h, approximately equivalent to 0.0265 mo

^{−1}at a glycemia of 5 mM and at a concentration of Hb A

_{1c}of 5% (i.e., at a concentration of Hb A

_{0}of 95%).

*K _{xa}* is the spontaneous elimination rate constant of Hb A

_{1c}. The elimination rate of Hb A

_{1c}is assumed to reflect the elimination of red blood cells. Assuming a normal lifespan of 126 days and zero-order kinetics for red blood cell clearance yields an elimination rate constant of 1/126 day

^{−1}(see 2,107, or ∼0.238 mo

^{−1}). Glycosylation of hemoglobin is assumed to be an irreversible reaction (31).

The two rates above must be reconciled if we believe that at young age and for a considerable number of years the level of Hb A_{1c} remains stable. From *Eq. 8*_{slow} [which reflects the assumption that Hb A_{1c} is continuously produced by glycosylation throughout the life of the erythrocyte (15)], under the assumption of early-life steady state, it is possible to derive the value of K_{ag} once the value of K_{xa} has been specified. Therefore, we have assumed a degradation rate K_{xa} = 0.238 mo^{−1}, computing the corresponding value of K_{ag} from initial equilibrium conditions.

Finally, the following assumptions have been made on the parameter values: at any slow time instant, G, I, and h must be at equilibrium, given the prevailing B and K_{xgI}; this means that, in the numerical integration of the system, a (nonlinear) set of algebraic equations has to be solved at every time point to yield the current G, I, and h values.

t_{0} Was assumed to be age 18 yr.

G_{λ} has been set equal to G_{0}, assuming the subjects to be in a perfect state of health at t_{0}.

From the assumption of steady state of A at t_{0} = age 18 yr,

From the assumption of steady state of B at t_{0}, hence from the assumption that λ(t_{0}) = 0, we have

From the assumption of steady state on η at t_{0}, we have

The half-life of η at 10 mM glucose may be computed as

By definition,

By definition,

From the assumption of fast equilibrium at *t*_{0},

By definition,

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

## RESULTS

### Mathematical Analysis

In the following, a number of statements on the qualitative properties of the model are made, which are either proven in appendix 2 or follow directly from the proofs in appendix 2 using standard computations.

#### Model behavior.

The qualitative behavior of the present model resembles that of the model developed by Topp et al. (61). Assuming positive values of glucose effectiveness at zero insulin [EG0 in Topp et al. (61), corresponding to K_{xg} here], both the βGI system and the present model may present three equilibrium points: a physiological steady state, which is locally asymptotically stable, consisting of the (normal) basal 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}).

The dynamic changes exhibited by the present model are consistent with physiological expectations. In hypoglycemic [G(t) < G_{b}] or moderate hyperglycemic conditions [G_{b} < G(t) < G_{s}], plasma glycemia tends to return to the physiological equilibrium point, G(t) → G_{b}. Imposing a moderate increase in glycemia (as can happen by prolonged inappropriate dietary behavior or by independently worsening insulin sensitivity) results in an increase in β-cell mass and insulinemia; when the saddle point threshold is exceeded [G(t) > G_{s}], glucose toxicity prevails and β-cell mass eventually declines to zero, leading to the pathological steady state with zero insulin production.

From a mathematical point of view, the same pathological steady state exists both in Topp et al. (61) and in the present model; B = 0 is a trivial equilibrium point for the β-cell mass dynamics, to which I = 0 and a positive finite hyperglycemia correspond, provided that the glucose effectiveness at zero insulin (K_{xg}) is strictly positive. The other two equilibrium points (the physiological steady state and the saddle point) are determined by the zeros of a particular function of G. In Topp et al. (61), such a function is a second-order polynomial function with a negative second-order coefficient, a negative value for G = 0 (zero-order coefficient), and a positive first derivative at G = 0 (first-order coefficient). This means that it may have a pair of zeros (corresponding to the physiological steady state and to the saddle point if the larger of the two values is <G_{s}). In the present paper, from the slow dynamics, a nonlinear function χ(G) (see appendix 2) that represents the asymptotic λ at a given level of glycemia is derived. This shares some properties with the second-order function of Topp et al. (61) in that it admits a negative value λ_{min} for G = 0, it increases monotonically for G > 0 until the (unique) maximum is reached, and then it decreases monotonically, approaching asymptotically the same negative value λ_{min} at G = 0. A pair of zeros (corresponding to the physiological steady state and to the saddle point if the greater of the 2 values is <G_{s}) may then be produced. The main mathematical difference between χ(G) and the parabolic function in Topp et al. (61) is that χ(G) is bounded, for all nonnegative G, by prespecified limits on the net rate of β-cell replication, whereas the parabolic function in Topp et al. (61) may reach unrealistically large negative values for high G.

#### Parameter sensitivity analysis.

Concerning the effects of parameter changes in model behavior, there is a substantial similarity between the present and Topp et al. (61) models in regard to the the parameters that have the same meaning. For these parameters, the following considerations hold for both models.

Decreasing the insulin sensitivity does not change the existence (and stability properties) of the three equilibrium points, nor does it change the relative equilibrium glycemias. Conversely, a decrease in insulin sensitivity produces an increase of both equilibrium β-cell mass and equilibrium insulinemia (needed to maintain the same equilibrium glycemia).

Increasing glucose production [R_{0} in Topp et al. (61), T_{gl} here] does not change the existence (and stability properties) of the three equilibrium points. Although an increase in glucose production does not increase the glycemia of the physiological steady state (which is stabilized by larger β-cell mass and corresponding insulinemia), it does increase the pathological (diabetic) glycemia (corresponding to the absence of β-cells and zero insulinemia). A decrease of glucose effectiveness at zero-glycemia K_{xg} produces the same effects on the equilibria as an increase in glucose production T_{gl}.

Increasing the insulin clearance rate does not change the existence (and stability properties) of the three equilibria. No effect is produced in the pathological steady-state levels (obviously, since there is no insulin to be cleared), whereas both physiological and saddle points maintain the same glycemia at an increased β-cell mass and corresponding insulinemia.

For the present model, changes in those parameters describing the dynamics of the pancreatic reserve of β-cell replicating ability do introduce the possibility of a variation in the glycemias of the physiological and saddle equilibria. By decreasing λ_{min} or T_{η} (and hence, the steady-state levels of η from *Eq. 2*_{slow}), G_{b} increases, whereas G_{s} decreases, so that the physiological steady state may be lost and the patient is fatally attracted to the diabetic steady state, whatever the dietary regimen or the insulin sensitivity.

In appendix 2, formal proofs of the main statements used for qualitative analysis are reported.

### Numerical Simulations

Dynamic model behavior concerning the long-term progression of the disease under four different physiological or pathophysiological scenarios is illustrated through simulations portrayed in Figs. 1–13.

Figures 1–7 show the time course of insulin sensitivity, β-cell mass, insulinemia, glycemia, pancreatic replication reserve, β-cell net replication rate, and Hb A_{1c}, respectively. They portray simulations for two individual subjects; the first subject is representative of a “physiological” decrease of insulin sensitivity (solid curves), such as usually occurs with aging (4), whereas the second subject is afflicted by a rather severe early development of insulin resistance and exemplifies one possible evolution toward frank diabetes.

The model predicts that modest and slow decrements in insulin sensitivity are compensated by increases in β-cell mass, insulin production, and fasting serum insulin concentration, with a modest increase in glycemia. These age-dependent trends in fasting insulinemia and glycemia are generally consistent with those reported by Iozzo et al. (35) for a nondiabetic population, although the decrease in fasting insulin reported for the oldest quartile is not predicted in the “normal” individual simulated here. On the other hand, a severe and rapid decrease of insulin sensitivity results in protracted hyperglycemia that accelerates late in life. The model behavior in this individual is consistent with a picture of glucose toxicity, where progressive hyperglycemia exerts a toxic effect on β-cell proliferation that prevents compensation, so that a picture of frank diabetes mellitus eventually develops.

In Fig. 1, the arbitrary time course of insulin sensitivity is portrayed, which serves as the driving function for the model. The time courses of the other state variables in Figs. 2–7 express the predicted reactions to variations in insulin sensitivity. In Fig. 2, the increase, then fall, of β-cell mass in the diseased subject depends on the accumulation of the effects of glucose toxicity on β-cell replication reserve. In Fig. 3, it can be appreciated that, in the diseased subject, insulinemia increases beyond the peak in β-cell mass; this is due to continued and increasing (glycemia-determined) stimulation of β-cell function, which fails only when β-cell mass is severely compromised. In Fig. 4, glycemia is maintained, despite early insulin resistance, at levels <7 mM until age 60 yr, and only then a rapid acceleration of glycemia increase is apparent. This effect does not depend on any external coincident phenomenon but is due to the eventual failure of the adaptation mechanisms that are due to prolonged glucose toxicity. In fact, β-cell replication reserve is seen to decrease as a function of protracted levels of hyperglycemia (>5 mM) in Fig. 5. For illustration purposes, complete absence of therapy is assumed. Net replication rate in the diseased individual responds to hyperglycemia within the limits of the existing pancreatic reserve (Fig. 6); it is thus positive until about age 70 yr, when it falls below zero, as pancreatic reserve is insufficient to even maintain existing β-cell numbers. In the decade between 60 and 70 yr, although positive, λ may be insufficient to produce the increase in β-cell mass that would be necessary to maintain euglycemia in the face of increasing insulin resistance. In evaluating the time course of predicted Hb A_{1c} (Fig. 7) and other state variables, it must be noticed that for illustrative purposes no therapy whatsoever is assumed to have been administered.

In Figs. 8–13, the modifications induced by a progressive decrease of pancreatic replicating ability are depicted. For comparison, in each figure four curves are reported; two (gray) correspond to the cases explored in Figs. 1–7 with normal β-cell replicating ability, indicating with a solid line the case of normal insulin sensitivity and with a dashed line the case of early severe insulin sensitivity. Two more curves (black) correspond to progressively impaired β-cell replicating reserve, again indicating with a solid line the case of normal insulin sensitivity and with a dashed line the case of early severe insulin resistance.

For purposes of illustration, an impairment in β-cell-replicating reserve has been simulated by forcing a decline in the term reflecting β-cell recovery from injury, T_{η}. In Fig. 8, the time course of the T_{η} parameter is shown. For the case of normal replicating ability, this parameter is kept fixed throughout life; otherwise the parameter is made to decrease linearly from a normal initial value to a somewhat decreased value at the end of the period of observation. Solid and dashed lines of the same color cannot be distinguished since they are made to vary in the same fashion. Conversely, Fig. 9 shows the time course of insulin resistance; this is the same as in Fig. 1, and in fact, black and gray lines of the same type (solid or dashed) cannot be distinguished since they are made to vary in the same way.

Figure 10 shows the result of the combination of insulin resistance and/or impaired β-cell replication on the maximal attainable replication rate at any given moment. Clearly, the black lines fall faster since an original defect in replication develops, but it must be realized that insulin resistance determines, at any level of original replicating ability, a further decrement in λ_{max} due to the appearance of progressively more important effects of glucose toxicity.

In Fig. 11, it can be appreciated how the superimposition of impaired β-cell replication onto a situation of near-normal insulin sensitivity leads to a decrease of β-cell mass (solid black vs. solid gray curves) that is slowly progressive. Conversely, the importance of maintaining β-cell replication is manifest when insulin resistance develops concurrently (dashed lines); in this case the coexistence of the two defects impairs the ability to mount an appropriate response to worsening hyperglycemia, with the maintenance of a stress response up to the point where sudden rapid pancreatic failure is apparent. This time course of β-cell mass determines the corresponding insulin levels (Fig. 12), which are essentially normal in the case of isolated β-cell replication impairment but fail to mount and indeed collapse suddenly in the case of coexistence of insulin resistance. The apparent time course of disease, as judged from glycemias, is shown in Fig. 13. It is important here to notice the glycemia scale; for the simulated subject with the particular isolated β-cell replication defect simulated here, glycemia essentially never rises above 7 mM, and the subject could go undiagnosed for the entire period of observation. On the other hand, when replication defects are superimposed on insulin sensitivity, the emergence of a clinical picture of frank diabetes mellitus may be anticipated by a substantial period of time, here by about a decade, and its progression is much accelerated.

## DISCUSSION

The need of a quantitative description of the likely evolution in the compensation state of the glucose-insulin system is acutely felt in diabetological research. Such a description would be very valuable, for instance, in the design of trials testing the efficacy of drugs delaying or preventing the occurrence of T2DM.

Obtaining accurate estimates of the parameters defining the evolution of the system is very difficult because of the very long time scale of the phenomenon, making longitudinal studies impractical, the ethical issues with observing subjects with increasing severity of the disease without intervening, and the practical impossibility of obtaining serial measures of β-cell mass, which is a key element of the compensation state. From this last point of view, imaging and functional tests hold promise and may deliver practical surrogate measures in the near future (54).

The information that we can use today to estimate the β-cell population response to chronic hyperglycemia stems in large part from chronic glucose infusion studies in rodents (30, 38, 39). Although it remains difficult to estimate the relative importance of β-cell apoptosis and replication rate changes in the overall adaptation to changes in (average) glycemia (7), it seems clear that, whereas short-term hyperglycemia induces an increased net replication rate of β-cells, at least in rodents (9), long-term hyperglycemia impairs the same replication rate (66). One possible interpretation of this apparently paradoxical effect, consistent with the mathematical formulation of the present model, could invoke two separate but possibly coexisting actions; on the one hand, hyperglycemia may well stimulate final differentiation of ductal cells into β-cells (6), thus determining a (short-term) stimulation of β-cell population growth, whereas the same hyperglycemia may impair the proliferation of the precursor ductal cells themselves, thus eventually negatively affecting the reserve of differentiable progenitor cells and determining a net loss of β-cells.

Whatever the postulated mechanism, it seems clear that the long-term action of sustained hyperglycemia is to diminish the residual ability of the pancreas to produce new β-cells, termed in the present work “pancreatic reserve.”

Whenever the long-term evolution of glycemia and insulinemia is considered, the exact meaning of blood or serum concentrations becomes unclear. We do not intend to discuss fasting values as such, because the driving stimulus for increased β-cell replication may well consist of the peaks of glycemia attained or of the integral average of the entire glucose time course over a given duration, and represents endogenous and exogenous sources of input, and also because it is clear that insulin effectiveness in bringing about a general reduction in glycemia is related to insulin concentrations attained during the elevations in glycemia induced mainly by meals. It would seem more correct to identify G and I in the model equations as some average quantities throughout the day or week were it not for the fact that the effects of glycemia and insulinemia may depend nonlinearly on the concentration and that, therefore, arithmetic averages may not accurately reflect the effects of the underlying oscillatory phenomenon. It is also clear that whatever experimental verification may be envisaged for this and similar models will in all likelihood make recourse to fasting concentrations sampled very infrequently (say every month) and that continuous measurement of concentrations with the attending numerical integration and averaging is, for the current state of sensor technology, essentially unattainable. Therefore, our point of view is that we are discussing “prevailing” concentrations, and so whenever comparison with data will be possible, these prevailing quantities will have to be approximated in some fashion by the corresponding fasting concentrations or by more sophisticated means, taking into account dietary intake and the like. We, of course, assume that the theoretical time course of prevailing quantities significantly reflects the effective time course of the disease.

The main structural difference between the present model and the theoretical analysis by Topp et al. (61) consists of the explicitly mechanistic modeling of the effect of hyperglycemia on replication on one hand (via *Eq. 6*) and on pancreatic reserve on the other (via the dynamics of η represented in *Eq. 2*). This gives rise to a mathematically consistent model (for whatever values of the state variables) of the effect of sustained hyperglycemia on β-cell net replication and attempts to link the expected replication rate to parameters that are directly physiologically meaningful (like the current state of pancreatic reserve or its spontaneous rate of recuperation). This approach advances the formulation of Topp et al. (61), who had followed a first-line approximation to both replication and apoptosis using segments of parabola indexed by ad hoc numerical parameters (r_{1r}, r_{2r}, r_{1a}, r_{2a}). This approximation has the disadvantages of producing a potentially negative value for either replication or apoptosis (contrary to physiological behavior) and of introducing parameters devoid of immediate physiological meaning.

One of the physiological implications of the present definition of pancreatic reserve (η) is that persisting hyperglycemia could in fact lead to such a decrease in η that, no matter what the prevailing glycemias, the net replication rate would be negative and pancreatic β-cell mass would decrease. This situation corresponds to a terminal state of pancreatic insufficiency and would naturally lead to a clinical picture of frank diabetes, given an adequate time frame. It must be noticed, however, that in this as well as in less dramatic prediabetic conditions, the effect of early therapy would be conducive, through imposed euglycemia, to the restoration of pancreatic reserve and to the delay or avoidance of the final, self-sustaining, positive-feedback deterioration of diabetes. Whether this mathematical interpretation of events actually well represents physiological and pharmacological changes is open to experimental verification, but as a first-line approximation this definition of pancreatic reserve seems to conform well to commonly accepted physiological notions.

A second difference with respect to the model in Topp et al. (61) concerns the spontaneous evolution of the insulin elimination rate K_{xi}. It has in fact been necessary to formally represent the observed decline in K_{xi} over age (35) in order for the model to be able to replicate another observed feature, the natural moderate decrease of insulin secretion in euglycemic subjects (29). It is true that, to obtain the same effect, one could also have considered that, as age advances, body size and in particular lean mass tend to decrease, which in turn would reduce the apparent volume of distribution for insulin and increase insulin concentrations (hence insulin effectiveness) for the same rate of insulin secretion. In order to not overly complicate, however, only the documented relationship of increasing age with decreasing insulin elimination rate was included.

Two additional contributions are proposed in the present work with respect to the published analysis by Topp et al. (61), the formal derivation of the slow model from a more general fast-slow system and the careful assessment of parameter values that would enable the application of the model to human studies.

In appendix 1, the entire derivation of the slow model from a more general fast-slow model, which accounts for both the fast dynamics of glucose and insulin in response to acute perturbations, and the slow evolution of the compensation over years is reported. From this derivation it can be appreciated that many different “fast” models could be in principle reconcilable with the same “slow” model. The only thing that matters, in the passage from different fast systems to the same slow system, is that the different fast processes, whatever they are, quickly reach the same equilibrium state.

This in turn means that the present slow model can be seen as a general description of long-term evolution robust to diverging hypotheses on the actual short-term functioning of the system. Therefore, it is possible, for example, that, in the work by Topp et al. (61), different assumptions on the fast dynamics were actually made (in that work, no explicit mention is made of these), although the end result is similar.

Finally, in the present work a substantial effort has been made to assess, if not statistically estimate, typical parameter values that are coherent with accepted physiology, as shown in the last segment of methods. To this end, the literature has been extensively consulted, and different sources of direct or indirect information on the same model parameter have been critically reconciled. The reasons for the final choice of parameter values have been reported, and although, clearly, different opinions are perfectly acceptable (indeed, ameliorations are to be expected), it is felt that the values obtained here represent, as far as possible, prudent and informed choices. In fact, the end result of such choices is that the model now reproduces naturally such time courses of glycemia and insulinemia over many years as are felt to be compatible with the observed clinical history in representative subjects.

Although we feel that our work represents an improvement in the direction of the approach by Topp et al. (61), we believe that there are substantial differences with respect to the model recently proposed by de Winter and colleagues (19, 20).

The stated modeling goal of both the present approach and that of de Winter and colleagues (19, 20) is to represent the long-term variations in the ability of the individual to adapt to glucose loads maintaining acceptable glycemias. Setting I = INS, G = FPG, and A = Hb A_{1c}, we may rewrite de Winter's equations as where IS stands for insulin sensitivity and BE for β-cell effectiveness (a measure of β-cell insulin secretion at a given glycemic stimulus). As these authors state (20), “ … both IS and BE were modeled as asymptotically decreasing functions of time … ”. In fact, for positive a and b, IS and BE are similar to decreasing exponentials in time. The time courses of both IS and of BE are determined by the chosen values of the respective a's and b's; i.e., they are considered input-forcing functions. Fixing these functions, what is modeled is then the glycemic equilibrium in terms of basal glycemia, basal insulinemia, and Hb A_{1c}. The situation remains qualitatively the same in a later work by the same authors (19), where the parameter IS is renamed S, BE is renamed B, EF_{IS} is changed into 1/E_{S}, and the functional forms of S and B are 1/[1 + exp(a + b × t)], with the same overall behavior as IS and BE.

Although these equations seem to extend previous mass equilibrium representations over a more extended period of time, they should not be taken as representing mass equilibria if I and G represent basal insulinemia and basal glycemia (sampled infrequently, say, once/day) over a period of many years. The automatic transposition of the formalism of mass balance, which is justified in the short term (e.g., 3 h), when transiently elevated glucose concentrations may themselves be responsible for their progressive decline due to insulin-dependent or insulin-independent mass transfer out of the central compartment, is unjustified when modeling the long run of basal glycemia. If we consider the model of de Winter and colleagues (19, 20) as one of the many possible models applicable to describe the short-term glucose regulation, where mass balance considerations do apply, it turns out that it is compatible with the presently introduced slow model. By considering this model at stability by setting dG/dt and dI/dt to zero, we derive that at equilibrium so that we may compute

Similarly,

If we take into account the fact that de Winter and colleagues (19, 20) consider effective glycemia to be its part >3.5, we find that equilibrium insulinemia is given as by (3_{slow}) and as by de Winter and colleagues' model (19, 20). These two expressions correspond term by term when the nonlinearity in the present model's action of glycemia on insulinemia is considered; in fact, the meaning of I_{maxB} and of is the same.

The two models [de Winter and colleagues (19, 20) and the present slow model] are therefore entirely compatible if one considers de Winter and colleagues' (19, 20) model as an instantaneous version of our slow model. In this framework, our model takes into account (and tries to predict) the slow adaptation of β-cell mass as the factor allowing a possible overall compensation. This is another instance of what has been mentioned previously, that our slow model is the slow equilibrium model for the fast models of a wide class. For instance, we might substitute the evolution equation of glucose in de Winter and colleagues' (19, 20) model with the same equation in the minimal model in Bergman et al. (5) and the resulting model would still tend, in the slow dynamics, toward our slow formulation.

We should remark that de Winter and colleagues (19, 20) introduce in a very natural way Hb A_{1c}. Although Hb A_{1c} is not relevant to the development of the overall regulation, it is, on the other hand, a readily observable descriptor of it, and for this reason it has also been taken into account in the present model. However, in the present model, *Eq. 8*_{slow} incorporates a form of the representation of the increase of Hb A_{1c}, where the contribution of the proportion of nonglycated hemoglobin is also taken into account.

A number of issues arose when trying to assess parameter values from information available in the literature. The solutions reached have been detailed in methods, but it is of interest to highlight some important dilemmas.

Regarding the estimates of β-cell turnover, which are key to the assessment of the rapidity with which the β-cell population is able to compensate variations in insulin sensitivity (e.g., produced by dietary or lifestyle variations), the estimated rate of β-cell production via neogenesis is substantially less than what would be predicted from β-cell replication (0.4%/day × 1.14 g of β-cells = 4.6 mg of β-cells/day from replication), and this goes against what several authors maintain, that neogenesis is more important than replication (12, 16, 67). On the other hand, we did not find information about possible attempts by these authors to integrate rate of replication, abundance of each cell type, and the likelihood of differentiating into a β-cell (100% probability for a β-cell but maybe only 2% for a duct cell). If our calculations are right, the estimated relative contribution of β-cell mass and neogenesis from ducts seems to be consistent with reports that ∼15% of β-cell mass is outside of islets and associated with ducts (11). These calculations in fact lead to a dilemma. We assume that rates of β-cell production and β-cell death should be very similar under most conditions (otherwise β-cell mass would typically be increasing or decreasing), but published information seems to suggest that β-cells apoptose at 10 times the rate of renewal under “typical” conditions (see discussion of this in methods). Unfortunately, only Meier et al. (47) provide a readily usable frequency of apoptosis, whereas the large human autopsy series from Butler et al. (16) provides observations as number of apoptotic cells per islet (or as cells·islet^{−1}·%β-cell area^{−1}), and this requires further assumptions to arrive at a frequency of apoptosis. Another potential source of significant error in these calculations could be the assumption that only 2% of replicating duct cells develop into β-cells (16). Further investigation should help decide whether the steady-state assumption is incorrect, the estimates of apoptosis rate are too high, or the estimates of replication rate are too low.

It must be underscored at this point that the model as suggested does not contemplate at all the fast development of (relative) insulin secretion impairment that is seen during acute hyperglycemia in conditions like, e.g., diabetic ketoacidosis. In these conditions, very high levels of glycemia (20 mM and higher), persisting for a few hours or a few days, may determine the inability of β-cells to secrete insulin in appropriate amounts, with the consequent establishment of a fast positive feedback phenomenon, the further rapid increase of glycemia, and the resulting devastating clinical condition. In this situation, however, the administration of insulin together with accompanying support measures like rehydration reduce glycemia to normal levels, with the corresponding quick return to normal of β-cell function, and the resolution of the acute episode. From the point of view of the present model, all of the above events take place over such a short time frame (say, a few days) that the general equilibrium conditions are essentially not affected. It could be an interesting matter for further investigation to consider whether an acute ketoacidotic episode would bring about a decrease in β-cell population, seen as essentially instantaneous in the time frame of the slow model but affecting the subsequent evolution of the compensation.

Similarly, according to a recent report (17), β-cell mass may increase even threefold during pregnancy. The upper bound of λ used in the present work is too small to be able to replicate this extreme rate of increase. Setting a higher upper bound for λ would not be prejudicial to the functioning of the model. In the consideration of the whole lifetime of a subject, a rapid increase in β-cell mass during pregnancy, followed presumably by a similarly rapid decrease of β-cell mass after delivery, would not change the overall time course of the equilibria (hence, its representation would not be necessary) unless pregnancy itself constituted (as may well be) a pathogenic insult or, at any rate, a significant event for the evolution of pancreatic compensation itself. This last phenomenon, however, is not accounted for by the present model and is one of its limitations.

The practical applications of a robust model of diabetes progression invest typically the design of clinical trials and the evaluation of therapeutic regimens. The most fundamental question that may be asked is whether there exists a possibility to intervene early enough in the history of the disease so that the eventual crisis and the development of the frank clinical picture of diabetes mellitus are averted. Indeed, the model we propose does predict aversion of diabetes by early intervention, and hence, judging from the model that we have built, the answer would be yes. If the interaction of glycemia, glucose toxicity, and β-cell replication is qualitatively similar to what has been expressed in the proposed equations, a sufficiently early and effective intervention would spare pancreatic reserve and actually prevent the positive feedback chain of events leading to the manifestation of the disease. To reliably quantify the impact of a given therapy on this course of events, the several parameters on which the model depends must be precisely determined. Furthermore, the interdependence of the parameters in producing clinical features (e.g., the increase of glycemia beyond 7 mM before a given age) must also be assessed, and sensitivity analysis in this direction should be coupled with epidemiologic assessment of the joint distributions of the parameters in the populations of interest. Eventually, simulation of large populations of subjects will be possible, making the present model of physiological response interact with reasonable models of both the development of the primary insult (say, insulin resistance) and of the therapeutic maneuvers commonly employed at the successive stages of the disease.

In conclusion, the present work presents and justifies a mathematical model of long-term diabetes progression. The current effort refines the previous model by Topp et al. (61) and makes the resulting model directly useful for clinical purposes through a careful assessment of the relevant parameters.

The availability of a robust quantitative description of the natural history of the development of type 2 diabetes mellitus should at this point prompt more specific and more easily assessable experimental questions and lead to the optimization of the resources needed to provide the corresponding answers.

## APPENDIX 1: DERIVATION OF THE SLOW MODEL

To clarify why the following equations are written in three different configurations (original, fast, and slow), assume to write the original model for β-cell mass (B), glucose (G), and insulin (I) as follows: (1) (2) (3) (4) (5) (6) where the meaning of the parameters and state variables is explained in Table 1.

Notice that G_{h} and ν_{h} in the function h have been fixed, on the basis of physiological considerations, to 9 mM and 4 [pure number, nondimensional quantity (#)], respectively. This choice produces a sigmoidal function h rising from 0 to 1 such that the central portion of the curve, between 5 (10%) and 10 mM (75%), is essentially linear, such that the portion <3.5 mM is essentially zero (<2%) and such that the curve essentially saturates >25 mM (>98%), with all three of these characteristics having been defined a priori on the basis of the current knowledge of the relevant physiology. Therefore, these two parameters are not free and not to be estimated from observations.

The extra parameter ε may be thought of as being a very small scaling factor. When ε is small, the variations of B and η over time also become very small, and therefore, B may be thought of as being slow in the fast time t, expressed, e.g., in minutes.

We may substitute the fast-scaled time t with a slow-scaled time τ expressed, e.g., in months, by setting τ = εt so that and writing (with a little abuse of notation)

The abuse of notation occurs because we should formally distinguish between the functions, say, B(t) and B*(τ), where B is the β-cell mass function of time in minutes and B* is the β-cell mass function of time in months, with B*(τ) = B(τ/ε) = B(t). However, in the autonomous case (like the present case for B, η, G, and I), this gives rise to differential equations of exactly the same form, and we continue using the same letter for expressions of the value of the same state variable as a function of different time units.

From the above equations, as ε becomes smaller, the left-hand sides of the differential equations defining the fast components tend toward zero, defining average equilibria for the fast variables in the slow time scale.

Now we may simplify our original system, given either by *Eqs. 1*–*4* or by *Eqs. 1*′–*4*′, depending on whether we want to concentrate on the fast or slow dynamics (notice that *Eqs. 5* and *6* always remain the same). As ε goes to zero, the ratio of the two time scales becomes extreme. Under these conditions, from *Eqs. 1*–*4* we obtain the fast subsystem whereas from *Eqs. 1*′–*4*′ we obtain the slow subsystem which, setting and may be rewritten as

It is to be noted that the final form (*Eqs. 1*_{slow}–*4*_{slow}) may be obtained from different formulations of the original problem; in this sense, the slow evolution expressed by *Eqs. 1*_{slow}–*4*_{slow} is compatible with different formalizations of the original and fast subsystems. In particular, delayed glucose and insulin terms in the fast subsystem will be compatible with the same slow system. Translated into physiological terms, this means that different models of the quick adaptation of glucose and insulin concentrations to perturbations (e.g., induced by meals or therapeutic or experimental glucose or insulin injections) are compatible with the same overall model of pancreatic β-cell adaptation over the many years of progression of the diabetic condition of the patient.

To simulate the progressive worsening of insulin resistance, the following very simple disease model has been used (parameters are detailed in Table 1):

Finally, to reproduce declining insulin elimination rate with advancing age, a linear time course of K_{xi} has been assumed:

## APPENDIX 2: PROOFS OF THE ANALYTICAL PROPERTIES OF THE MODEL

It is shown in the following that the model admits a locally stable equilibrium value, and this is independent of the choice of attributing a value of zero or nonzero to the parameter K_{xg} (thereby leaving open the possibility of admitting insulin-independent first-order linear glucose elimination or not).

As a preliminary step, *Eqs. 1*_{slow} and *2*_{slow} will be rewritten by taking explicitly into account the instantaneous relationship existing between B and G.

### Lemma 1

*Eqs. 3*_{slow} and *4*_{slow} provide a bijection between B and G in the ranges (0, +∞) and (0, γ/ρ), respectively.

#### Proof.

Consider the relationship between B and G coming from *Eqs. 3*_{slow} and *4*_{slow}: (A2.1)

For elementary physiological reasons, and as can be easily verified from a mathematical point of view, the time evolution of B is nonnegative so that, according to the constraint (*Eq. A2.1*), a reasonable domain for G is (0, γ/ρ), to which the range (0, +∞) corresponds for B. The bijection is readily verified by explicitly computing the Jacobian dB/dG of *Eq. A2.1*, which is strictly negative for any choice of G in (0, γ/ρ). In case of K_{xg} = 0, the constraint (*Eq. A2.1*) changes as ρ = 0 so that the ranges of B and G both change to (0, +∞). Moreover, the relationship between B and G is still a bijection.

According to Lemma 1, we define Γ: [0, +∞) → (0, γ/ρ], the bijection from B to G such that:

Then, *Eqs. 1*_{slow} and *2*_{slow} may be rewritten as (A2.2)

Let us consider the stability analysis of the equilibrium points of *Eq. A2.2*. If K_{xg} > 0, then B = 0 provides an equilibrium point for η = T_{η} ρ/(γ K_{ηG}). Other equilibrium points have to be found from the solutions of the system: (A2.3) with G in (0, γ/ρ). After standard computation, it becomes readily apparent that and that the Jacobian dχ/dG admits a unique root for G* = G_{λ} ^{3}/x_{0}, corresponding to the unique maximum of χ.

Then, iff (A2.4) there is only one further equilibrium point; iff (A2.5) there are two further equilibrium points; otherwise there are no other equilibrium points.

In the case of K_{xg} = 0, the equilibrium points have to be found among the roots of χ(G) = 0, with G in (0,+ ∞) so that there are two equilibrium points iff (A2.6)

Otherwise, there are no equilibrium points at all.

### Theorem 1

Assume K_{xg} > 0. Then, the equilibrium point [0, T_{η}ρ/(γK_{ηG})] is locally asymptotically stable iff (A2.7)

#### Proof.

Computing the Jacobian matrix, (A2.8) we obtain that

The equilibrium point is locally asymptotically stable iff all eigenvalues of the Jacobian have strictly negative real parts, that is if hypothesis (*Eq. A2.7*) is fulfilled.

### Theorem 2

Assume K_{xg} > 0. If condition (*Eq. A2.4*) is fulfilled, then the equilibrium point with B_{eq} ≠ 0 is locally asymptotically stable. If condition (*Eq. A2.5*) is fulfilled, then the equilibrium point corresponding to the lower value of G_{eq} is locally asymptotically stable; the other equilibrium point is unstable.

#### Proof.

Denote (B_{1}, η_{1}) the equilibrium point. In both cases of conditions (*Eq. A2.4*–*A2.5*), the equilibrium point corresponds to a value of G (namely G_{1}) that is a root of χ(G), with λ(G_{1}, η_{1}) = 0. Then, the Jacobian (*Eq. A2.8*) becomes

from which it follows that the eigenvalues of J(B_{1}, η_{1}) are the roots of the polynomial:

Recall that ∂λ/∂G > 0, ∂λ/∂η > 0, and dΓ/dB < 0 so that the equilibrium point is locally asymptotically stable iff (A2.9)

By exploiting the computation of the Jacobians, condition (*Eq. A2.9*) becomes from which the thesis readily follows.

### Corollary 1

Assume K_{xg} = 0. Then, if condition (*Eq. A2.6*) is fulfilled, the equilibrium point corresponding to the lower value of G_{eq} is locally asymptotically stable; the other equilibrium point is unstable.

#### Proof.

The proof follows the same steps as the proof of *theorem 2*.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2008 by American Physiological Society

## REFERENCES

- 1.↵
- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 17a.↵
- 18.↵
- 18a.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 65.↵
- 66.↵
- 67.↵
- 68.↵