A Mechanistic Model of Intermittent Gastric Emptying and Glucose-Insulin Dynamics following a Meal Containing Milk Components

To support decision-making around diet selection choices to manage glycemia following a meal, a novel mechanistic model of intermittent gastric emptying and plasma glucose-insulin dynamics was developed. Model development was guided by postprandial timecourses of plasma glucose, insulin and the gastric emptying marker acetaminophen in infant calves fed meals of 2 or 4 L milk replacer. Assigning a fast, slow or zero first-order gastric emptying rate to each interval between plasma samples fit acetaminophen curves with prediction errors equal to 9% of the mean observed acetaminophen concentration. Those gastric emptying parameters were applied to glucose appearance in conjunction with minimal models of glucose disposal and insulin dynamics to describe postprandial glycemia and insulinemia. The final model contains 20 parameters, 8 of which can be obtained by direct measurement and 12 by fitting to observations. The minimal model of intestinal glucose delivery contains 2 gastric emptying parameters and a third parameter describing the time lag between emptying and appearance of glucose in plasma. Sensitivity analysis of the aggregate model revealed that gastric emptying rate influences area under the plasma insulin curve but has little effect on area under the plasma glucose curve. This result indicates that pancreatic responsiveness is influenced by gastric emptying rate as a consequence of the quasi-exponential relationship between plasma glucose concentration and pancreatic insulin release. The fitted aggregate model was able to reproduce the multiple postprandial rises and falls in plasma glucose concentration observed in calves consuming a normal-sized meal containing milk components.


Introduction
The mathematical simulation of glucose-insulin dynamics in response to a meal is of great value for decision support related to management of plasma glucose in animals under our care, including domestic species and human patients. The classification of foods according to their effect on the incremental area under the curve of plasma glucose concentrations after ingestion, the so-called glycemic index, was developed to assist in glycemia management [1]. The glycemic response is not just a characteristic of foods, however. Upon consumption of a meal, the ability to dispose of the absorbed glucose, and thus minimize postprandial hyperglycemia and its potentially negative consequences, is dependent on the subject's pancreatic responsiveness, insulin sensitivity and glucose effectiveness. Absorbed glucose has a varying ability to stimulate pancreatic insulin secretion across individuals, with type 1 diabetes being an extreme case in which insulin is only minimally secreted or not at all. The circulating insulin acts on the liver, muscle and adipose tissues to decrease hepatic glucose production and increase peripheral glucose disposal, respectively. In addition, glucose can stimulate its own disposal via a mass-action effect on influx into various tissues. Impairments in any of these three factors of pancreatic responsiveness, insulin sensitivity or glucose effectiveness can exacerbate postprandial hyperglycemia, a hallmark of metabolic syndrome and type 2 diabetes. The minimal glucose and insulin models of Bergman et al. [2] and Toffolo et al. [3] describe all three of these contributions to glucose disposal following an intravenous glucose dose.
Aside from the prominent role of insulin, appearance of glucose in plasma following ingestion of a carbohydrate laden meal is dependent on the rate of gastric emptying. In the early phase of meal consumption, the proximal gastric wall distends and relaxes to accommodate the increased volume. This is quickly followed by tonic contractions of the proximal stomach to push contents towards the distal end where large particles are degraded, aided by the action of peristaltic contractions, and from which small particles and liquids can be pushed into the duodenum through the pyloric sphincter [4]. Liquids display exponential gastric emptying proportional to the gastric volume and without an initial lag, while solids show biphasic gastric emptying with a lag [4]. The intestine responds to early sampling of stomach contents following emptying by altering gut peptide secretion, dependent on the characteristics of the meal. This altered gut peptide profile influences neuro-hormonal control of gastrointestinal function and gastric emptying. Overall, postprandial gastric emptying is a function of meal composition, volume, osmolality and caloric load [4,5].
Equations to simulate gastric emptying have been incorporated into models to interpret plasma and glucose curves following an oral glucose tolerance test [6,7] but we are aware of only one effort to include gastric emptying equations in the simulation of glycemic responses to a normal-sized meal [8,9]. To simulate the pattern of plasma glucose concentrations following a meal, Dalla Man et al. [9] combined a gastric emptying/absorption model of 9 parameters with a glucose-insulin model of 26 parameters, where the 4 components of gastric emptying/ absorption, glucose utilization, endogenous glucose production, and pancreatic insulin secretion were each fitted separately from tracer flux and plasma concentration data. Part of the difficulty with modelling gastric emptying is that outflow from the stomach is intermittent, which produces an erratic timecourse. We have previously used models of 2 to 4 parameters to describe the intermittent appearance of acetaminophen (Ac) in plasma following an oral dose [10]. Acetaminophen is very slowly absorbed from the stomach but rapidly absorbed in the proximal small intestine [11] so its appearance in plasma can be used to estimate gastric emptying rate [10].
In this paper, we present the development and evaluation of a novel mechanistic model that incorporates glucose-insulin dynamics of MINMOD [12] with an intermittent gastric emptying model [10] to predict glycemic responses to a milk-based meal in infant cattle. Prior to development of the ruminant habit, the abomasum of pre-weaned calves functions like the glandular stomach of monogastric animals [13] The carbohydrate in milk is the liquid-associated disaccharide lactose, made up of galactose and glucose moieties. It is hydrolyzed by intestinal lactase (EC 3.2.1.108) and approximately 90% of the absorbed galactose is converted to glycogen in the liver [14,15] while approximately 90% of the absorbed glucose enters the peripheral circulation as free glucose [7,9]. We used the aggregate model to evaluate the importance of gastric emptying to postprandial glycemia, relative to the parameters of insulin release and glucose disposal.

Model Structure Database
Four Holstein-Friesian female calves, housed individually in wheat-straw bedded hutches, at the Trouw Nutrition Ruminant Research facility (Boxmeer, The Netherlands) were maintained on either 4 or 8L of milk replacer (150 g/L dry matter; 24% crude protein, 18% crude fat, and 45.2% lactose on a dry basis; Trouw Nutrition, Deventer, The Netherlands) given twice daily from day 8 to 7 weeks of age. At 4 and 7 weeks of age, calves were provided their morning meal of milk replacer containing 150 mg/kg BW 0.75 Ac as a gastric emptying marker. Blood samples were taken from a jugular vein catheter at -30, 30, 60, 90, 120, 150, 180, 210, 240, 300, 360 and 420 minutes relative to the meal and plasma was analyzed for cAc P , cGl P and cIn P . The -30-min time-point was used as the 0-min time-point in the model (S1 file). The meals of 2 and 4 L contained 136 and 271 g glucose, respectively, producing a range of 1.8 to 5.1 g/kg BW. Postprandial timecourses of plasma concentrations of acetaminophen (cAc P ), glucose (cGl P ) and insulin (cIn P ) in calves at 4 and 7 weeks of age from MacPherson et al. [16] were used to guide model development. Procedures complied with the Dutch Law on Experimental Animals, and the ETS123 (Council of Europe 1985 and the 86/609/EEC Directive) and were approved by the Animal Care and Use Committee from Utrecht University.
The animal datasets are shown in Fig 1. The timecourses of cGl P and cIn P do not all exhibit the ideal behaviour of an oral glucose tolerance test, where there is typically a rapid rise to a peak followed by a sustained decline back to baseline. Rather, the concentrations tend to fluctuate up and down throughout the 420 min, in delayed synchrony with changes in cAc P that reflect gastric emptying rates, and may even sink below baseline concentrations prior to the meal. These erratic behaviours must be accommodated in the simulation model. All model parameters and descriptions are listed in Table 1.

Gastric Emptying
To simulate gastric emptying of Ac in horses, Cant et al. [10] assumed that flow of digesta out of the stomach was either on or off during intervals between successive blood samples and multiplied the first-order rate constant for gastric emptying (k SP ) by a Z-value of either 1 or 0 to indicate whether gastric emptying was occurring (ΔcAc P,i > 0 in sampling interval i) or not (ΔcAc P,i 0), respectively. We used the same approach, where initial mass of Ac in the stomach (iAc S = initial Ac S ) is given by the dose of Ac administered with the meal, and the rate of disappearance to blood follows first order kinetics (Fig 2): Plasma Ac (Ac P ) arises from gastric emptying and disappears according to the first-order elimination constant (k Ac,UAc ): Volume of distribution of Ac P was fixed at 0.9 L/kg BW [10], so that Best-fit parameters of analytical solutions of the gastric emptying equations were estimated with the Solver function of Microsoft 1 Office Excel 1 2007 to minimize residual sums of squares between predicted and observed cAc P . Curve fits were evaluated based on the root mean square prediction error (rMSPE) as a percentage of the mean, calculated as: where pred i is the i-th prediction, obs i is the i-th observation, and n is the number of observations. The model predicted cAc P with an average rMSPE across all 4 calves of 12% of the mean observed cAc P . However, the model did not capture the nuances of apparently intermediate outflows between k SP and 0 (arrows in Fig 3). Such intermediate flows were considered important to reproduce because of the changes in plasma glucose and insulin concentrations with which they were associated (Fig 1). An intermediate rate of Ac P appearance indicates reduced or discontinuous outflow from the stomach during the sampling interval. This intermediate flow was accommodated with a second, non-zero k SP value based on the slope (ΔcAc P, ) between successive time points (i). Thus, k SP = 0 when ΔcAc P,i < -0.05 mg L -1 min -1 , k SP,2 when -0.05 mg L -1 min -1 ΔcAc P,i 0.05 mg L -1 min -1 , and k SP,3 when ΔcAc P,i > 0.05 mg L -1 min -1 . The threshold value of 0.05 mg L -1 min -1 was chosen to minimize rMSPE. The parameter Z in eqs 1 and 2 was replaced with an array of Z-values equal to 0, 1 or 2 to denote k SP = 0, k SP,2 or k SP,3 , respectively, for each sampling interval (Fig 3). The prediction errors decreased from 12 to 9%, on average, when an intermediate gastric outflow was allowed (Fig 3).

Glucose-Insulin Dynamics
Plasma glucose dynamics (Fig 2) are due to exogenous appearance (PGl ex ) and endogenous production (PGl end ) and utilization of glucose (UGl), so that Gastric emptying rates k SP,2 and k SP,3 were used to represent emptying of carbohydrate in the meal from the stomach to the small intestine, as Acetaminophen appearance in plasma during the first 30 min after the meal was typically as rapid as its appearance in the next 30 min (Fig 1), indicating uninterrupted gastric emptying of the meal during the first 60 min. However, plasma glucose concentrations did not typically increase until after the 30-min sample (Fig 1), which could be due to retention of lactose in the stomach or an additional time period for intestinal carbohydrate hydrolysis and absorption. To simulate the delay between the deliveries of acetaminophen and glucose into the plasma, we assigned an absorption time-lag of T lag,SP minutes to the gastric emptying, so that To simulate basal cGl P before the meal, when PGl ex = 0, consideration was given to PGl end . It is known that PGl end decreases as PGl ex increases following a meal, due in part to the effects  Table 1 for parameter descriptions); dashed lines represent stimulatory effects; the dotted line represents inhibitory effect; boxes represent state variables; Ac S = stomach acetaminophen (mg), Ac P = plasma acetaminophen (mg), Gl S = stomach glucose (mmol), Gl P = plasma glucose (mmol), In P = plasma insulin (μg), k SP = gastric emptying rate constant (min -1 ), k Ac,UAc = first-order acetaminophen utilization rate constant (min -1 ), PGl ex = rate of exogenous glucose appearance (mmol min -1 ),), iPGl end = rate of endogenous glucose appearance (mmol min -1 ), PIn = rate of pancreatic insulin release (μg/min), UIn = rate of insulin utilization (μg/min).  of insulin on gluconeogenesis and glycogenolysis. However, the PGl end -suppressing effects of insulin on cGl P are indistinguishable from the UGl-stimulating effects of insulin on cGl P so we chose to represent insulin effects on cGl P via UGl, according to the established equations of the Bergman et al. [2] minimal model, and not via PGl end . Thus, PGl end is maintained at a constant, zero-order flux equal to an initial PGl end (iPGl end ), and UGl is a function of cGl P and the insulin signal (Is) according to a glucose effectiveness constant (k Gl,UGl ) and an insulin sensitivity constant (k Is,UGl ): The Is in MINMOD is a first-order delay of cIn P [12]. This approach is adequate when cIn P exhibits a single peak in the timecourse, but if cIn P rises and falls multiple times following a meal, due to intermittent gastric emptying, the first-order delay removes much of the temporal variation in Is. To retain this temporal variation in Is in the current model, we represented the delay as a time lag of T lag,IS min, so that Prior to T lag,IS , Is(t) is given an initial value (iIs) equal to the baseline cIn P . The differential equation for In (Fig 2) contains the difference between insulin production (PIn) and utilization (UIn): In MINMOD, PIn is a linear function of cGl P above a certain threshold [12]. This structure accommodates, with 2 parameters, the quasi-exponential relation between cGl P and pancreatic insulin release at the lower end of its sigmoidal relationship [17] but it introduces a breakpoint around the threshold cGl P value that causes the first derivative to be discontinuous, which we considered undesirable for continuous simulations that may cross this threshold several times in one run. We chose to simulate sigmoidal kinetics of PIn relative to cGl P continuously with the Hill equation of 3 parameters (V PIn , K Gl,PIn and exp PIn ): UIn remains, as in MINMOD, a mass-action effect of cIn P according to the first-order rate constant k In,UIn : Concentrations of Gl P and In P are calculated assuming a volume of distribution equal to 0.251 L/kg BW, previously estimated in calves [18]: and cIn P ¼

Behaviour and Sensitivity Analyses
Differential equations of the model were written in ACSLX (Aegis Technologies Group, Inc., Orlando, USA) and solved with a 4 th -order Runge-Kutta algorithm using an integration step size of 0.002 min. The final model contains 9 variables and 20 parameters (Table 1, S2 file), 8 of which can be obtained by direct measurement (BW, Z, iAc S , iAc P , iGl S , iGl P , iIn P and iIs, where i represents the initial value) and 12 by fitting to observations (k SP,2 , k SP,3 , k Ac,UAc , k Gl,UGl , k Is, UGl , iPGl end , T lag,SP , K Gl,PIn , k In,UIn , V PIn , exp PIn and T lag,IS ). The model reproduces the multiple rises and falls in cGl P and cIn P that occur in calves during the 420 min following consumption of a normal-sized meal (Figs 4-6). Depending on parameter values, predicted cGl P and cIn P can fall below their respective baseline values prior to the meal, which is an important behaviour to reproduce as it is not uncommon (Fig 1). Values of the parameters can also affect the timing and degree of the cGl P response to changes in gastric emptying and cIn P . In order to understand which characteristics of the postprandial glycemic response are affected by each of the parameters, they were each perturbed above and below reference values to an extent that allowed changes in the cGl P and cIn P curves to be examined. Reference values (Table 1) were set for a 60-kg calf consuming an Ac dose of 3234 mg and carbohydrate load of 396 mmol hexose-equivalents, representing a 2-L meal of milk replacer. The array of Z-values indicating fast, slow or zero gastric emptying was set to a typical pattern of fast initial emptying followed by intermittent, slow gushes. It was assumed that digestibility of lactose was 100% [19], 10% of intestinal galactose entered the circulation as free glucose [14,15] and 10% of intestinal glucose was removed by the splanchnic bed during absorption [7,9]. The remaining parameters of glucose-insulin dynamics were set to generate typical patterns of Ac, Gl and In appearance in plasma over the course of 420 min. Parameter assignments were subject to the constraint that differential eqs 5 and 10 equal zero at t = 0, so that the predicted non-steady, post-prandial state arises from a steady, pre-prandial state. Thus, when V In , K Gl,PIn , or exp In were perturbed in the sensitivity analysis, k In,UIn was set to according to the steady-state constraint and eqs 10 to 14. Likewise, when k Gl,UGl or k Is,UGl were  Mechanistic Model of Gastric Emptying and Glucose-Insulin Dynamics perturbed, the basal steady-state of Gl P was achieved by setting according to eqs 5, 8, 9, 13 and 14. The predicted reference postprandial patterns of cGl P and cIn P each exhibited 3 peaks in association with the gastric emptying profile (Figs 4 to 6). Altering T lag,IS , K Gl,PIn and exp PIn changed the number and amplitude of peaks in cGl P and cIn P , while k SP , k Gl,UGl , k Is,UGl and V PIn only affected the amplitude of peaks (Figs 4 and 6). Of these latter parameters, k Gl,UGl and V PIn exerted relatively small effects on the amplitude of cGl P peaks. In addition to affecting wave amplitudes, T lag,IS , k Is,UGl and V PIn shifted the times at which cGl P or cIn P peaks occurred. T lag,SP also affected peak times without altering other characteristics of the glycemic response (Fig 5). According to these simulations, each of the parameters exerted unique effects on the postprandial responses, except for similarities between K Gl,PIn and exp PIn , and k Is,UGl and V PIn .
Gastric emptying is under neural and hormonal controls that regulate delivery of nutrients to the periphery for metabolism. How much of the glycemic response to a meal is due to gastric emptying rate versus pancreatic responsiveness and tissue insulin sensitivity remains an interesting question for the maintenance of normoglycemia and the attendant therapeutic implications. Because the model contains all three of these control elements, we used it to evaluate the relative role of each in the glycemic response, as indicated by areas under 420-min cGl P and cIn P curves (AUC Gl and AUC In , respectively). Results of changing each parameter from 0.5 to 1.5X its reference value are presented as a sensitivity coefficient (SC A ), equal to the fractional change in AUC A relative to the fractional change in parameter value: Parameters that affected AUC Gl the most ( Table 2) were those related to pancreatic response (K Gl,PIn and exp PIn ), followed by insulin sensitivity (k Is,UGl , T lag,IS ) and then gastric emptying (k SP ). The sensitivity to K Gl,PIn was 50X that to k SP . However, sensitivity of AUC In was 20X greater than sensitivity of AUC Gl to k SP (Table 2). Pancreatic responsiveness and insulin sensitivity parameters still exhibited stronger effects than the gastric emptying parameter on AUC In , but only by 2X instead of 50X. The difference between insulin and glucose responses to k SP is interesting because one might intuitively presume that both would respond in a similar fashion, given the reciprocal nature of the control paradigm, in which cGl P affects cIn P , and cIn P affects cGl P . A large effect of k SP on AUC In , with a much smaller effect on AUC Gl , suggests that pancreatic responsiveness is stimulated by faster gastric emptying, so that rapid insulin release (high AUC In ) prevents hyperglycemia, leading to low AUC Gl . When gastric emptying was delayed in humans by the amylin analog pramlintide, AUC In over the first 120 min decreased significantly while AUC Gl over the same time period was affected little [20], similar to our simulation results. This decrease was not accompanied by a large change in pancreatic responsiveness parameters [20]. Our simulations show that the pancreatic response to k SP can be independent of changes in the parameters of pancreatic responsiveness V In , K Gl,PIn and exp PIn , because we found an increase in AUC In without altering those parameters. The cause of the higher AUC In as k SP increased in our simulations was the quasi-exponential nature of the effect of cGl P on insulin release, where small changes in cGl P exert larger effects on PIn when cGl P is at high values compared to low. When k SP was varied in the sensitivity analysis, amplitudes of peaks in cGl P and cIn P curves were affected but the number and timing of peaks did not change. Slower gastric emptying might be expected to delay the time to peak cGl P and cIn P , as in the pramlintide experiment of Hinshaw et al. [20]. When gastric emptying was set at a continuous rate, by setting all instances of the parameter array Z = 2, then a slower k SP prolonged the time to peak cIn P (data not shown). Thus, time to peak is a consequence of the rate of gastric emptying (kSP × cGlS), and the number and sequence of fast versus slow or negligible gastric emptying bouts, which we did not present because of the large number of permutations possible.

Parameter estimation
Estimation of k SP,2 , k SP,3 and k Ac,UAc from the 12 cAc P data points has already been described. The remaining parameters of glucose-insulin dynamics (k Gl,UGl , k In,UGl , T lag,SP and T lag,IS in eqs 7 to 9, and V PIn , K Gl,PIn , exp PIn , and iIn P in eqs 11, 12 and 14) were estimated with a differential evolution algorithm [21] to minimize the sum of residual sums of squares between predicted and observed cGl P , and predicted and observed cIn P . The differential evolution procedure is designed to search for optimal solutions within a large parameter space so that global rather than local minima in RSS are reached. A steady basal state was enforced with eqs 15 and 16. The evolutionary algorithm was run with 80 sets of parameter values in each generation. After 200 generations, means and standard errors of each parameter in the 30 best-fit sets were calculated.
Parameter estimates and their standard errors for each set of animal data are presented in S1 table. According to standard errors, parameter values were significantly different from 0 and were highly identifiable. The fitted aggregate model was able to reproduce postprandial glycemia in the 4 test subjects following the consumption of a milk-based meal with rMSPE percentages ranging from 5 to 15% for glucose (Fig 7) and from 19 to 58% for insulin (Fig 8). Curves with multiple rises and falls in cGlP and cInP were not worse fit than less erratic curves.
Entry rates of glucose into the plasma from intestinal lactose were based on fixed parameters of 100% absorption, 10% conversion of intestinal galactose to GlP, and 90% conversion of intestinal glucose to Gl P . To test the effect of these parameters on model outputs, we assessed fits to the 4 test subjects with 50% higher iGl S values than what the fixed parameters predicted. Percentage rMSPE for cGl P curves were 9.1% on average with the original iGl S values and 9.5% with the higher iGl S values, indicating little effect on curve fits.

Discussion
We have presented a novel model for the simulation of postprandial glycemia based on intermittent, exponential gastric emptying with a time lag for appearance of intestinal glucose in plasma, and the minimal models of glucose and insulin dynamics of MINMOD [12]. Dalla Man et al. [8] accommodated intermittent gastric emptying in a model of postprandial glycemia as a convolution of three pulses with a second-order decay function for a total of 6 parameters. Other glycemia models with gastric emptying have assumed continuous gastric outflow, where the delay between emptying and appearance in plasma is acommodated with an intermediate intestinal pool from which absorption proceeds according to a first-order rate constant [7,9]. With 2 k SP parameters and a lag-time (T lag,SP ), our intermittent gastric emptying model of 3 parameters can be considered a minimal model. We previously developed and assessed mathematical models of gastric emptying using data from horses given various dietary treatments containing Ac [10]. The best fit was achieved by a first-order model of plasma Ac appearance with parameters to describe the duration of periodic gushing of gastric contents as well as the quiescence between gushes. However, the use of two non-zero rate constants to describe gastric emptying rate in the current model provided a simpler solution for improvement of the prediction accuracy of Ac P appearance. A comparison of Ac and Gl profiles in plasma indicated there was a delay between appearance of the liquid marker in the meal and appearance of glucose. Casein in the meal forms a clot in the stomach that slows its rate of emptying into the small intestine [22,23]. It is possible that the clot retains some of the liquid-associated lactose so that lactose entry into the small intestine is also slowed. However, cumulative bihourly samples of duodenal digesta collected from calves fed clotting or non-clotting milk replacers were not different in lactose content, indicating no effect of the clot on lactose emptying [22]. More frequent samples, as in our datasets [16], may be subject to a clot effect because the reference T lag,SP values we found were less than 30 min. The contribution of lactose hydrolysis to the delay time does not appear to be significant because appearances of 13 C label from intact and hydolyzed lactose in a meal consumed by lactose-tolerant humans were identical [24]. This leaves glucose absorption as a potential cause of the delay.
As a percentage of mean observations, prediction errors for cIn P were higher than for cAc P or cGl P . Part of the reason for higher relative errors is the large fold-changes in cIn P that occur during the 420-min timecourse, compared to fold-changes in cAc P and cGl P . However, the lower goodness of insulin fits suggests that improvements are possible in the insulin simulations, possibly through consideration of incretin dynamics. Although the model was developed using datasets from calves fed milk, it has utility for predicting post-meal glucose and insulin kinetics in other animal models as well as human subjects.

Conclusion
The combination of minimal models of gastric emptying, plasma glucose dynamics and plasma insulin dynamics suitably describes the erratic postprandial glycemia following a milk-based meal, including depressions below the baseline cGl P prior to the meal. Sensitivity analysis of the model indicates that faster gastric emptying increases pancreatic responsiveness and keeps plasma glucose concentrations low, independent of the parameters of pancreatic response, simply through the quasi-exponential nature of the relationship between cGl P and pancreatic insulin release. The model has the potential to be used in the evaluation of dietary treatments for their net effects on pancreatic responsiveness, insulin sensitivity and glucose effectiveness, as well as to predict glycemic responses to various normal-sized meals. Supporting Information S1 File. Postprandial plasma acetaminophen (mg/L), glucose (mM) and insulin (ng/ml) of four datasets used for model development.
(XLSX) S2 File. The acslX program with reference animal parameter values. (PDF) S1 Table. Parameters values of four test datasets. Best-fit parameters are shown as mean ± standard error. (PDF)