## Abstract

Noninvasive breath tests for gastric emptying are important techniques for understanding the changes in gastric motility that occur in disease or in response to drugs. Mice are often used as an animal model; however, the gamma variate model currently used for data analysis does not always fit the data appropriately. The aim of this study was to determine appropriate mathematical models to better fit mouse gastric emptying data including when two peaks are present in the gastric emptying curve. We fitted 175 gastric emptying data sets with two standard models (gamma variate and power exponential), with a gamma variate model that includes stretched exponential and with a proposed two-component model. The appropriateness of the fit was assessed by the Akaike Information Criterion. We found that extension of the gamma variate model to include a stretched exponential improves the fit, which allows for a better estimation of *T*_{1/2} and *T*_{lag}. When two distinct peaks in gastric emptying are present, a two-component model is required for the most appropriate fit. We conclude that use of a stretched exponential gamma variate model and when appropriate a two-component model will result in a better estimate of physiologically relevant parameters when analyzing mouse gastric emptying data.

- motility
- smooth muscle
- gastrointestinal tract
- empirical mathematical models

the measurement of gastric emptying in human subjects and animal models of disease is an important technique for understanding the changes in gastric motility that occur under normal physiological conditions, in disease, or in response to drugs. Acute physiological changes in gastric emptying can be identified in response to the content and nature of a meal including the effects of calories, osmolarity, lipid, protein, and carbohydrate content (9). Determination of gastric emptying has clinical relevance since the delay of gastric emptying of solids is required for positive diagnosis of gastroparesis. Gastroparesis can occur secondary to drugs or systemic disease as a complication of diabetes or remain idiopathic (21). Accelerated gastric emptying is also observed in subsets of patients with symptoms suggestive of functional dyspepsia, or dumping, as well as in a subset of patients with diabetes particularly early after diagnosis (12, 25). Accelerated gastric emptying is also observed in animal models of diabetes (7). Prokinetic agents have been identified and studied by using gastric emptying tests in experimental animals and candidate drugs have been tested for unwanted effects on gastric emptying (24).

Several noninvasive tests have been developed to follow gastric emptying in humans including scintigraphy and breath tests (5, 29), and some of these tests have been applied for animal studies (3, 6, 10, 28). The ability to do repeated measurements in the same, unrestrained, diabetic mouse has been useful to identify the factors that contribute to the development of diabetic gastroparesis and to show that treatment of mice can reverse the symptoms and underlying pathological changes (7, 8, 15).

Analysis of the data from gastric emptying tests has typically used models that gave a readout of the time to half emptying (*T*_{1/2}) value or reported the proportion of a meal that remained in the stomach at fixed time points. Time points at 1, 2, and 4 h after feeding can be effectively used to screen for fast or slow gastric emptying using scintigraphy (29). A 4-h time point is required to make an accurate diagnosis of delayed gastric emptying. Breath tests allow samples to be collected at multiple time points with detailed kinetic information available. These breath tests depend on the principle that the rate of appearance of ingested isotope-enriched CO_{2} in the breath is limited by the rate of gastric emptying. Assays based on detection of radioactive ^{14}C or stable, spectroscopically detectable, ^{13}C have been developed (11). For example, we use a system that measures ^{13}CO_{2} from air samples at a rate of 1 Hz for at least 4 h and that automatically switches sampling between up to 12 mouse chambers so that individual chambers are sampled at every 4–6 min (10). We have typically applied a gamma variate model to these data. This model was first applied to gastric emptying in mice by Symonds et al. (28) to obtain *T*_{1/2} (time to half emptying) values for emptying of a solid meal. Our approach was originally developed for human solid gastric emptying studies where the gamma variate model works relatively well for fitting the data. However, when we use our system that tracks gastric emptying at high time and signal resolution, we have observed that the time course of gastric emptying is not well described by a single gamma variate and that occasionally there is evidence for more than one phase to gastric emptying, as previously reported by Symonds et al. (28); the additional phase showing up as a second peak in the plotted data. The existing models were relatively weak with respect to fitting these two important characteristics of the data, namely the fit to the later time points in the study and the handling of data with two peaks of ^{13}CO_{2} in the breath.

The outputs of the gamma variate model fit are the *T*_{1/2} value, the gastric emptying coefficient (GEC), which is a measure of the overall gastric emptying rate, and the *T*_{lag} value (11, 13). In the literature two definitions of *T*_{lag} have been used, the time from eating until the first food particles have been triturated to a size to pass into the duodenum, and the time to the peak rate of gastric emptying, which correlates with the time where all the food has been triturated to a size where it can be emptied (32). We have not detected in mice a clear lag between feeding and the appearance of content in the duodenum. However, *T*_{lag}, defined as the time to peak emptying, has been reported in some studies on gastric emptying in mice (28).

The accuracy and reproducibility of the *T*_{1/2} values are important for identifying changes in gastric emptying that are physiologically or clinically significant. Investigators carry out costly, longitudinal studies on mice, and the goal was to improve the power and reliability of gastric emptying studies to make more accurate decisions on disease pathophysiology and effectiveness of potential therapies. The ability to fit the data to two peaks is important because although this phenomenon has been previously reported (28), none of the currently used models take it into account and therefore cannot be used to test the underlying physiological basis for the double peaks. Therefore, a key objective of developing new models is to allow studies to elucidate the physiological or pathological factors that affect the presence and timing of the second peak.

To accommodate gastric emptying patterns that do not fit to simple models and to allow investigation of how these variations in gastric emptying might reflect variations in gastric motility, we proposed and tested a gamma variate model with stretched exponential and a two-component gamma variate model with stretched exponentials, designed to describe data with two peaks. From the best fits by these models we calculate *T*_{1/2}, and the time to the peak rate of emptying (*T*_{lag}) (32), which is a variable that is not always well estimated by the conventional gamma variate model. The gastric emptying coefficient is well defined for a single gamma variate model, but it becomes difficult to define it for a two-component model. Therefore, to characterize the gastric emptying rate, we estimate the maximal value for the rate of ^{13}CO_{2} production (*R*).

To evaluate the merits of the proposed models, we fitted 175 gastric emptying data sets with two standard models (gamma variate and power exponential), with a gamma variate model that includes stretched exponential, and with a proposed two-component model. The appropriateness of these models in fitting the data with respect to the quality of the fit and the number of free parameters introduced was assessed by Akaike Information Criterion (AIC) (14, 17, 27). Our results suggest that the extension of the gamma variate model to include stretched exponential improves the fit, which allows better estimation of characteristic parameters (*T*_{1/2}, *T*_{lag}, *R*). Also, it became evident that the two-component model is required for an appropriate fit in some cases.

## MATERIALS AND METHODS

#### Experimental methods and data preprocessing.

The gastric emptying tests were done as previously reported (6) using female NOD/LtJ mice (Jackson Laboratories, Bar Harbor, ME). Protocols were approved by the Mayo Clinic Institutional Animal Care and Use Committee. Mice were habituated to the testing chambers and received egg meals for 2 wk before testing. Mice were fasted for 12–18 h with free access to water before testing. Metal mesh flooring was inserted into the cages to prevent coprophagia. At the time of testing, mice are placed in chambers that allow the freedom to turn around and the chambers are connected to air flow that originates in a zero CO_{2} air generator. The flow rate of each chamber is controlled by individual flow meters and adjusted to be between 0.5 and 1.5 l/min so that exhaled breath samples measured between 1,500 and 2,000 ppm CO_{2}. After a baseline reading, mice were fed 200 mg cooked egg yolk containing 2.5 μmol ^{13}[C]octanoic acid. Because the mice were fasted and trained, they generally ate the test meal within 2 min. A computer-controlled solenoid valve system, called the multiple input unit, controlled which animal chamber is sampled and analyzed by a laser spectrometer CO_{2} analyzer (Los Gatos Research, Mountain View, CA), which measured the ^{13}CO_{2}, ^{12}CO_{2}, and H_{2}O content of the breath. The instrument utilizes off-axis integrated cavity output spectroscopy to mitigate the need to resolve CO_{2} spectra via traditional laser 7 resonance spectroscopy. The concentrations of ^{13}CO_{2}, ^{12}CO_{2}, and H_{2}O are calculated via integration of defined peaks within the spectra generated by the rapidly tunable laser.

Samples containing exhaled breath were analyzed at a rate of 1 Hz and the valves switched between up to 12 mouse chambers every 25–45 s so that each chamber was sampled at least every 4–6 min for up to 4 h. It means that the system collects a rapid series of data points from say chamber 1 for a period of time and then the valves switch to collect data from the adjacent chambers in turn before returning to chamber 1. This implies that for a given mouse there are groups of data points resolved on a time scale of seconds, and those groups are resolved on a time scale of 4–6 min. The enrichment of ^{13}CO_{2} to ^{12}CO_{2} was determined and expressed as a percent change in ^{13}C calculated from the sample.

The measurements from the first 10 s of each time group of data were removed, because the preliminary experiments using the multiple input unit to switch between standard tanks of 500 ppm CO_{2} to 2,000 ppm CO_{2} indicated that 10 s was sufficient to completely flush the previous sample. In the remaining 15–35 data points from a given time group, there still can be outliers due to misinterpretation of the spectrogram that is usually the consequence of laser alignment. Based on inspection, we chose to remove at most 10 outliers per time group using the generalized extreme studentized deviate many-outlier procedure (22).

After the outliers were removed, data were adjusted for the baseline, which is the level of ^{13}CO_{2} that is normally exhaled by a given mouse. The baseline is approximated by the average of the measurements recorded before the egg is administered. It is noteworthy that the baseline value is much larger (at the level of ∼1.1) than the signal specific to enhancement of ^{13}CO_{2}, which is in the range [0,0.04]. This implies that the small errors in the estimation of the baseline may influence the analysis. Therefore, we introduced an additive constant in the model as a correction to the baseline, which is then determined by fitting. Data sets from mice with accelerated, normal, and delayed gastric emptying values were analyzed using the four empirical models described below.

#### Empirical models.

In previous studies fitting data for ^{13}CO_{2} enrichment in the breath vs. time was accomplished by using gamma variate function (6, 11, 28). According to this model, the rate *y*(*t*) of *C*^{13} excretion in breath is written in the form:
(1)
where *a*, *b*, and τ are positive constants and *t* is time. The function exhibits one peak with *y*(0) = *y*(∞) = 0. The gamma variate function can be considered useful as a simple approximation of the linear combinations of exponential functions that occur in linear multi-compartment models. The data from high-resolution gastric emptying measurements are a reflection of more than one compartment in the complex physiological processes underlying ^{13}C kinetics. In the following we will refer to the model given by *Eq. 1* as the gamma variate (GV) model.

Another empirical model used to analyze ^{13}CO_{2} data (11, 28, 23) can be written in the form:
(2)
where *a*, *k* > 0, and β > 1. This is also unimodal function with *y*(0) = *y*(∞) = 0. In the following, we will refer to this model as power exponential (PE) model according to Ref. 23.

In this article we propose a simple generalization of the gamma variate model that includes a stretched exponential function:
(3)
Here *a*, *b*, τ, *c* > 0, and the function is also unimodal with *y*(0) = *y*(∞) = 0. This model is motivated by the idea that by generalizing the gamma variate function to include the stretched exponential, one can only better approximate the multiexponential function related to multicompartment models. On the other hand the stretched exponential contains the power function (*t*/τ)^{c}, which is a signature of fractal kinetics (1, 2, 26). It is likely that fractal kinetics plays a role in physiological processes involved, as they are happening in spaces exhibiting fractal structures. In our experience the fits to gastric emptying data by the model given by *Eq. 3* are almost always preferred by the AIC for model selection (14, 17, 27) compared with the models given by *Eqs. 1* and *2*. In the following we will refer to this model as the gamma variate model with stretched exponential (GVS).

In a considerable number of cases, gastric emptying data exhibited patterns of visually distinguished two peaks or one peak and a shoulder (in the present study that was 19.4% of cases). This was the motivation to introduce a model capable of describing two peaks, which is not the case for models given by *Eqs. 1–3*. The model we propose is the linear combination of two gamma variates with stretched exponentials, which is a phenomenological model sufficiently flexible to describe either two peaks or a peak and a shoulder. It is given by:
(4)
where *b*_{j}, τ_{j}, *c*_{j} > 0. In principle *a*_{1} or *a*_{2} can achieve a negative value as a result of fitting the model to the data. In these cases it is possible that in some time intervals the best-fit curve may achieve, usually small, negative values, while still approaching zero for infinite times. Since a gastric emptying curve by definition is nonnegative, we define our model by:
(5)
In the following we will refer to this model as the two-component gamma variate model with stretched exponentials (GVS2).

To compare the adequacy of the models considered (*Eqs. 1*, *2*, *3*, and *5*), we performed least-squares fits to the data by each of them and then used the AIC (14, 17, 27) to determine which should be selected as preferable. This criterion takes into account how good the fit is and how many free parameters were used. Assuming that the errors in data points are distributed according to the same normal distribution, which is reasonable for measurements of ^{13}CO_{2}, the corrected Akaike Information Criterion (AIC_{c}) is based on the statistic
(6)
where
(7)
is the least sum of squares of errors, *n* is the number of data points (*t*_{i}, *y*_{i}), and *ν* is the number of free parameters in the model. *y*(*t*_{i}) Represents the rate given by any of the considered models, while δ is a correction to the baseline. Namely, the baseline is initially only roughly determined from the initial recordings (see above). The minimization is performed with respect to all free parameters including δ. The model with the lowest value of AIC_{c} is considered the preferred model.

The AIC criterion does not provide a good sense to what extent one model is preferred over the other models (4, 30). Therefore, we have calculated Akaike weights for the models we compared, defined by (4, 30):
(8)
where *M* is the number of models compared, *A*_{1} is the lowest AIC_{c} for the considered models, and *A*_{k}, *k* = 2, . . ., *M* are AIC_{c} values for the rest of the models. From the weights we calculated the evidence ratio *w*_{k}/*w*_{j} for models *k* and *j*, which shows how many times model *k* is more likely the better model in terms of Kullback-Leibler discrepancy than model *j* (30).

#### Evaluation of characteristic parameters.

In previous studies it was established that emptying time *T*_{1/2} and the initial delay in gastric emptying *T*_{lag} are of interest in studying gastric emptying. The value of *T*_{1/2} is defined through the cumulative excretion curve *C*(*t*), as the time when one-half of its maximum *K* = *C*(∞) is achieved:
(9)
In the case of the PE model (*Eq. 2*) one gets the following analytic expression for *T*_{1/2}
(10)
and for the GVS model (*Eq. 3*), it is given by the solution of the nonlinear equation:
(11)
where γ(α, *x*) = ∫_{0}^{x}*e*^{−t}*t*^{α−1}*dt* is the incomplete gamma function and Γ(*x*) is the gamma function (19). The solution of this equation can be found numerically, and we used the bracketing and bisection method (20). Obviously, when *c* = 1 we get the equation for the classical GV model (*Eq. 1*). It is seen that our GVS model (*Eq. 3*) requires essentially the same effort to calculate *T*_{1/2} as GV model does.

To obtain *T*_{1/2} in two-component model GVS2, one solves *Eq. 9* with *y*(*t*) given by *Eq. 5*. This is accomplished by using the definition of incomplete gamma function for the intervals in which *y*(*t*) is not zero, whereby ∫_{u}^{v}*e*^{−t}*t*^{α−1}*dt* = γ(α, *v*) − γ(α, *u*).

Another characteristic parameter previously proposed by Symonds et al. (28) is *T*_{lag}, the time required for the stomach to grind the meal into particles fine enough to pass through the pylorus. Mathematically, it is the time when the excretion curve achieves its first maximum. For models PE and GVS, there are analytical expressions for *T*_{lag}:
(12)
In the case of GVS2 we can have two peaks and we define *T*_{lag} as time when the first peak is achieved. It should be calculated numerically. We found *T*_{lag} by numerically solving the equation *y*′(*T*_{lag}) = 0 in the neighborhood where the first peak is found just by calculating *y*(*t*) on sufficiently fine grid. The equation is solved by bracketing and bisection method (20).

Symonds et al. (28) have also considered the GEC = ln*a* based on the GV model as yet another quantitative characterization of gastric emptying rate. Clearly, it can be likewise evaluated for our GVS model. However, it is not clear how this index would be evaluated for our two-component GVS2 model. We believe that a more intuitive characterization of gastric emptying rate is its maximal value, which can be easily calculated for the models considered here. In case of PE, GV, or GVS we define the maximal rate *R* = *y*(*T*_{lag})/*K*. The function *y*(*t*) is divided by *K* = *C*(∞) to ensure that *R* does not depend on the coefficient *a* which may vary from measurement to measurement due to possibly different scaling and losses of ^{13}CO_{2}. For PE, GV, and GVS models, there are analytical expressions for *R*:
(13)
where *c* = 1 corresponds to GV model. For GVS2 we can have two peaks. Assuming they have occurred at times *T*_{1} and *T*_{2}, *R* is evaluated as greater of the two, i.e., *R* = max[*y*(*T*_{1}),*y*(*T*_{2})]/*K*.

#### Estimation of characteristic parameters by least-squares fitting.

We use the least-squares fitting to estimate the model parameters from which *T*_{1/2}, *T*_{lag}, and *R* are calculated based on the above equations. The numerical minimization in the least-squares fitting is found to proceed faster when the model parameters δ, *a*, *a*_{1}, *a*_{2} are expressed through the rest of parameters, which is possible because δ, *a*, *a*_{1}, *a*_{2} appear linearly in the model. The details are explained in the appendix. Our custom developed software uses the SIH algorithm for minimization (18), which likely finds the global minimum. The uncertainties in characteristic parameters are estimated by the bootstrap method (31), which is briefly explained in the appendix.

## RESULTS

We analyzed 175 gastric emptying data sets using the considered empirical models, namely the PE model given by *Eq. 2*, the standard GV model given by *Eq. 1*, the model with a gamma variate function that includes stretched exponential (GVS) given by *Eq. 3*, and the model with a linear combination of two gamma variate functions with stretched exponentials (GVS2) given by *Eq. 4*. The models were fitted to data according to a least-squares method as described in materials and methods and in appendix.

The most appropriate model in each case is chosen based on the AIC_{c} and corresponding Akaike weights and evidence ratios. In the majority of the gastric emptying data sets analyzed, GVS was preferred by this criterion compared with GV and PE. This is the case for example shown in Fig. 1, where GVS is preferred compared with PE. On the other hand the PE model is preferred compared with GV (Table 1, *case 1*). GVS is evidently the preferred model in terms of Akaike weights and evidence ratio (see Table 1, *case 1*). Furthermore, GVS leads to significantly different values for *T*_{1/2} and *R* than those obtained by PE or GV (Table 1, *case 1*). This example illustrates our rationale for introducing the generalization of the standard GV model. The GVS model was preferred in 90% of the analyzed cases, compared with the GV model, and in 91% of the analyzed cases compared with the PE model. Interestingly, GV was preferred in 76% cases compared with PE. This indicates that GV model is perhaps a more suitable empirical model for gastric emptying data than the PE model and a posteriori justifies our idea to generalize the GV model.

The question arises whether in the considered example (shown in Fig. 1) GVS2 model would be preferred by AIC_{c}, although data do not suggest more than one peak. When the fitting by GVS2 is performed (Fig. 2), one gets a lower value for AIC_{c}, but at the same time an obvious artifactual peak is produced, which is essentially the result of one deviant group of data. In analyzing data we adopted the criterion that when GVS2 fit leads to a sharp spike (or dip) as a result of few deviant groups of data, we consider this to be an overfitting and chose a simpler model with the second lowest AIC_{c} value as the most appropriate description of the data. While often it is visually pretty obvious from the fit when the peak (dip) is spurious, we introduced an objective criterion for spurious peaks (dips). We first overlay the best-fit curves obtained by GVS2 and GVS as in Fig. 2. Then, we estimate the width of the spurious peak (dip) from the two intersections of GVS2 curve by GVS curve at the location of the peak, which provide a measure of the width of the peak (dip) as indicated in Fig. 2 by two green lines. The intersections are numerically determined by an interactive computer program, which uses bisection method (20). The user must provide bracketing for each intersection. This is done by visually observing the intersections from the graph similar to Fig. 2 that the program displays. The times at which intersections occur (obtained by the bisection method) determine the time interval that represents the width of the peak (dip). If this interval is less than 24 min (which is generally the interval between 4 group of points), then we considered the peak (dip) to be spurious.

In the example shown in Fig. 3, GVS fit is preferred compared with GV and PE fits. It leads to different values for characteristic parameters and especially for *T*_{lag}. (Table 1, *case 2*). Data do not suggest more than one peak and indeed the fit with GVS2 leads to a sharp artifactual spike according to our criterion described above. Thus, the GVS2 fit is not considered appropriate. A similar example is shown in Fig. 4, where the differences in values of *T*_{1/2}, *T*_{lag}, and *R* for GVS and GV fits are all significant, while the GVS fit is considered preferred by AIC_{c} compared with GV and PE fits. The GVS2 fit was inappropriate with a spurious dip according to our criterion.

Such a situation is not the case for all data we analyzed. In 34 out of 175 cases we visually observed two distinct peaks (or peak and a shoulder), which implied the necessary application of the two-component model GVS2. An example is shown in Fig. 5 where the second peak is observable and by AIC_{c} GVS2 is preferable compared with GVS, GV, and PE. GVS2 leads to estimation of *T*_{1/2} which is clearly less uncertain than the estimation with GVS (see Table 1, *case 4*).

In the example shown in Fig. 6 the second peak can be barely observed, but it has considerable effect on the shape of the best-fit curve. Such shape can be only obtained by GVS2, which is also preferable by AIC_{c} (together with Akaike weights and evidence ratio) compared with GVS, GV, and PE. The evaluated *T*_{1/2} is significantly lower compared with the value obtained by fit with GVS (Table 1, *case 5*). Also, maximal rate *R* is significantly higher, while *T*_{lag} is essentially the same. This example further illustrates that in some cases GVS2 should be used to properly determine the emptying half-time and maximal rate.

We encountered cases where two peaks are distinctly observed (Fig. 7), but both models GVS2 and GVS lead to not so different value for *T*_{1/2}, while the values for lag time and maximal rate are significantly different (Table 1, case 6). In this case AIC_{c} provided preference to GVS2, as one would expect, while GVS model was preferred over GV and PE models.

In our final illustrative example we show the case where data are not suggesting two peaks, while GVS2 yields best-fit curve with one peak and unusual shape, which cannot be obtained by GVS, GV, or PE models (Fig. 8). GVS2 is the preferred model in this case, and interestingly, PE model is preferred compared with GVS and GV models. The value obtained by GVS2 for lag time is significantly different from the value obtained by PE (see Table 1, *case 7*).

Altogether the GVS2 model was preferred by AIC_{c} in all considered cases. However, according to the criterion for spurious spikes (dips), it was selected only 41.1% of the times as the most appropriate. GVS was selected 51.4% times, GV was selected only 2.9% times, and PE 4.6% times. In the cases when PE was selected as the preferred model, the best-fit values for characteristic parameters differed on average 1.5% (with maximum of 4.6%) from the corresponding parameters obtained by GVS.

## DISCUSSION

The GV model has been previously used for the analysis of ^{13}CO_{2} gastric emptying data for mice and in particular to estimate *T*_{1/2}, (11, 28, 6). These studies have to some extent validated the GV model. The objective of our study was to improve the analysis based on this model. Therefore, we have introduced a simple generalization of the GV model, which consists of replacing the exponential function with the stretched exponential function (see *Eq. 3*). Such an enhanced model has more flexibility to describe the data and consequently provides better fits that generally lead to more accurate estimation of characteristic parameters. As we demonstrated (see Table 1), it can yield values for emptying time (*T*_{1/2}) and the initial delay in gastric emptying (*T*_{lag}) that significantly differ from the values obtained by classical gamma variate model GV. The same is true when our model with stretched exponential GVS is compared with power exponential model PE (see *Eq.2*), which was previously used to analyze gastric emptying data (28, 23).

The ability of classical GV and PE models to fit observed gastric emptying data from humans has been questioned in the past. A multicomponent mathematical treatment (based on convolution integral) has been developed to try to separate postgastric processing of the label from gastric emptying (16). The presence of two peaks in the gastric emptying curves for mice has been reported but not analyzed in the past (28). Our study is the first to analyze complex gastric emptying patterns in mice. This now appears possible as data can be and are collected at high temporal resolution. Instead of including complexities introduced by the model of Maes et al. (16), we found it appropriate to introduce a two-component empirical model defined by a linear combination of two gamma variates with stretched exponentials (see *Eq. 4*). This model (GVS2) was capable of representing data with two peaks (Figs. 4–6) and more intricate one-peak data, which were better described by the a two-component model (Fig. 7) than by a one-component model (GVS). The introduction of a two-component model did have an effect on the values of the characteristic parameters (Table 1): *T*_{1/2}, *T*_{lag}, and normalized maximal rate of gastric emptying *R*. We introduced characteristic parameter *R* as perhaps a better characterization of gastric emptying rate than the GEC (see Ref. 28), which can be reasonably defined only for one-component models.

We have developed our own software for ^{13}CO_{2} data analysis, partially in R (data preprocessing) and in FORTRAN (fitting). In it we implemented some features that appear required when working with ^{13}CO_{2} measurements from mice, especially younger mice. These are the ability to remove dropped points that are clear artifacts and accurate correction of the background. Both issues arise due to the low signal obtained from an animal weighing less than 20 g, which generates small amounts of CO_{2}. Occasionally points need to be dropped from the curve due to blocks in the tubing by condensation or instability in the calibration of the machine. Both of these blockages result in a failure to discriminate the spectroscopic signals from H_{2}O and ^{13}CO_{2} and generate a clear artifact in the curve at a single time point. The second issue is the difficulty of accurately setting the baseline when the signal is very small in relation to the noise in the system. This issue has been addressed by introducing an additive constant to the model to ensure that the data are fitted without distortion by an inaccurately set baseline.

A correction for failure to recover all the ^{13}C fed to the animal is implicitly incorporated in coefficient *a* (*Eqs. 1*, *2*, and *3*) or in *a*_{1} and *a*_{2} (*Eq. 4*). This is important since a small but significant amount of ^{13}C is incorporated into the bone and/or excreted in the urine and feces (16). Incomplete fasting of the mice can also be detected and can be a confounding factor in this type of study due to the higher than background levels of ^{13}C in standard mouse chow (28).

To investigate the merit of the considered four models (PE, GV, GVS, and GVS2) we applied them to 175 data sets and used AIC to select the preferable model for each data set. In most cases the one-component GV model with stretched exponential GVS was selected as the optimal representation of gastric emptying data with respect to goodness of fit and the number of free parameters. Of course this was not the case when two peaks were observable or when the shape of the gastric emptying curve was considerably different from the shapes that the GV model with stretched exponential can provide (Fig. 7).

The two classical models GV and PE were rarely selected as a preferred model, i.e., 7.5% times. One may pose the question whether in analyzing gastric emptying data one has to try these models as well. Our analysis indicates that this is not necessary, because the GVS provides essentially the same values for the characteristic parameters. This is, of course expected when the GV is selected, because this model is nested in GVS, i.e., GV is obtained from GVS by fixing parameter *c* to 1. In the case where PE was selected the difference was also insignificant.

Using a model that fits data well has a significant impact on the interpretation of the gastric emptying data. For example, in *case 3* of Table 1 (corresponding to Fig. 4) *T*_{1/2} provided by GVS model is 153 min with large uncertainty of 43 min, while GVS2 value is 109 min with uncertainty of 14 min. In our studies 153 min is considered delayed, whereas 109 min is considered normal. Since delayed gastric emptying is associated with defined anatomical and biochemical changes (6, 7), then the accuracy of our determination as to whether the mice have delayed or normal gastric emptying affects the interpretation of these associations.

In summary, our study points out that two factors necessary to fit the data more effectively are the stretching of the GV model and the use of a two-component model to fit data with two clear peaks. Stretching the GV model may better reflect the influence of metabolic differences at later time points when gastric emptying contributes less to the overall levels of ^{13}CO_{2} released in the breath as reported in analysis of human data (16). Future studies can be directed towards testing the proposed models in studies on gastric emptying in humans where they might be useful in detecting changes in gastric emptying that do not show up with the existing approaches to analyzing breath test data.

## GRANTS

This work was supported by National Institute of Diabetes and Digestive and Kidney Diseases Grants DK-68055 and DK-57061.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## AUTHOR CONTRIBUTIONS

Author contributions: Ž.B., S.J.G., D.R.L., and G.F. conception and design of research; Ž.B., S.J.G., H.D.C., D.R.L., and G.F. analyzed data; Ž.B., S.J.G., H.D.C., D.R.L., and G.F. interpreted results of experiments; Ž.B. and S.J.G. prepared figures; Ž.B., S.J.G., H.D.C., D.R.L., and G.F. drafted manuscript; Ž.B., S.J.G., H.D.C., D.R.L., and G.F. edited and revised manuscript; Ž.B., S.J.G., H.D.C., D.R.L., and G.F. approved final version of manuscript.

## ACKNOWLEDGMENTS

We thank Kristy Zodrow for help with the preparation and submission of the manuscript and Jessica Phillips and Stephanie Hein for work carrying out the gastric emptying studies utilized in this manuscript. We also thank Chetan Offord for assistance in preparing the figures. Finally, we thank the reviewers for comments that helped to considerably improve the manuscript.

## Appendix

Here we first describe in some detail how the least-squares fitting was performed. Then, we explain how the uncertainties in characteristic parameters were estimated.

The models involved contain free parameters that are arranged linearly in the model function. For example in case of GVS model the sum of squares of errors (SSE) is given by:
(14)
where δ is the baseline correction. The minimization should be performed with respect to *a*, *b*, *c*, and δ. It was found useful to solve for linear parameters *a* and δ, which at the minimum should satisfy:
(15)

yielding
(16)
where
(17)
Thus the least-squares method is reduced to minimizing SSE given by *Eq. 14* with respect *b*, *c*, τ while δ and *a* are given by *Eq. 16*. In calculating *f*(*t*) as given in *Eq. 14*, we found that numerical overflow is better controlled if *f*(*t*) is expressed as exp[−(*t*/τ)^{c} + *b* ln *t*]. In some cases, the fitting leads to baseline correction δ that is unacceptably large in absolute value compared with the mean *ȳ* of all data (*Eq. 16*). In this case we introduce the penalty function which requires that |δ/ȳ| < 0.5.

We used the same procedure of solving for *a* and δ for GV and PE models. The fitting with GVS2 model is more challenging due to higher number of free parameters. However, it is still possible to reduce the number of fitting parameters by analytic method outlined above for the GVS model. In this case we search the minimum of
(18)
with respect to *b*_{j}, τ_{j}, *c*_{j}, *j* = 1, 2. Zero partial derivatives with respect to *a*_{j}, *j* =, 1, 2 and δ provide the system of three linear equations that can be solved yielding *a*_{1}, *a*_{2} and δ as functions of *b*_{1}, *b*_{2}, τ_{1}, τ_{2}, *c*_{1}, *c*_{2}. Then, the numerical minimization with respect to these parameters is performed via SIH algorithm. In relatively rare cases when the fit function becomes negative for some *t* it is set to zero according to *Eq. 5*. This in effect restricts allowable parameters *b*_{1}, *b*_{2}, τ_{1}, τ_{2}, *c*_{1}, *c*_{2}. We found this acceptable, as the number of free parameters involved appears more than sufficient to describe various patterns of excretion curves.

Once the best-fit model parameters are found we can calculate gastric emptying characteristic parameters *T*_{1/2}, *T*_{lag}, and *R*, according to formulas discussed in materials and methods. However, it is also of interest to estimate the uncertainties in these parameters. We used the bootstrap method (31) assuming that measurement errors in data are distributed with an identical but otherwise not necessarily known distribution. This assumption is consistent with the nature of the measurements involved. In brief, the residuals *r*_{i} = *y*_{i} − *y*(*t*_{i}), *i* = 1, . . ., *n* obtained by fitting the original data are considered to be a sample of errors, which are then used to construct synthetic data according to the prescription:
(19)
where *l* is an integer picked at random from the set {1, 2, . . .,}. Synthetic data are then fit by the model and characteristic parameters obtained as in the case of actually measured data. We used 200 synthetic data sets and obtained a population of T_{1/2}, *T*_{lag} and *R*. The uncertainty in parameters is then estimated by the standard deviation for that population.

- Copyright © 2015 the American Physiological Society