## Abstract

A mathematical model that represents the dynamics of intracellular insulin granules in β-cells is proposed. Granule translocation and exocytosis are controlled by signals assumed to be essentially related to ATP-to-ADP ratio and cytosolic Ca^{2+} concentration. The model provides an interpretation of the roles of the triggering and amplifying pathways of glucose-stimulated insulin secretion. Values of most of the model parameters were inferred from available experimental data. The numerical simulations represent a variety of experimental conditions, such as the stimulation by high K^{+} and by different time courses of extracellular glucose, and the predicted responses agree with published experimental data. Model capacity to represent data measured in a hyperglycemic clamp was also tested. Model parameter changes that may reflect alterations of β-cell function present in type 2 diabetes are investigated, and the action of pharmacological agents that bind to sulfonylurea receptors is simulated.

- β-cell
- insulin secretion
- insulin granule dynamics
- type 2 diabetes

in the past few years, our knowledge of the events involved in the insulin secretion by pancreatic β-cells has advanced considerably. New techniques emerged, such as the use of fluorescent proteins that can be targeted to secretory granules, the real-time imaging of granule trafficking in living cells, and the cell capacitance measurements that allow monitoring of the changes in the cell surface area that result from the fusion of granules with the cell membrane (41).

Since the early observations of the first and second phase of insulin secretion in response to a step in extracellular glucose (11, 13), pathways of stimulus-secretion coupling have been identified. In the triggering pathway, the glucose stimulus causes an increase in the ATP-to-ADP ratio. This leads to closure of ATP-dependent K^{+} (K_{ATP}) channels, cell membrane depolarization, and Ca^{2+} influx through voltage-sensitive Ca^{2+} channels. The consequent increase in cytosolic Ca^{2+} concentration ([Ca^{2+}]_{i}) stimulates exocytosis of insulin granules. Still not completely characterized are other signaling pathways, and several mediators may be involved in the K_{ATP}-independent Ca^{2+}-dependent (amplifying) pathway (43, 21, 31, 30). ATP and Ca^{2+} cover a major role in the machinery that regulates granule translocation from the *trans*-Golgi network to the cell membrane, granule priming, fusion with the cell membrane, and release of contents in the extracellular space (41). Other important aspects of islet behavior are as follows: *1*) the heterogeneity in the responsiveness to glucose of the β-cell population (26) and *2*) the presence of spontaneous sustained oscillations in the concentration of signaling molecules and in the insulin release, revealed both in isolated cells and islets, and in animal models as well as in human patients (19, 40).

Complex mathematical models of glucose-stimulated insulin secretion were proposed previously by Grodsky and coworkers (20, 37) (storage-limited model) and by Cerasi and coworkers (12, 36) (signal-limited model). Subsequently, models of insulin secretion were proposed that provide estimates of β-cell function by using C-peptide data (49, 8, 33) and that are more easily utilizable in clinical practice. Complex models of the Hodgkin-Huxley type that represent the bursting behavior of the β-cell electrical activity and the ATP and Ca^{2+} oscillations were also presented (46, 4, 38).

The mathematical model proposed in the present paper focuses on the dynamics of formation, translocation to cell membrane, and exocytosis of insulin granules in β-cells. We have not attempted to model in detail the biochemical events that lead to changes in ATP concentration ([ATP]) and cytoplasmic Ca^{2+} concentration (or other signaling molecules) in response to stimulation by extracellular glucose or to the experimental manipulations that cause cell membrane depolarization. Instead, these processes are represented as time changes of rate coefficients regarded as “control” signals that regulate granule trafficking. The model provides a preliminary framework in which a number of experimental observations can be accommodated and analyzed in quantitative terms. The capacity of the model to analyze experimental data obtained in a hyperglycemic clamp was ascertained. In addition, to our knowledge, this is the first model that tries to represent the pharmacological action of drugs that stimulate insulin secretion.

## Glossary

### State Variables (Equation Number Where Quantities First Appear)

- I
- Pool of proinsulin aggregates (
*Eq. 1*) - V
- Pool of granule membrane material (
*Eq. 1*) - R
- Reserve pool (
*Eq. 3*) - D
- Pool of docked granules (
*Eq. 4*) - D
_{IR} - Pool of immediately releasable granules (
*Eq. 4*) - F
- Pool of granules fused with plasma membrane (
*Eq. 6*) - γ
- Rate coefficient of granule externalization and priming, related to the ATP-to-ADP ratio (min
^{−1}) (*Eq. 3*) - ρ
- Rate coefficient of granule fusion with cell membrane, related to [Ca
^{2+}]_{i}(min^{−1}) (*Eq. 5*)

### Additional Variables

- C
- Pool of unbound L-type Ca
^{2+}channels - G
- Extracellular glucose concentration (input variable) (mmol/l) (
*Eq. 7*) - ISR
- Insulin secretion rate (pmol/min) (
*Eq. 11*) - IS
- Insulin amount secreted in a time
*T*during the first-phase response (pmol) (*Eq. A4*) - G
_{m} - Measured plasma glucose concentration (mmol/l) (
*Eq. C1*) - I
_{p} - Measured plasma insulin concentration (pmol/l) (
*Eq. C1*) - G
_{a} - Estimated arterial glucose concentration (mmol/l) (
*Eq. C1*)

### Parameters in the Equations of Granule Dynamics

*k*- Rate constant of formation of proinsulin-containing granules in
*trans*-Golgi network (min^{−1}) (*Eq. 1*) - b
_{I} - Biosynthesis rate of proinsulin aggregates at a given glucose concentration (min
^{−1}) (*Eq. 1*) - α
_{I} - Rate constant of degradation of proinsulin aggregates (min
^{−1}) (*Eq. 1*) - b
_{V} - Rate of biosynthesis of granule membrane material (min
^{−1}) (*Eq. 2*) - α
_{V} - Rate constant of degradation of granule membrane material (min
^{−1}) (*Eq. 2*) - τ
_{V} - Time delay related to recycling of granule membrane material (min) (
*Eq. 2*) - C
_{T} - Pool of total Ca
^{2+}channels (*Eq. 4*) *k*_{1}^{+}- Rate constant of association for the binding between granule and Ca
^{2+}channel (min^{−1}) (*Eq. 4*) *k*_{1}^{−}- Rate constant of dissociation for the binding between granule and Ca
^{2+}channel (min^{−1}) (*Eq. 4*) - σ
- Rate constant of insulin release from granules fused with cell membrane (min
^{−1}) (*Eq. 2*)

### Parameters and Functions in the Equations Representing Stimulus-secretion Coupling

- η
- Rate constant in the equation for γ (min
^{−1}) (*Eq. 7*) - γ
_{b} - Basal value of γ (min
^{−1}) (*Eq. 7*) - ψ
- Oscillatory function that represents events inducing oscillations in γ (min
^{−1}) (*Eq. 7*) - h
_{γ} - Function of G representing the activatory action of glucose on γ (min
^{−1}) (*Eq. 7*) - τ
_{G} - Time delay related to time required by glucose metabolism for activation of γ (min) (
*Eq. 7*) - G*
- Glucose concentration threshold for the activation of γ (mmol/l) (
*Eq. 8*) - ĥ
- Maximal value of h
_{γ}(min^{−1}) (*Eq. 8*) - Ĝ
- Glucose concentration over which h
_{γ}remains constant and equal to ĥ (mmol/l) (*Eq. 8*) - ζ
- Rate constant in the equation for ρ (min
^{−1}) (*Eq. 9*) - ρ
_{b} - Basal value of ρ (min
^{−1}) (*Eq. 9*) - h
_{ρ} - Function of γ representing the activatory action of γ on ρ (min
^{−1}) (*Eq. 9*) - k
_{ρ} - Parameter representing the sensitivity of ρ on the activatory action of γ (
*Eq. 10*)

### Additional Parameters and Functions

- I
_{0} - Insulin amount contained in a granule (amol) (
*Eq. 11*) *N*_{c}- Average number of β-cells in an islet (
*Eq. 11*) *N*_{i}- Number of islets in pancreas (
*Eq. 11*) *N*- Total number of β-cells in pancreas (
*Eq. 17*) - f
- Fraction of β-cells responding to glucose stimulus (function of glucose concentration G) (
*Eq. 11*) - f
_{b} - Basal value of the fraction f (
*Eq. A5*) *K*_{f}- Parameter in the function f(G) (mmol/l) (
*Eq. A5*) - V
_{G} - Glucose distribution volume (liters) (
*Eq. C1*) - V
_{C} - Volume of C-peptide accessible compartment (liters) (
*Eq. 17)* - Q
- Cardiac output (l/min) (
*Eq. C1*) - S
_{I} - Insulin sensitivity index [l·min
^{−1}(pmol/l)^{−1}] (*Eq. C1*) - S
_{d} - Dynamic responsivity index (
*Eq. 17*) - S
_{s} - Static responsivity index (min
^{−1}) (*Eq. 18*) *a*and*b*- Parameters in the approximate expression of the first-phase ISR (min
^{−1}) (*Eq. A1*) *x*_{b}- Denotes the basal value of the generic variable
*x* *x̄*- Denotes the steady-state value of the variable
*x*

## DEVELOPMENT OF THE MATHEMATICAL MODEL

### Granule Dynamics in a β-Cell

Let us consider for the moment a single β-cell. A schematic diagram of the formation and trafficking of insulin-containing granules is given in Fig. 1. Essentially, the model describes the kinetics of four intracellular pools of insulin granules: the reserve pool, the pool of docked granules, the pool of immediately releasable granules, and the pool of granules fused with cell membrane. In the following, the model of Fig. 1 will be illustrated and formulated as a mathematical model.

Let us denote by R(*t*) the number of insulin granules in the reserve pool at time *t*. Although it is still unclear how (pro)insulin aggregates are sorted and segregated into granule membranes (35), we assume that the rate of formation of insulin-containing granules in the *trans*-Golgi network may be expressed as *k*I(*t*)V(*t*), where *k* is a segregation rate constant, I(*t*) is the pool of “free” (i.e., not yet segregated into granule membranes) proinsulin aggregates, and V(*t*) is the pool of available free granule membrane material, not yet enclosing proinsulin, at time *t*. A conservation equation for I is written as follows: (1) where the overdot indicates derivative with respect to time, α_{I} is a degradation rate constant, and b_{I} is a production rate that describes proinsulin biosynthesis and aggregation. Changes in proinsulin biosynthesis rate do not appear to be involved in the first and second phase of insulin secretion elicited by a step increase in ambient glucose concentration (42, 22). However, proinsulin synthesis has been reported to increase with glucose concentration (45), so b_{I} is actually a function of the extracellular glucose concentration, G. The governing equation for V is (2) where α_{V} is a degradation rate constant and b_{V} is the rate of production of the granule membranes. The last term in the right-hand side of *Eq. 2* accounts, in simplified way, for the contribution to granule membrane formation resulting from the (complete) recycling of membrane material. This material becomes part of the plasma membrane at exocytosis and is successively removed and returned to the *trans*-Golgi network. We have denoted here by F the number of granules that are fused with cell membrane, by σ the rate constant of the fusion process, and by τ_{V} the constant time interval required for recycling. Incomplete recycling may easily be accounted for by means of a multiplicative coefficient smaller than 1.

To write the governing equation for the reserve pool, R, we assume that granules formed at the *trans*-Golgi network immediately enter the reserve pool and that granules leave this pool to become docked granules with a rate constant γ. If granules are substantially stable, that is, they cannot be disaggregated or degraded, we can write: (3)

The mechanism by which granules are translocated from *trans*-Golgi to plasma membrane in β-cells has not yet been completely elucidated, although it appears to involve unidirectional movement along microtubules by means of motor proteins (52). Experimental data from mouse and rat β-cells show that only 1–2 granules/min per β-cell are released at extracellular glucose concentration around 3 mmol/l, and around 20–30 granules/min per cell are released at the peak of the first phase of insulin secretion after a step increase in glucose (2, 7, 48). Thus granule transport must involve only a small fraction of the reserve pool that has been estimated to contain 9,000–13,000 granules (48). To keep the model simple, we include in the translocation also granule priming, so a unique rate coefficient, γ, will represent the processes that produce docked granules that are competent for release. As discussed in the next section, γ has to be considered a rate coefficient that is actually dependent on the glucose concentration G and represents a main rate-limiting step in the response to glucose.

According to a well-established view (2, 48), within the pool of docked and primed granules here denoted by D, a pool of immediately releasable granules, D_{IR}, must be distinguished. The latter granules are likely to be tightly associated with L-type Ca^{2+} channels (2). Thus we assume (see Fig. 1) that the pool D_{IR} is formed through the binding of docked granules with unbound L-type Ca^{2+} channels, C. For simplicity, we will also assume that the Ca^{2+} channels that become unbound at granule fusion are immediately available to bind new granules. Denoting by C_{T} the constant pool of total Ca^{2+} channels (bound plus unbound), we have C = C_{T} − D_{IR} and we may write: (4) (5) where *k*_{1}^{+} and *k*_{1}^{−} are the rate constants of association and dissociation, respectively, for the binding between docked granule and Ca^{2+} channel, and ρ is a rate coefficient that accounts for the factors that promote the fusion of granules with cell membrane. Actually, granule fusion is strictly related to changes in [Ca^{2+}]_{i}. As seen in the following sections, we consider ρ a rate coefficient that is affected by changes in the ATP-to-ADP ratio or by experimental manipulations, such as the increase in extracellular K^{+} concentration ([K^{+}]), that cause depolarization of the cell membrane and thus an increase in [Ca^{2+}]_{i}.

The rate of granule membrane fusion with cell membrane is given by the last term in the right-hand side of *Eq. 5*, i.e., ρD_{IR}(*t*), so the equation for the pool of granules undergoing fusion with the cell membrane, F(*t*), is given by (6) where σ is the rate constant that regulates the duration of the fusion event. Although insulin cargo release follows the formation of the fusion pore with some delay (41), for simplicity, we lump the two events together by considering a unique rate constant. So σF(*t*) is the number of granules releasing insulin in the unit time at time *t* in a single β-cell. When granules have completed the fusion process and released their content in the extracellular space, the material that forms granule membrane is recycled to the *trans*-Golgi network. This flow of internalized material, delayed by a time τ_{V}, appears thus in *Eq. 2*.

### Stimulation by Changes in Extracellular Glucose

The above model of granule dynamics must be complemented by equations that relate the glucose stimulus to the quantities that govern the machinery of granule trafficking. This has been done on the basis of the scheme in Fig. 2, modified from Fig. 1 in Ref. 21. The ATP-to-ADP ratio, and the other glucose-dependent factors that enhance granule supply from the reserve pool to the docked pool (21, 31, 30), are assumed to be globally represented by the phenomenological parameter γ that appears in *Eqs. 3–4*. [Ca^{2+}]_{i} is represented by the rate constant ρ in *Eqs. 5–6*. The fast oscillations in [Ca^{2+}], related to changes in cell membrane potential, have been neglected in the present study. More comprehensive models of ATP and Ca^{2+} kinetics have been proposed (4, 38, 5).

Experimental evidence suggests that the glucose-induced increase in cytosolic [ATP] stimulates granule movement along microtubules by enhancing kinesin activity (52). Granules are recruited from the reserve pool to the plasma membrane and, in fact, sustained insulin release requires kinesin-dependent granule transport. Granule priming also appears to be accelerated by [ATP] increase (2, 41). Ca^{2+} contributes to granule translocation because kinesin motor activity is increased by the Ca^{2+}-dependent dephosphorylation of kinesin heavy chain (16), and R-type Ca^{2+} channels are implicated in granule mobilization and priming (25). However, for simplicity, we will not include in the model this regulatory role of Ca^{2+} on granule transport. Moreover, it has been found that [ATP] undergoes oscillations, possibly originating in the glycolytic pathway, with a period of 5–10 min (19). By ignoring the actual molecular mechanisms underlying the glucose-induced supply of granules to the docked pool, we assume for γ a first-order kinetics. In writing this equation, we account for: *1*) the existence of a basal value of [ATP], and thus of γ, *2*) the possible presence of spontaneous oscillations, and *3*) the ATP response elicited by the glucose stimulus. Thus we have: (7) where η is a rate constant, γ_{b} is the basal value at low glucose, and ψ is an oscillatory forcing function that represents the events inducing [ATP] oscillations. Glucose activation is modeled by the function h_{γ}, and the time delay τ_{G} is the time required by glucose metabolism (41, 48). Data on isolated β-cells and cell clusters stimulated by glucose (26) show that cell response is maximal and almost constant for G >10–12 mmol/l and much lower below 6 mmol/l. So h_{γ} has been chosen equal to zero for G smaller than a threshold G*, linearly increasing for G in the glucose concentration range (G*, Ĝ), and then constant. We have: (8) where ĥ is the maximal value of h_{γ} and is reached at G = Ĝ.

The increase in the ATP-to-ADP ratio leads to the closure of K_{ATP} channels and increase in [Ca^{2+}]_{i} (see Fig. 2). So, the rate coefficient ρ will depend on γ. In writing the governing equation for ρ, *1*) we account for the existence of a basal value ρ_{b}, and *2*) we assume linearity to model the action of [ATP] on [Ca^{2+}]_{i}. We obtain the following equation: (9) where ζ is a rate constant, and the function h_{ρ}, which embodies the action of ATP on Ca^{2+}, is simply defined as (10) We remark that, for G below the threshold G* and ψ = 0, the rate coefficients γ and ρ stay at their basal values, and the ISR may be thought to be determined by the “constitutive pathway” of insulin secretion. When G exceeds the threshold, h_{γ} is no longer equal to zero, and the ISR becomes determined by the “regulated pathway.”

In summary, as shown in Fig. 2, extracellular glucose modulates the mitochondrial signal γ according to *Eqs. 7–8*, and γ regulates granule translocation and priming. In turn, γ affects ρ, and thus the fusion of D_{IR} granules with cell membrane, according to *Eqs. 9–10*. So, in the present model, it is the parameter γ that has a major role both in the triggering pathway and in the amplifying pathway (43, 21). *Equations 7–10* represent in a very rough way the coupling of glucose metabolism to the signals that activate insulin secretion. In particular, we neglected the dependence of production and degradation rates of ATP on [Ca^{2+}]_{i}, as well as the kinetics of Ca^{2+} stores in the endoplasmic reticulum, which may cause oscillations in the concentrations with a phase shift between γ and ρ (1, 15).

### The β-Cell Population

According to *Eq. 6*, the insulin amount released in the unit time by a single β-cell in the extracellular space is given by I_{0}σF, where I_{0} is the insulin amount contained in a granule, if complete emptying of granules that undergo fusion is assumed. Actually, by targeting fluorescent probes to granule membrane or granule lumen (50), it has been found that membrane fusion is not always associated with insulin release. Thus only a fraction (15–40%) of granule insulin amount, with this fraction being possibly modulated by the stimulus intensity, may be released.

To compute the insulin amount released in the unit time by the whole β-cell population (ISR), we must account for both the size of the population (total number of β-cells) and the fraction of cells that actually release insulin in the given experimental situation. Denoting by *N*_{c} the average number of β-cells per islet and by *N*_{i} the number of islets, the total number *N* of β-cells may be written as *N* = *N*_{c} × *N*_{i}. In the experiments with isolated islets, if a single islet is considered, we have *N*_{i} = 1.

When the intact pancreas in vivo is considered, the number of active β-cells is regulated by several factors. The β-cell number is known to undergo time changes in relation to the balance among cell replication, cell death by apoptosis, and neogenesis by differentiation from ductal epithelium. Regulation of these processes appears to be altered in obese and type 2 diabetic subjects (9, 14). Estimates of the total number *N* of pancreatic β-cells have been reported for mouse (1.06 × 10^{6}) and rat (2.76 × 10^{6}) (6). In studies on human pancreata (14), the number of islet equivalents (IEq, islet with average size of 150 μm diameter) per kilogram of body weight was found to be 6,640 IEq/kg (decreased to 2,641 IEq/kg in type 2 diabetes).

Another important factor is the heterogeneous responsiveness of the cell population to glucose (neglecting other secretagogues). The fraction of responsive islet cells and clusters increases with the extracellular glucose (26), and cell subpopulations with different intrinsic sensitivity to glucose are likely to exist. In addition, because glucose diffuses in the islet volume and is consumed by β-cells, even in isolated islets, the cells will be exposed to different glucose concentrations both in the basal state and during changes in glucose concentration (3). In spite of the rich vascularization of pancreas, extracellular glucose concentration will decrease as the distance from vessels increases, and a wide range of glucose concentrations is likely to be found. Moreover, it is known that insulin itself, acting as a paracrine signal, sensitizes β-cells to glucose stimuli and may thus participate in cell recruitment (47).

To keep the model as simple as possible, we have adopted a rough description of β-cell recruitment by glucose stimulus, by assuming that only a fraction, f, of the total cell population responds to glucose. This fraction will increase as the glucose concentration G increases, with a time delay for recruitment equal to τ_{G}, so we have f(*t*) = f[G(*t* − τ_{G})]. Under the further assumption that β-cells respond to stimulation in a synchronous manner, which may be realistic on a not-too-fast time scale thanks to cell-to-cell electrical and chemical interactions (38), we may finally write for the insulin secretion rate of the whole cell population at time *t* the following expression: (11)

Given the time course of glucose concentration G, *Eqs. 1–11* predict the time evolution of the control rate coefficients γ and ρ, the evolution of the quantities that describe granule dynamics (I, V, R, D, D_{IR}, F), and the evolution of the recruited cell fraction f and hence of the ISR. The following parameters may be considered constant, independent of the degree of stimulation of the β-cell: *k*, α_{I}, α_{V},τ_{V}, C_{T}, *k*_{1}^{+}, *k*_{1}^{−}, σ, τ_{G}, η,γ_{b}, G*, Ĝ, ĥ, ζ, ρ_{b}, k_{ρ}, and I_{0} (parameters of granule dynamics and signaling model); *N*_{c} and *N*_{i} (parameters of the β-cell population). By contrast, other quantities in the model, namely b_{I}, b_{V}, ψ, and f, may be affected by changes in the concentration of glucose (and other secretagogues) or as a consequence of changes in the experimental conditions. An equation for f as a function of G is given in appendix a, and b_{V} will be taken time and glucose independent.

### Steady-state Conditions

To represent the basal condition (glucose concentration G_{b}), or a situation in which glucose concentration is kept constant for a long time at a basal or suprabasal value, we may consider the steady-state solution of the model in which b_{I}, γ, ρ, and f take values corresponding to the assigned value of G, ψ = 0, and the pools I, V, R, D, D_{IR}, and F are constant. These constant values are obtained by setting to zero the time derivatives in model equations. We observe that, in the steady state, we have the following chain of equalities: (12) Taking *Eq. 12* into account, from *Eqs. 1–3* we see that I, V, and R are given in the steady state by the following expressions: (13) Thus we have (14) Using *Eqs. 4–6*, D, D_{IR}, and F can be conveniently expressed in terms of the steady-state value of R: (15) To guarantee the positivity of D in the steady state, the constraint ρC_{T} − γR > 0 has to be satisfied.

For the ISR in the steady state, *Eq. 14* gives showing that, as intuitive, the rate of insulin secretion in the long time is essentially determined by the proinsulin biosynthesis rate, b_{I}. Because b_{I} and f depend on G, the above equation is the “dose-response relationship” predicted by the model when glucose levels (and all other quantities) are kept constant for a long time. As we will see in the following, this relationship may not reflect the ISR data measured from an early part of the second phase of insulin secretion.

### Assessment of Model Parameter Values

For some of the parameters of the model of *Eqs. 1–11*, we may use experimental measurements or estimates, e.g., I_{0} = 1.6 amol in rat granules (41) and *N*_{c} = 1,000–2,000 cells for an islet (2). The granule fusion event plus the insulin cargo release has been reported to require a time of the order of a few seconds (41), so we may take an estimate for σ in the range 6–30 min^{−1}. Moreover, C_{T} = 500 was reported as an upper estimate for the number of Ca^{2+} channels (2). We will take τ_{G} = 1 min (2).

Other indications about the values of model parameters can be obtained in view of the available experimental data concerning the basal state in isolated islets. At the glucose concentration of 3 mmol/l, only 1–2 granules·min^{−1}·cell^{−1} are released in the mouse or rat (48). We will assume that, at the basal state, σF (number of granules released per minute by a β-cell) and the other quantities equal to σF that appear in *Eq. 12* are equal to 1. With this choice, assuming *N*_{c} = 10^{3} active cells/islet and f = 1 (that is, 100% of islet cells release insulin), the ISR of a single islet predicted by *Eq. 11* at the basal state is equal to 1.6 × 10^{−3} pmol/min (or 9.28 pg/min), which is in the range of the values experimentally observed (22, 44, 47).

To determine the baseline values of model parameters, we proceed by taking σF = 1. The pool of immediately releasable granules, D_{IR}, is reported to contain ∼50 granules in the basal state (2, 7, 48). By taking D_{IR} = 50, from *Eq. 12* we estimate that the value of ρ in the basal state, denoted as ρ_{b}, is ρ_{b} = 0.02 min^{−1}. A plausible estimate is that the reserve pool, R, contains 10,000 granules, and that 1,000 granules are docked with the plasma membrane: thus we take D + D_{IR} =1,000. The chosen R value leads to the estimate for γ in the basal state (denoted by γ_{b}) as γ_{b} =10^{−4} min^{−1}. To proceed further, we take as plausible the value of 10 for both I and V at the basal state (it may be seen that this assumption does not affect the general pattern of the first and second phase of insulin secretion). So *k*IV = 1 (from *Eq. 12*) implies *k* = 10^{−2} min^{−1}. Also, from *Eq. 13* we have b_{I}*k*b_{V}/(*k*b_{V} + α_{I}α_{V}) = 1. The quantity *k*b_{V}/(*k*b_{V} + α_{I}α_{V}) ranges from 0 to 1 and increases with b_{V}, so we may look at it as a factor that takes into account the variability in the vesicle membrane production rate related, for instance, to variability in lipid metabolism within the β-cell. Taking for this factor the value 0.25 at the basal state, we deduce for the basal b_{I}, denoted as b_{I,b}, the estimate b_{I,b} = 4 min^{−1} and, using *Eq. 1* written in the basal state, we also obtain α_{I} = 0.3 min^{−1}. Let us now assume that the degradation rate constant for V is larger than the degradation rate constant for I, e.g., α_{V} = 2α_{I}, and use a similar reasoning as above. This leads to the following estimates: α_{V} = 0.6 min^{−1} and b_{V,b} = 6 min^{−1}.

Concerning τ_{V}, *k*_{1}^{+}, *k*_{1}^{−}, and C_{T}, we have tentatively assumed, on the basis of simulation results, the following baseline values: τ_{V} = 5 min, *k*_{1}^{+} = 1.447 × 10^{−5} min^{−1}, *k*_{1}^{−} = 0.10375 min^{−1}, C_{T} = 500. With these parameter values, the (basal) steady-state values of model variables are thus as follows: I = 10, V = 10, R = 10,000, D = 950, D_{IR} = 50, F = 3.33 × 10^{−2}, and ISR = 9.28 pg·min^{−1}·islet^{−1} (f = 1). The values of the other parameters were chosen as discussed in appendix a.

## RESULTS

### Stimulation by High K^{+}

Data of islet response to stimulation by high extracellular [K^{+}] in low glucose and in the presence of diazoxide have been given by various authors (see Ref. 21). Henquin et al. (22) report time courses in mouse islets of [Ca^{2+}]_{i} and of the insulin secretion rate in response to a step increase in extracellular K^{+} and to a sequence of impulsive K^{+} increments with pulses lasting 6 min and interpulse intervals of 6 min. The continuous increase in K^{+} caused a sustained [Ca^{2+}]_{i} rise, whereas the insulin secretion rate exhibited a large initial peak followed by a slow return to the initial value. The impulsive increments in K^{+} caused an intermittent rise of [Ca^{2+}]_{i} and a sequence of decreasing peaks of the insulin secretion rate. The pattern of the response, but not its basal value, was influenced by the glucose level (0 or 3 mmol/l).

This type of stimulation is modeled here by assuming that the time course of ρ reflects the time course of [Ca^{2+}]_{i} reported in Ref. 22, whereas γ retains its basal value γ_{b}. Thus, starting from the basal value ρ_{b}, the time course of ρ in the case of continuous stimulation has been obtained as the response of a linear first-order system with rate constant ζ to a step function (the step increase in K^{+}). A component increasing linearly with slope *s*_{ρ} has been added to mimic more closely the measured [Ca^{2+}]_{i} in the time interval 0–50 min. So, if the step increase in K^{+} occurs at *t* = *t*_{1}, we have ρ(*t*) = ρ_{b} for *t* < *t*_{1} and (16) see the plot in Fig. 3*A*. The values of ρ̂ and *s*_{ρ} were chosen on the basis of simulation results, for a K^{+} increase from 4.8 to 30 mmol/l (22).

The function ρ(*t*) is in this case the input signal for the model of *Eqs. 1–6*. The model response to this stimulation (taking 1,000 β-cells in an islet and f = 1) is shown in Fig. 3, to be compared with data reported in Ref. 22. Figure 3*A* gives the time course of pools D and D_{IR}. The latter pool is seen to be rapidly depleted upon stimulation. The response of I and V (not reported) shows that these pools are affected after the time lag, τ_{V}, required for the recycling from the plasma membrane to the pool V, which is in fact transiently increased. The pool R remains virtually unchanged. The predicted ISR is shown in Fig. 3*B* (solid line). The number of granules released per β-cell, given by the integral of σF over the time, is equal to 67 granules in the first 5 min of the response. According to this model, this type of stimulation does not affect the steady-state values of I, V, R, F, and ISR, which are recovered in the long time: only the steady-state values of D and D_{IR} are decreased (see *Eqs. 13* and *15*). If a Ca^{2+}-mediated enhancement of microtubule-based granule transport (16) were included in the model, high K^{+} stimulation would produce a component of granule translocation to plasma membrane with a decrease in R.

For the intermittent stimulation, Fig. 3*C* shows the time course of ρ (a sequence of 4 pulses lasting 6 min with values of ρ̂, ζ, and *s*_{ρ} being the same as for the continuous stimulation). The predicted response is illustrated in Fig. 3, *C* and *D*. Note in Fig. 3*C* the refilling of pool D_{IR} during the interpulse interval. When C_{T}, *k*, and *k* were set to alternative values that produce the same amount of granules in the various pools at the steady state, the pattern of the response was qualitatively similar. The ISR is reported in Fig. 3, *B* and *D*, by the dotted line (83 granules released in the first 5 min). The simulated ISR of the islet in Fig. 3*D* reproduces the declining amplitude of the peaks experimentally observed at zero glucose (22). The data show a slight broadening of the peaks, possibly related to the incomplete synchrony of β-cells in an islet.

### Stimulation by Extracellular Glucose

As seen in the previous section, insulin secretion induced by K^{+} rise can be simulated in the model by an increase in the rate coefficient ρ and shows a first-phase response only. In that case, the ISR pattern is monophasic, and, after the peak, the secretion rate declines toward the basal value. The secretory response to a step increase in extracellular glucose is, instead, biphasic. The second phase of glucose-stimulated insulin secretion had been initially considered to be sustained by proinsulin biosynthesis. Subsequently, this interpretation was refuted (18, 42). Confirmed by model simulations (data not shown), an increase in the rate of proinsulin biosynthesis b_{I}, concomitant with the increase in ρ, does not produce a secondary response if γ is kept at γ_{b}. As previously discussed, a critical role is played by the increase in the rate coefficient γ modulated by extracellular glucose, because γ intervenes in both the triggering and the amplifying pathway (21, 43).

To compute the model response to glucose stimuli, the functions h_{γ}(G) and f(G) have to be specified. To this aim, we rely on data of insulin secretion by the perfused rat pancreas, reported by Grodsky (20). In particular, we used the data of insulin amount released during the first phase of the response to different values of a step change in glucose concentration (IS, see Fig. 4). An approximate expression of the first-phase ISR and of the insulin amount secreted, as predicted by the present model, is derived in appendix a. appendix a also gives the form chosen for the function f (*Eq. A5*) and describes the procedure followed to estimate the parameters of h_{γ}(G) and f(G) from Grodsky's data. Figure 4 shows the IS data (20) with the fitting curve of *Eq. A4* (parameter values reported in appendix a). The number of β-cells used in the simulations in this and the following section is equal to the estimate reported for the rat.

Figure 5 shows examples of the model response to changes in glucose concentration. Figure 5*A* shows the simulated step change of G from 1 to 16.7 mmol/l (300 mg/dl in Grodsky's data) with the corresponding time course of D and D_{IR}. Figure 5*B* reports the corresponding ISR (solid line), together with the ISR predicted with step in G equal to 150 and 500 mg/dl. In the 16.7 mmol/l glucose step, 570 granules per β-cell are released in 50 min, to be compared with the value of 680 granules reported in Ref. 48 for rat pancreatic islets that exhibit a large second phase. The ISR predicted by the model exhibits both the first and the second phase of insulin secretion. However, the amplifying effect can be “deleted” in the model by letting γ and ρ evolve according to *Eqs. 7–10* but imposing that γ is kept to its basal value γ_{b} in the *Eqs. 3–4* that govern granule translocation. In this case, the response only shows a first phase, followed by a decline in the ISR (Fig. 5*B*, dashed-dotted line). The difference between the solid and dashed-dotted line in Fig. 5*B* may thus be related to the amplifying action of glucose on insulin secretion. Figure 5, *C* and *D*, shows the response to two other patterns of glucose stimulation. The resensitization effect resulting from a previous prolonged stimulation is shown in Fig. 5*C*, and the model prediction agrees fairly well with experimental data (20). The ISR patterns obtained when G has a fast or a slow ramp rise are shown in Fig. 5*D*, where the delay in the response is because of the time required to reach G*. As also found experimentally (20, 37), the ramp stimulus greatly reduces or virtually abolishes the first phase of insulin secretion.

A point to stress is that neither the increase in the proinsulin biosynthesis rate b_{I}, which has been observed to follow glucose increase (45), nor the recycling of granule membrane material appreciably affect the second phase of the response in the short time. A large increase in b_{I}, concomitant with the increase in G, does not modify appreciably the ISR response in a time interval of 1–2 h after the glucose step. By contrast, as expected from *Eq. 14*, a rise in b_{I} modifies the ISR in the long term. Figure 6*A* reports the response to a step increase in G computed up to 12 h. After a phase of ∼3 h in which the secondary response increases, the ISR tends slowly to the initial value with a marked depletion of the pool R. If, upon glucose stimulation, the biosynthesis rate is increased, the ISR tends to a larger value with a diminished depletion of the pool R, which is now refilled at a greater rate by proinsulin biosynthesis (see *Eq. 13*). The b_{I} increase with G will be maintained in the simulations that follow. A large value of τ_{V}, which slows down granule membrane recycling, tends to delay this effect.

The simulations of Fig. 6*A* show that the ISR measured at 1–3 h after the glucose step, which may be assumed to represent the stationary response and reported in a dose-response curve ISR vs. G, might be different from the actual stationary response of the ISR. In fact, because the long response time of the second phase related to the small value of γ (see appendix b), the stationary state may be attained only after a very long time.

Figure 6*B* shows the effect on ρ and ISR of a sinusoidal oscillation possibly arising in the glycolytic pathway, represented by the function ψ in *Eq. 7*. Because the oscillations of β-cell signaling molecules appear to be enhanced by extracellular glucose (19), we tentatively assumed the amplitude of ψ equal to (1/2)h_{γ}(G) (G changed stepwise from 1 to 16.7 mmol/l). The oscillation period was set to 10 min. When G exceeds the threshold G*, the oscillation becomes clearly visible in both ρ and ISR. Depending on the amplitude, frequency, and site (oscillations may also be present in the fraction f of active β-cells), the oscillatory component may represent a substantial or a major part of the ISR (40). In Fig. 6*B*, we also simulate the effect of the experimental manipulation consisting in the increase of extracellular K^{+}, in the presence of diazoxide, when glucose is at a stimulatory concentration (19). As in the preceding section, K^{+} rise is modeled by an increase in ρ (we set ρ = 2 at *t* = 4 h). When ρ is kept at a high value, the second-phase response is preserved since γ remains at a high level and oscillates, but the oscillations in the ISR disappear. This behavior is consistent with the observations (19).

### Simulation of Derangements of Physiological Mechanisms and Effect of Drugs

As seen in appendix a, the product k_{ρ}ĥ, but not the values of the two factors, was estimated from Grodsky's data. Figure 7*A* shows that, when ĥ is increased and k_{ρ} decreased with k_{ρ}ĥ kept constant, the first-phase response to a step increase in G is unchanged, whereas the second phase has a larger amplitude and a more prompt decay. Pool R is more completely depleted as ĥ increases. This behavior agrees with the analysis in appendix b and may be seen as the result of an increase in the sensitivity of second-phase ISR with respect to glucose. The pattern of the response at large ĥ, where the second phase may largely exceed the peak of the first phase, is more reminiscent of the ISR observed in the rat than of that observed in the mouse (22). For a given G step, a larger ĥ causes a larger rise in γ, and this agrees with the idea that the amplifying action may be stronger in the rat than in the mouse (22). Figure 7*A* also shows the ISR obtained with the alternative values of C_{T}, *k*_{1}^{+}, and *k*_{1}^{−}, given in the legend of Fig. 3. Because the product *k*_{1}^{+}C_{T} is larger with the alternative parameters than with baseline parameters, the second-phase ISR has a faster increase, as predicted by the initial value of the derivative (see appendix b).

In the simulations of Fig. 7*A*, the parameter changes modified only the second phase of the response. Figure 7*B* shows the effect of changes that affect both the first and the second phase. The dotted-line plot reports the ISR response predicted with b_{V} = 1.5 b_{V,b}, which may represent an increased synthesis rate of granule membranes. The basal values of granule pools, as well as the basal ISR and the whole ISR response, are increased. In contrast, the solid-line plot was obtained by assuming a decreased effectiveness of stimulus-secretion coupling. A decreased ĥ may simulate an impaired effectiveness of processes that regulate ATP metabolism, such as a loss of function of glucokinase, and decrease in k_{ρ} simulates impairment of Ca^{2+} transport in the cell related, for instance, to a defective closure of K_{ATP} channels. Moreover, a smaller σ may represent impairment of processes involved in the exocytotic machinery (41). With ĥ, k_{ρ}, and σ decreased with respect to the baseline values (Fig. 7*B*, solid-line plot up to *t* = 4 h), the glucose-stimulated ISR almost completely lacks the first phase, and both the second phase and the oscillation amplitude are markedly reduced. These parameter changes may reflect alterations in β-cell function present in type 2 diabetes (41).

In the solid-line plot of Fig. 7*B*, we also simulate, within the extremely simplified representation of the stimulus-secretion coupling adopted in the present model, the action of pharmacological agents that bind to sulfonylurea receptors (27, 24) and cause closure of K_{ATP} channels, leading to increased conductance of voltage-sensitive Ca^{2+} channels and increased Ca^{2+} influx. An increase in the model parameter k_{ρ} translates into an increase in ρ, that is, in the [Ca^{2+}]_{i}. Thus, the drug action may be simulated by an increase in k_{ρ}. The solid-line plot of Fig. 7*B* shows, up to *t* = 4 h, the response to a glucose step of a β-cell with impaired stimulus-secretion coupling. At *t* = 4 h, to simulate the drug effect, k_{ρ} is exponentially increased to the baseline value. This change causes a rise in ρ and thus a transient increase in the ISR, as it would occur with a rise in extracellular K^{+}. The amplitude of the oscillations is restored.

In general, the effect of antidiabetic drugs should guarantee that the signal ρ is large enough to allow exocytosis of the docked granules that are fed by the reserve pool at rate γR, as remarked in Ref. 21. If the quantity ρC_{T} − γR in *Eq. 15* for D becomes small (and positive), the pool D may in fact become exceedingly large.

### C-peptide Secretion During Hyperglycemic Clamp in Humans

As a preliminary validation of the present model, we simulated the C-peptide response to a hyperglycemic clamp in healthy humans. Previously obtained data of glucose, insulin, and C-peptide were used (32). The C-peptide secretion rate was computed according to the model of *Eqs. 1–11*, and the whole body C-peptide kinetics were represented by the usual two-compartment model (17, 39) with the standard kinetic parameters computed for each subject (51).

To determine the time course of glucose concentration that actually stimulates the in vivo pancreatic response, we estimated the time course of arterial glucose concentration, G_{a}, as shown in appendix c. In the insulin secretion model, the number *N* of β-cells in human pancreas was set according to the estimate in Ref. 14. A parameter identification was not attempted; instead, a few parameters were modified with respect to the baseline values previously reported. In particular, we used the alternative values for *k*_{1}^{+}, *k*_{1}^{−}, and C_{T}, and we set ĥ = 1.38 × 10^{−3} min^{−1}, a value smaller than that used in the simulations of rat pancreatic response. Other parameters were chosen in the following ranges: ρ_{b} = 0.01–0.02 min^{−1}, k_{ρ} was two to four times the value in the rat, and b_{V} = 6–8 min^{−1}. We found that C-peptide data were adequately approximated by assuming f constant at a low value (f = 0.06–0.07). This may be related to a very slow recruitment of the β-cells during the course of the experiment (fasting subjects) or to inhibitory signals acting on the islets in vivo. To obtain the observed value of C-peptide preceding the clamp, we added to the ISR a “constitutive” insulin secretion rate given by I_{0}σFf*N* with σF = 1 and f = 0.05.

Figure 8 shows, for one of the subjects, the experimental data of glucose and C-peptide together with the time course of G_{a}, reconstructed according to *Eq. C1*. Figure 8*B*. reports the time courses of the ISR and of the C-peptide concentration as predicted by the model when the glucose stimulus in *Eq. 7* is G_{a} (solid lines) or G_{m} (dotted lines). The initial more rapid increase of G_{a} with respect to G_{m} may be a factor that accounts for the occurrence of an appreciable first-phase ISR and the observed rapid increase of the C-peptide concentration. No significant differences in the C-peptide concentration calculated by the model with G_{a} or G_{m} were observed, likely because of the long time constant of plasma C-peptide kinetics (51).

### Interpretation of β-Cell Responsivity Indexes

Using the approximate expressions of the first- and second-phase ISR given in appendixes a and b, we may derive the expressions for the dynamic and the static responsivity index of the β-cell, denoted as Φ_{d} and Φ_{s} in Refs 8, 49, and 10, that are predicted by the present model for a step increase in glucose concentration.

According to Ref. 10, we define a “dynamic responsivity index,” S_{d}, as the ratio where IS is the insulin amount secreted during the first-phase response, given by *Eq. A4*, and V_{C} is the volume of the C-peptide-accessible compartment. Note that G* has a meaning similar to the threshold parameter h in Ref. 10. The index S_{d} is dimensionless, like the index Φ_{d}. As seen from *Eq. A4*, the index S_{d} depends on G, so we consider the case of a small G step (a different expression may be derived for large G). By approximating exp(−*aT*) in *Eq. A3* as 1 − *aT*, we get *Z*(G) = ρ̄D_{IR,b}*T*. Furthermore, neglecting ρ_{b} in ρ̄, we obtain from *Eq. A4* the equation for S_{d}: (17) The above expression shows that S_{d} essentially depends on the size of the pool of immediately releasable granules and on the product of the [ATP] sensitivity to glucose, ĥ/(Ĝ − G*), times the [Ca^{2+}]_{i} sensitivity to [ATP], k_{ρ}. These parameters are defined in *Eqs. 8* and *10*.

The “static responsivity index,” in the spirit of the index Φ_{s} in Ref. 10, takes into account the dynamics of the second-phase ISR. Thus we have obtained a dynamic index by computing the derivative at *t* = 0 of the second-phase ISR and dividing by the rate constant, *k* C_{T}, of ISR increase. From *Eqs. B3* and *B4*, we get for this ratio the quantity I_{0}γ̄R_{b}f*N*. Proceeding as for the dynamic index, the static responsivity index S_{s} (with dimension of time^{−1}) has the expression (18) that depends on the [ATP] sensitivity to glucose and the size of the reserve pool. As noted in appendix b, *Eq. 18* is obtained under the assumption of large G (free [Ca^{2+}]_{i} at its maximum level).

With the values of model parameters used in the present study to simulate the response to the hyperglycemic clamp in a normal human subject, we obtain S_{d} ≈ 670 × 10^{−9} and S_{s} ≈ 26 × 10^{−9} min^{−1}. These values are in the range of the Φ_{d} and Φ_{s} values estimated in Refs. 8, 49, and 10 under different experimental conditions and by means of a different model.

## DISCUSSION

In this study, a model of insulin granule trafficking in pancreatic β-cells is proposed. Essentially, the model describes the formation of insulin granules and the dynamics of four intracellular granule pools: the reserve pool, the pool of docked granules, the pool of immediately releasable granules, and the pool of granules fused with cell membrane. Granule translocation from the *trans*-Golgi network to the cell membrane and exocytosis of granules fused with the cell membrane are regulated by two rate constants, γ and ρ, assumed to be related to the ATP-to-ADP ratio and other mitochondrial signals, and to the free [Ca^{2+}]_{i}, respectively. As anticipated by Grodsky (20), the “labile” insulin may be grossly identified with the insulin contained in the D_{IR} pool. The “stable” insulin is the R pool, with the D pool, that feeds the D_{IR} pool as insulin is secreted, being in an intermediate position. The model links in a quantitative, although simplistic, way the dynamics of granule release by the single β-cell, with the recruitment of a subpopulation of glucose-responsive cells within the whole β-cell population.

The different roles of the triggering pathway and of the amplifying pathway in the glucose-stimulated insulin secretion are modeled according to Ref. 21. The rate coefficient γ, which represents events leading from glucose influx into the cell to ATP production, has an essential role in promoting the granule externalization and priming that supports the second-phase ISR. This rate coefficient, in particular the parameters η and ĥ in its governing equation, is a rate-limiting step in the glucose-stimulated insulin secretion. The representative parameters of the part of the triggering pathway from the ATP-to-ADP ratio to [Ca^{2+}]_{i} are ζ and k_{ρ}, which thus represent another rate-limiting step in the insulin secretion (48). We stress that the model does not include the regulatory role of Ca^{2+} in granule translocation and thus in the second-phase insulin release (16, 25), nor Ca^{2+} implication in ATP metabolism (1, 15). Moreover, the events on the fast time scale of the rapid fluctuations in β-cell membrane potential (38, 46) are not accounted for.

The pattern of the pancreatic response to glucose is shown in the present study by the expression of the ISR in the first and the second phase (see appendixes a and b), by the simulation tests, and by the responsivity indexes of *Eqs. 17–18*. The first phase mainly depends on the insulin amount in the D_{IR} pool, and on the product of [ATP] sensitivity to glucose, ĥ/(Ĝ − G*), times the [Ca^{2+}]_{i} sensitivity to [ATP], k_{ρ}. Accordingly, the dynamic responsivity index S_{d} accounts for the fusion of D_{IR} granules (release of labile insulin), essentially regulated by the increase in [Ca^{2+}]_{i}. Notably, agents that produce membrane depolarization, e.g., through closure of K_{ATP} channels, or potentiate the triggering pathway at some point, such as gastric inhibitory peptide (GIP) and glucagon-like peptide (GLP)-1, may result in an increased S_{d}. This action is, however, transient if it is not accompanied by adequate increase in the proinsulin synthesis rate and provision of granules competent for release. In contrast, the static responsivity index S_{s} mainly accounts for granule provision, which is regulated by the size of the reserve pool and the increase in the ATP-to-ADP ratio. The rate constant γ is a synthetic descriptor of the signals that regulate these events, and its increase, together with the increase in b_{I}, is able to produce a durable release of insulin without degranulation of the β-cell.

The secretory response is also determined by the number, f*N*, of responsive β-cells. Given the total number *N* of β-cells in pancreas, possibly subject to slow but considerable changes in pathological conditions (14), the fraction f of recruited cells will change with the time in response to a variety of stimuli (glucose, insulin, GIP, GLP-1, etc.). The complex phenomena and interactions among the cells distributed in pancreas volume [for instance, the diffusion and consumption of glucose in the islets with consequent heterogeneity of glucose concentration (3) and electrical and chemical interactions (38)] are not truly represented in the present model. These phenomena are described by a dependence of the fraction f of active cells on G, at most with a time delay. It may be speculated that cell recruitment has a slower dynamics, also related to the concentration of insulin and other substances in the interstitial fluid that surrounds the β-cells. In vivo, because of inhibitory signals, the maximal value attained by f is likely to be substantially smaller than the unity, and was estimated to be of the order of 6–7%.

Although the present model does not allow to discriminate among the various pathways involved in the pathophysiology of type 2 diabetes, it might provide some insight into the coupling between energy changes and insulin secretion. Therefore, should the mitochondrial derangement described in type 2 diabetes for skeletal muscle (29) occur also in β-cell, then the parameters regulating the kinetics of the signals γ and ρ would be altered with a subsequent reduction in insulin secretion. The problem also arises of understanding how the various events involved in the trafficking of intracellular granules are coordinated and how the physiological amount of granules in the various pools is maintained over time. There may be control mechanisms devoted to “sensing” the size of granule pools to avoid granule depletion or excessive crowding and devoted to check the status of the network of microtubules that vehicle granule translocation.

The model proposed allows a good fitting of experimental data obtained during a hyperglycemic clamp, suggesting that it may be applied to in vivo experimental conditions where glucose is stably increased over time. However, the experimental pattern of the ISR and the value of the responsivity indexes, determined for instance in an oral glucose tolerance test (8, 10), will also be determined by factors other than glucose, such as the signals arising in the enteroinsular axis. The present model can represent the effect of these signals through changes in the equations for the control variables γ and ρ. For instance, the incretin effect may be translated into an increase in k_{ρ}. The k_{ρ} increase may also be related to an increase in the Ca^{2+} efficacy on exocytosis.

A similar consideration is also true for the effect of pharmacological agents on insulin secretion. The mechanism of action of different drugs that stimulate insulin secretion can be modeled, at least partially, by taking into account the drug effect at different levels, for instance on the K_{ATP} channels or on the metabolic energy level as it happens for acetylcarnitine. According to the observations in Ref. 21, it is important that drugs stimulating insulin secretion act not only at the level of the amplifying pathway but also of the triggering pathway; otherwise, the final effect would be that of an imbalance in the equilibrium between these pathways. Therefore, tailor-made drugs in diabetes should be designed to obtain the maximum effect.

A full description of these events cannot be given by the present model because the representation of the signaling pathways involved in the insulin secretion is still oversimplified. For instance, the role of long-chain acyl-CoAs and protein kinases (43, 53) is not taken into account. A more detailed description of the kinetics of signaling molecules should improve the capability of the model to investigate the events that lead to β-cell dysfunction and predict the effects of agents stimulating insulin secretion.

## APPENDIX A: ANALYSIS OF FIRST-PHASE RESPONSE

We derive here an approximate expression for the insulin amount released in the first-phase response following a step in glucose concentration. The parameters appearing in this expression will be estimated from data reported by Grodsky (20) (see Fig. 4). We assume, for the sake of simplicity, that the dynamics of γ and ρ are much faster than granule dynamics (large η and ζ values). Although the measured [ATP] and free [Ca^{2+}] in β-cell cytosol after a glucose step show rise times of a few minutes (1), glucose appears to provoke the formation of subplasma membrane ATP microdomains, where ATP rapidly responds to a glucose increase with a half-time = 20 s (28). The state of nearby K_{ATP} and L-type Ca^{2+} channels is likely to rapidly change, leading to fast Ca^{2+} influx at sites close to the granules of the D_{IR} pool. The transient increase in [Ca^{2+}]_{i} at these sites, and thus in ρ, may be completed in 1 min. We will take as baseline values η = ζ = 4 min^{−1}, so the transient of γ and ρ will require ∼1 min. Simulations of a recently proposed model for mitochondrial ATP production (5) show that mitochondrial ATP and [Ca^{2+}] rapidly respond to a step increase in fructose bisphosphate, with a transient that is completed in <1 min.

Following the step change in glucose concentration to a suprabasal value G, and after the time delay τ_{G}, γ and ρ will rapidly attain their steady-state values, denoted as γ̄ and ρ̄, respectively. From *Eqs. 7* [with ψ(*t*) = 0] and *9*, we have γ̄ = γ_{b} + h_{γ}(G) and ρ̄ = ρ_{b} + k_{ρ}h_{γ}(G). Let us take the time origin so that γ(*t*) = γ̄ and ρ(*t*) = ρ̄ for *t* > 0. Assuming that pool D does not change appreciably during the first-phase response, we solve *Eq. 5* with D(*t*) held constant at its initial (basal) value D_{b}, ρ equal to ρ̄, and initial condition D_{IR}(0) = D_{IR,b}. The pool D_{IR} will be depleted according to (A1) where *a* = *k*_{1}^{+} D_{b} +*k*_{1}^{−} + ρ̄, with ρ̄ being a function of G, and *b* = *k*_{1}^{+}D_{b}C_{T}. Then we solve *Eq. 6* for F, with initial condition at the basal value F_{b}, ρ = ρ̄, and D_{IR}(*t*) given by *Eq. A1*, obtaining (A2)

To compare the prediction of the present model with the data of insulin amount released by the perfused rat pancreas during the first-phase response, we must compute the integral of the ISR over a time interval, say of length *T*, that just contains the first ISR peak. Recalling *Eq. 11*, we compute the total number of insulin granules released by a single β-cell in the interval *T*, that is, the integral of σF(*t*) from *t* = 0 to *t* = *T*. Assuming *a* ≪ σ and F_{b} small, this integral, denoted as *Z*(G), is given by (A3) Thus the insulin amount, IS(G), released in the first-phase response to a glucose step G can be written as (A4) where *N* here is the number of β-cells in rat pancreas. The function f was chosen of the form (A5) where f_{b} is the fraction of responding cells at glucose concentration smaller than the threshold G* (for simplicity assumed equal to that in *Eq. 8*), and *K*_{f} is a constant related to the effectiveness of recruitment. According to *Eq. A5*, ≈100% of β-cells may be recruited for a large enough G.

*Eq. A4* was used to fit Grodsky's data. Taking the baseline values of model parameters, we have still to estimate the quantities G*, Ĝ, the product k_{ρ}ĥ, f_{b}, and *K*_{f}. To simplify the estimation problem, we set Ĝ = 10 mmol/l according to data in Ref. 26 and f_{b} = 0.05. Least-squares fitting, used to estimate the remaining parameters, provided the following values: G* = 4.58 mmol/l, k_{ρ}ĥ = 1.38 min^{−1}, and *K*_{f} = 3.43 mmol/l (see Fig. 4). By model simulation we chose, as baseline parameters, ĥ = 3.93 × 10^{−3} min^{−1} and k_{ρ} = 350. Finally, we note that a simple expression for the first-phase ISR can no longer be obtained if the dynamics of γ and ρ are not much faster than granule dynamics. Model simulations show that the peak ISR becomes smaller and is delayed in time when η and ζ decrease, but the peak values in Fig. 5*B* can be recovered by increasing the product k_{ρ}ĥ. If, for instance, η and ζ are decreased from the baseline value of 4 min^{−1} to η = ζ = 1 min^{−1}, the baseline value of k_{ρ}ĥ must be multiplied by 2.5.

## APPENDIX B: ANALYSIS OF SECOND-PHASE RESPONSE

By a similar analysis as in appendix a, an approximate expression for the second-phase ISR can be obtained as follows. Let us assume that, during this phase, the pool R decays exponentially from the basal value R_{b} to a final value R̄, where R̄ is expressed as in *Eq. 13* with the final values of γ and b_{I},γ̄ and b̄_{1}, respectively, and the other parameters unchanged. The rate constant of the exponential is γ̄, so we have (initial time *t* = 0) (B1) Using the quasistationary approximation Ḋ_{IR} = 0 and summing the equations for D and D_{IR}, we may write (B2) For ρ̄ >> *k* D + *k*(valid when G is large and the concentration of free Ca^{2+} is at its maximum level), the differential equation in *B2* may be written as a linear differential equation for D that can be solved using *Eq. B1* and with D(0) = D_{b}, obtaining (B3) By the further quasistationary approximation Ḟ = 0, we finally get (B4) that may be used with *Eq. B3* to compute the approximate expression of the second-phase ISR. The derivative of σF(*t*) at *t*= 0 is given by *k* C_{T} (γ̄R_{b} − *k* C_{T}D_{b}).

If *k* C_{T} = 7.24 × 10^{−3} min^{−1}, and γ̄ ≤ γ_{b} + ĥ = 3.93 × 10^{−3} min^{−1} (baseline parameters), the descending front of second-phase ISR has a large time constant (1/γ̄ > 4 h). The time constant becomes smaller if ĥ is larger, as shown in Fig. 7*A*. With the baseline parameters, the negative term in the above derivative may be neglected for G >10 mmol/l.

## APPENDIX C: ESTIMATION OF ARTERIAL GLUCOSE CONCENTRATION

An estimate of arterial glucose concentration was obtained by using the measured plasma glucose concentration, G_{m}, and by computing the net tissue glucose uptake from the insulin sensitivity index and the plasma insulin concentration. Neglecting non-insulin-dependent uptake, the following equation can be written: (C1) where V_{G} is the glucose distribution volume, Q is the cardiac output, S_{I} is the insulin sensitivity, and I_{p} and I_{p,b} denote the current and the basal insulin concentration in plasma. The time course of G_{a} was thus obtained by computing the derivative Ġ_{m} from a smooth approximation of G_{m}. The parameters were assessed for each subject as follows: V_{G} was obtained by using the estimate of 0.16 l/kg body wt (23), and cardiac output was given by Q = 0.2 × kg body wt^{0.75} (Brody's formula); S_{I} was available from Ref. 32; I_{p,b} was computed as the mean I_{p} value in the last 30 min preceding the clamp. An improved estimate of G_{a} might be obtained by means of a more complex circulatory model (34).

## 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 © 2007 by American Physiological Society