Building a Glucose-Insulin Model for Mixed-Meal Data

From a Mixed Meal to Concentration Predictions

A mixed-meal test produces linked time courses: carbohydrate leaves the stomach, glucose appears in the circulation, plasma glucose changes, and beta cells respond by secreting insulin. A quantitative model makes those links explicit. It also makes clear which observations inform each process and which assumptions are needed before a predicted curve can be interpreted.

The model below is deliberately compact. It describes meal-derived glucose appearance, plasma glucose concentration, and plasma insulin concentration with coupled ordinary differential equations. It is a useful starting model for data analysis, not a complete account of meal physiology: it uses a fixed endogenous glucose input and a lumped glucose-removal term, and it does not separately represent incretin effects or the nutrient effects of meal fat and protein.

Definition: Let \(M(t)\) be the meal-derived glucose-equivalent amount remaining in the stomach and \(A(t)\) the amount in an intestinal absorption compartment, both expressed in mg/dL-equivalent. Let \(G(t)\) be plasma glucose in mg/dL and \(I(t)\) plasma insulin in \(\mu\text{U}/\text{mL}\). The meal-derived glucose appearance rate is \(R_a(t)\), in mg/dL per hour. The model parameters and initial conditions are defined below.

Meal Compartments and the Appearance Rate

The meal is represented by two sequential first-order processes. Gastric emptying transfers material from \(M\) to \(A\) at rate \(k_gM\); intestinal absorption transfers material from \(A\) into the circulation. The fraction \(F\) of the meal-equivalent amount that reaches the circulation is the bioavailability fraction, with \(0\leq F\leq1\). Set \(M(0)=D\) and \(A(0)=0\), where \(D\) is the ingested carbohydrate amount expressed in the model's concentration-equivalent units.

\frac{dM}{dt}=-k_gM,\qquad \frac{dA}{dt}=k_gM-k_aA,\qquad R_a(t)=F k_aA(t).

Here \(k_g\) and \(k_a\) have units of inverse hours. \(R_a(t)\) is an instantaneous rate, whereas \(FD\) is an amount. They must not be compared as though they had the same units. The separate physical conditions are \(R_a(t)\geq0\) and \(\int_0^t R_a(s)\,ds\leq FD\): cumulative appearance by time \(t\) cannot exceed the bioavailable meal amount.

The compartment equations give a useful closed form. For \(k_a\ne k_g\), solving the first equation gives \(M(t)=De^{-k_gt}\). Substituting this into the equation for \(A\), using the initial condition \(A(0)=0\), gives the expressions below.

A(t)=\frac{Dk_g}{k_a-k_g}\left(e^{-k_gt}-e^{-k_at}\right),\qquad R_a(t)=\frac{FDk_gk_a}{k_a-k_g}\left(e^{-k_gt}-e^{-k_at}\right).

For positive \(k_g\) and \(k_a\), the difference of exponentials has the sign needed to make \(A(t)\) and \(R_a(t)\) nonnegative, regardless of which rate constant is larger. The limiting case \(k_a=k_g=k\) is obtained by taking the limit: \(A(t)=Dkt e^{-kt}\). To verify the cumulative constraint, integrate the appearance rate, or use the compartment balance \(d(M+A)/dt=-k_aA\). Since \(M+A\) is nonnegative and starts at \(D\), cumulative absorption cannot exceed \(D\), and cumulative appearance, which is \(F\) times cumulative absorption, cannot exceed \(FD\).

Worked Example: Appearance Rate and Cumulative Appearance

Suppose a hypothetical meal has \(D=40\) mg/dL-equivalent, \(F=0.80\), \(k_g=0.50\ \text{h}^{-1}\), and \(k_a=1.00\ \text{h}^{-1}\). At one hour, calculate the instantaneous appearance rate and cumulative appearance.

The appearance-rate formula becomes \(R_a(t)=32(e^{-0.50t}-e^{-1.00t})\) mg/dL per hour, because \(FDk_gk_a/(k_a-k_g)=0.80(40)(0.50)(1.00)/(0.50)=32\). Therefore, \(R_a(1)=32(0.6065-0.3679)\approx7.637\) mg/dL per hour.

Integrating the rate from zero to one hour gives \(\int_0^1R_a(s)\,ds=32[(1-e^{-0.50})/0.50-(1-e^{-1.00})]=4.95\) mg/dL-equivalent, to three significant figures. This cumulative amount is below \(FD=0.80(40)=32\) mg/dL-equivalent. The comparison is dimensionally appropriate: both quantities are amounts. The rate at one hour, by contrast, has units per hour.

Coupling Meal Appearance to Glucose and Insulin

To represent the glucose response, let \(G_b\) and \(I_b\) be the fasting baseline concentrations. The glucose equation includes meal appearance, a constant endogenous input \(P_0\), and two lumped removal terms: glucose-dependent removal \(k_GG\) and insulin-associated removal \(S_II\). The constant input is chosen so that the fasting baseline is an equilibrium. Thus \(P_0=k_GG_b+S_II_b\). This is a bookkeeping device for the baseline balance, not a claim that endogenous glucose production is truly constant during a meal.

\frac{dG}{dt}=R_a(t)+P_0-k_GG-S_II,\qquad P_0=k_GG_b+S_II_b.

The units are \(k_G\) in \(\text{h}^{-1}\), \(S_I\) in \((\text{mg}/\text{dL}/\text{h})/(\mu\text{U}/\text{mL})\), and \(P_0\) in mg/dL per hour. At \(G=G_b\), \(I=I_b\), and \(R_a=0\), the equation gives \(dG/dt=0\). The insulin term represents one aggregate insulin-associated effect on glucose removal; adding another term for the same effect would change the model rather than merely make it more detailed.

A minimal beta-cell and insulin-clearance model uses a basal secretion rate plus secretion that increases with glucose above baseline. Insulin is removed by first-order clearance. This simple relationship is a phenomenological approximation: it does not separately describe secretion driven by incretins, amino acids, or other meal components.

\frac{dI}{dt}=S_0+\beta\max(G-G_b,0)-k_eI,\qquad S_0=k_eI_b.

Here \(S_0\) is basal insulin secretion in \((\mu\text{U}/\text{mL})/\text{h}\), \(\beta\) is the glucose-stimulated secretion coefficient in \((\mu\text{U}/\text{mL}/\text{h})/(\text{mg}/\text{dL})\), and \(k_e\) is the insulin clearance rate constant in \(\text{h}^{-1}\). The positive-part function \(\max(G-G_b,0)\) prevents this particular meal-stimulation term from becoming negative when glucose is below baseline. At baseline, \(S_0-k_eI_b=0\), so insulin is also at equilibrium.

Assumptions: The meal compartments are initially empty apart from \(M(0)=D\); transfer and absorption follow constant first-order rate constants; \(F\) is constant; and the glucose and insulin parameters do not change over the simulated interval. The model treats its glucose-equivalent meal amount and plasma concentration as commensurate after scaling. It assumes a fixed endogenous glucose input, a lumped removal relationship, and glucose-driven secretion above \(G_b\). These assumptions should be checked against the intended use and the observed time range.

Worked Prediction of Glucose and Insulin Time Courses

For illustration, use the meal parameters from the first example and set \(G_b=90\) mg/dL, \(I_b=10\ \mu\text{U}/\text{mL}\), \(k_G=0.08\ \text{h}^{-1}\), \(S_I=0.03\ (\text{mg}/\text{dL}/\text{h})/(\mu\text{U}/\text{mL})\), \(k_e=0.50\ \text{h}^{-1}\), and \(\beta=0.04\ (\mu\text{U}/\text{mL}/\text{h})/(\text{mg}/\text{dL})\). Basal balance gives \(P_0=0.08(90)+0.03(10)=7.50\) mg/dL per hour and \(S_0=0.50(10)=5.00\ \mu\text{U}/\text{mL}/\text{h}\). Begin at \(G(0)=90\), \(I(0)=10\), \(M(0)=40\), and \(A(0)=0\).

To show the calculation transparently, use forward Euler integration with step \(h=0.50\) h: each state at the next time is its current value plus \(h\) times its current derivative. These values are approximate numerical predictions, not exact solutions of the coupled glucose-insulin equations.

Worked Example: One Euler Step and a Short Predicted Course

At time zero, \(dM/dt=-0.50(40)=-20\) mg/dL-equivalent per hour and \(dA/dt=0.50(40)-1.00(0)=20\) mg/dL-equivalent per hour. The initial glucose and insulin derivatives are both zero because the system is at baseline and \(R_a(0)=0\). After \(h=0.50\) h, the Euler updates give \(M=40+0.50(-20)=30\), \(A=0+0.50(20)=10\), \(G=90\), and \(I=10\).

At 0.50 h, \(R_a=0.80(1.00)(10)=8.00\) mg/dL per hour. The glucose derivative is \(8.00+7.50-0.08(90)-0.03(10)=8.00\) mg/dL per hour, and the insulin derivative is \(5.00+0.04\max(90-90,0)-0.50(10)=0\). The next Euler update gives \(G(1.00)=90+0.50(8.00)=94.00\) mg/dL and \(I(1.00)=10.00\ \mu\text{U}/\text{mL}\).

Continuing the same update gives \(M=22.50\), \(A=12.50\), \(G=94.00\), and \(I=10.00\) at 1.00 h; then \(M=16.875\), \(A=11.875\), \(G=98.84\), and \(I=10.08\) at 1.50 h. At 2.00 h the predicted concentrations are approximately \(G=103.235\) mg/dL and \(I=10.237\ \mu\text{U}/\text{mL}\). For example, at 1.00 h the current appearance rate is \(0.80(12.50)=10.00\) mg/dL per hour; the glucose derivative is \(10.00+7.50-0.08(94)-0.03(10)=9.68\) mg/dL per hour, and the insulin derivative is \(5.00+0.04(94-90)-0.50(10)=0.16\ \mu\text{U}/\text{mL}/\text{h}\). These derivatives produce the stated 1.50 h values after multiplication by 0.50 h. At 1.50 h, \(R_a=0.80(11.875)=9.50\), \(dG/dt=9.50+7.50-0.08(98.84)-0.03(10.08)=8.7904\), and \(dI/dt=5.00+0.04(98.84-90)-0.50(10.08)=0.3136\), giving \(G=103.2352\) and \(I=10.2368\) at 2.00 h.

The rising glucose concentration stimulates insulin secretion, while clearance removes insulin continuously. These illustrative values depend on the chosen parameters and the relatively coarse step size. A smaller step or an adaptive numerical solver should be used to check that predictions are not artifacts of the integration step.

Comparing Predictions with Mixed-Meal Data

Observed measurements are samples of the model's continuous-time concentration curves. A basic observation model is \(G_{\text{obs},j}=G(t_j)+\epsilon_{G,j}\) and \(I_{\text{obs},j}=I(t_j)+\epsilon_{I,j}\), where measurement errors have mean zero under the assumed observation model. The errors may differ in scale between glucose and insulin; residuals should therefore be interpreted in their own units or appropriately scaled. The ODEs specify the predicted time courses, while an estimation procedure would choose parameters to bring those predictions into agreement with measurements.

Worked Example: Residuals at a Sampling Time

Suppose a hypothetical mixed-meal data set reports glucose of 101.0 mg/dL and insulin of \(10.40\ \mu\text{U}/\text{mL}\) at 1.00 h. The Euler predictions above are 94.00 mg/dL and \(10.00\ \mu\text{U}/\text{mL}\). Define residual as observed minus predicted. The glucose residual is \(101.0-94.00=7.00\) mg/dL, and the insulin residual is \(10.40-10.00=0.40\ \mu\text{U}/\text{mL}\).

If, for illustration, the assumed measurement-error standard deviations are 3.00 mg/dL for glucose and \(1.00\ \mu\text{U}/\text{mL}\) for insulin, the sum of squared standardized residuals for these two observations is \((7.00/3.00)^2+(0.40/1.00)^2=5.6044\). This is a measure of discrepancy under those assumed error scales; it is not by itself a test of model adequacy or evidence that any particular parameter is correct. The example also illustrates why residuals should not be compared directly across analytes with different units and variability.

Diagnostics, Alternatives, and Limits

A prediction should be checked as a mathematical object and against the data. First check units in every differential equation: each term on the right-hand side must be a concentration rate for \(G\) or \(I\), or a compartment amount rate for \(M\) or \(A\). Next check the meal balance: \(R_a(t)\) should be nonnegative, and its integral should not exceed \(FD\). Check baseline equilibrium by setting \(R_a=0\), \(G=G_b\), and \(I=I_b\). For numerical integration, repeat the calculation with a smaller step and compare the predicted curves; substantial changes indicate that the original step was too large.

In data, systematic residual patterns are informative. Residuals that are mostly positive around the glucose peak may indicate that the model underpredicts the peak, while consistently late predictions may point to an inadequate timing description. However, several parameter combinations can produce similar curves. A meal appearance profile and a glucose removal parameter, for example, can compensate for one another over a sparse sampling schedule. Insulin observations help constrain secretion and clearance, but do not automatically identify them separately without adequate timing and informative data.

A two-compartment gastric-absorption model is more mechanistic than prescribing an arbitrary appearance curve, but it still simplifies meal processing. An empirical appearance function can be useful when the measured data do not support separate gastric and intestinal rate constants. Similarly, forward Euler is transparent for hand calculations, whereas higher-order or adaptive solvers are generally preferable for fitting data and checking numerical accuracy. No solver can repair a structural model that omits a process essential to the question.

The equations should not be extrapolated without scrutiny to extreme glucose or insulin values, meals with substantially different composition, or conditions where endogenous production and insulin sensitivity change over the observation window. The fixed \(P_0\), linear insulin-associated removal, and thresholded linear secretion are approximations. If residual patterns, physiology, or the study question demand more detail, revise the structure and justify the added components rather than interpreting every fitted parameter as a direct physiological measurement.

Common pitfall: Saying that the instantaneous rate \(R_a(t)\) cannot exceed the total available meal amount compares mg/dL per hour with mg/dL and is dimensionally invalid. The correct separate statements are \(R_a(t)\geq0\) and \(\int_0^tR_a(s)\,ds\leq FD\), the latter comparing cumulative amount with available amount.

Check Your Understanding

Answer each question, then check the concise solutions below.

What units should \(R_a(t)\) and \(\int_0^tR_a(s)\,ds\) have, and which can be compared with \(FD\)?
If \(G=G_b\), \(I=I_b\), and \(R_a=0\), show that the glucose derivative is zero.
For \(D=30\) mg/dL-equivalent and \(F=0.70\), what is the maximum cumulative meal-derived appearance allowed by this model?
A model predicts glucose of 96 mg/dL when the observed value is 100 mg/dL. What is the residual under the observed-minus-predicted convention, including units?

Solutions. (1) \(R_a(t)\) has units mg/dL per hour; its integral has units mg/dL and is the quantity compared with \(FD\). (2) Substitution gives \(0+k_GG_b+S_II_b-k_GG_b-S_II_b=0\). (3) The bound is \(FD=0.70(30)=21\) mg/dL-equivalent. (4) The residual is \(100-96=4\) mg/dL, positive because the observation exceeds the prediction.