## Abstract

The detection of hormone secretion episodes is important for understanding normal and abnormal endocrine functioning, but pulse identification from hormones measured with short interval sampling is challenging. Furthermore, to obtain useable results, the model underlying hormone secretion and clearance must be augmented with restrictions based on biologically acceptable assumptions. Here, using the assumption that there are only a few time points at which a hormone is secreted, we used a modern penalized nonlinear least-squares setup to select the number of secretion events. We did not assume a particular shape or frequency distribution for the secretion pulses. Our pulse identfication method, VisPulse, worked well with luteinizing hormone (LH), cortisol, growth hormone, or testosterone. In particular, applying our modeling strategy to previous LH data revealed a good correlation between the modeled and measured LH hormone concentrations, the estimated secretion pattern was sparse, and the small and structureless residuals indicated a proper model with a good fit. We benchmarked our method to AutoDecon, a commonly used hormone secretion model, and performed releasing hormone infusion experiments. The results of these experiments confirmed that our method is accurate and outperforms AutoDecon, especially for detecting silent periods and small secretion events, suggesting a high-secretion event resolution. Method validation using (releasing hormone) infusion data revealed sensitivities and selectivities of 0.88 and 0.95 and of 0.69 and 0.91 for VisPulse and AutoDecon, respectively.

- secretion
- pulsatility
- AutoDecon
- model complexity
- model selection
- nonlinear fitting
- infusion study
- simulation study

use of radioimmunoassays for hormone measurements in the 1960s revealed that many hormones are not released continuously into the blood stream but rather in specific hormone-dependent patterns characterized by diurnal rhythms and episodic secretion. Further studies have revealed that the secretion profiles of different hormones, such as those of the somatotrope axis and gonadotropic system, are dependent on age and/or sex. In addition, diseases such as acromegaly and infertility are associated with remarkable changes in secretion patterns (26, 30). Such results indicate that secretion patterns carry useful information on normal and abnormal functioning. The first reports on hormone secretion patterns were purely descriptive, but several methods have since been developed for quantifying pulses and other secretion parameters.

One class of such methods does not explicitly use models but characterizes time series by summary statistics. For instance, secretion regularity has been investigated with the approximate entropy (ApEn) statistic (23, 24), which allows quantification of the degree of regularity or reproducibility of patterns in a time series. The ApEn is a model-free, translation- and scale-invariant, asymptotically Gaussian-distributed statistic that distinguishes the orderliness of data series of 30 or more points in length. Furthermore, ApEn evaluates both dominant and subordinate patterns in data; notably, it can detect changes in the underlying episodic behavior not reflected in the peak events or amplitudes. The Pulsar and Cluster programs (19), as well as those described by Santen and Bardin (25), are additional model-free methods.

Summary statistics are useful for describing the particular properties of a process; however, such processes are often too complex to be described by a single statistic. A second class of methods considers explicitly the underlying models for hormone concentration time series (7, 8, 19). These methods include Detect and the multiparameter methods developed by Johnson et al. (5) and Veldhuis et al. (28). A model for secretion and irreversible metabolic clearance is also well accepted. The burst-like secretion episodes are often modeled with Gaussian, Weibull, or other distribution shapes, and irreversible concentration-dependent loss is facilitated by several routes, including liver metabolism and kidney excretion. Dilution into the blood compartment, receptor-mediated signaling in target organs, and intracellular metabolism are not modeled, because much remains unknown about their quantitative aspects.

Various decay-secretion models exist that differ largely in how the secretion or decay processes are formulated. For example, secretion may include basal secretion, and the exponential decay phase may consist of multiple decay rates. It is unclear which combination of secretion and decay functions most accurately reflects the underlying biological process. A problem with models based on combinations of these functions is their mathematical unidentifiability, such that unique solutions cannot emerge from mathematics alone.

Here, we present a model framework for endocrine pulse identification that centers on assuming episodic hormone secretion and model simplicity. Model simplicity was tuned by making a tradeoff between the number of parameters and the fit. This framework allowed for pulsatility analysis without assumptions about the pulse shape or frequency distribution being made, thus generating a versatile method for endocrine pulse analysis.

## MATERIALS AND METHODS

### The Model

The experimental data included time series of 145 samples of serum hormone concentrations obtained at equidistant 10-min intervals over 24 h. These time series were obtained for different individuals and hormones. The data came from an experiment that compared the hormone dynamics of patients with narcolepsy and matched healthy controls (12, 13, 14, 15, 21) from a dopamine D2 agonist intervention in obese women (10, 11) and from an infusion study (9).

Current models for serum concentration time series are based on episodic secretion and exponential hormone decay (8, 26, 28). The decay, or metabolic clearance, is the fraction by which the hormone concentration is reduced over a certain time period. A pseudo-model that captures decay and secretion is given by

This model assumes a constant decay during the period considered and equidistant sampling. In words, *Eq. 1* indicates that the concentration *x*_{j} equals the preceding concentration multiplied by the decay (*x*_{j−1}*Decay*) plus the secretion and residual error (including both measurement and model errors) at *time point j* (see appendix for a more detailed explanation of the notation). If the decay function is assumed to be monoexponential (*e*^{−kj}) and the parameter *k* is expressed in units of 10 min, then the model (*Eq. 1*) can be specified further as
*u*_{j} represents the residual error and *φ*_{j} is the (presumed nonnegative) secretion at *time index j*. The implicitly defined (*x̂*) is the modeled part of the data.

*Equation 2* is an example of *Eq. 1*, where the *Decay* term is given by *e*^{−k} and *Secretion* by *φ*_{j}. This is the simplest description of the base model (*Eq. 1*), but more complex instances are conceivable. For instance, sometimes there are reasons to believe that hormone decay is not monoexponential but bi- or triexponential (20). The base model can easily be adapted to describe this behavior. However, often the fast decay rate is hard to estimate from short time series, because measurement errors make discrimination from noise difficult. In a later section, we explore whether adding exponents and/or baseline secretion to the model is useful.

### The Problem

Using the model described above, luteinizing hormone (LH), which is known for its episodic secretion (26), was used to investigate the pulse identification problem in greater detail. *Equation 2* was fitted to the data such that the sum of the squared residuals *φ*_{j} was nonnegative.

Figure 1 shows the result for a 24-h LH profile, with estimated secretion profiles and decay. This figure reveals two fundamental problems of the modeling procedure. The first is overfit; the modeled concentration matched the measured hormone concentrations exactly. This implies that the residuals are zero, which is unrealistic considering the measurement noise. The second problem is nonuniqueness; many solutions could be found that gave the same (zero) residual error. These problems are explained in more detail in the appendix. From a biological perspective, the secretion is expected to be episodic rather than continuous. Therefore, solutions of this unconstrained model (*Eq. 2*) have no biological meaning. Considering the large body of evidence in favor of this model (3, 5, 8, 20, 22, 26, 29, 31), the model specification is unlikely to be the problem. Therefore, to obtain sensible results, biological knowledge must be used to apply additional constraints in the fitting process of *Eq. 2*. The constraints on the base model differ between methods, and our model constraint limits the number of secreting pulse times.

### The Solution

Although there is consensus in the literature about the basic model, a divergence of opinion exists in the ways suggested to impose these extra constraints. Because hormone secretion is purportedly episodic, at many time points no secretion should be found. Several researchers have investigated various ways to arrive at such a solution, thereby devising different approaches (3, 20, 29, 31). A few studies (3, 20) used *Eq. 2* without simultaneously solving for the decay rate and secretion episodes, instead identifying the secretion episodes for a given decay rate. Yang et al. (31) proposed a spline-based fitting, variable (diurnal) baseline secretion and parameterized secretion form. However, it is unclear how to choose between variable baseline secretion and variable decay rates, because they are interdependent. Johnson and colleagues (5, 6) proposed a model similar to our own, with (optional) basal secretion and Gaussian-shaped secretion events. This latter model is currently the most commonly applied in the field.

In this study, only the number of secretory episodes was constrained. We searched for an optimal balance between the fit of the data and a minimum of secretory events. To strike this balance, the signal-to-noise-corrected Akaike's information criterion (AICu; *Eq. 16* in the appendix) was used. The AICu makes a tradeoff between model complexity (the number of parameters used) and the model fit. In our case, the model complexity was directly related to the number of secretory events, *φ*_{j}, because allowing for more episodes made the model more complex as it carried more parameters to estimate.

The attractiveness of this approach is that it is free from assumptions about the shape of secretion. This aspect makes it a general pulse identification model, which suits the exploration of secretion patterns. The model explicitly sets the number of silent secretory times, thereby avoiding too many false positives. More details on the fitting and selection procedures are given in the appendix. The newly developed fitting procedure was coined VisPulse.

The AICu can also be used to choose between common model specifications, e.g., *1*) monoexponential decay with episodic nonbasal secretion (our reference model), *2*) monoexponential decay with episodic and basal secretion, or *3*) biexponential decay with episodic secretion. The use of a generic information criterion thus enables us to select a model specification.

### Method Validation

#### Method validation using synthetic data.

Method validation is not simple, since the timing of the secretory episodes and the decay rate constant are generally unknown. In such cases, computer simulations on synthetic data can be used and the outcomes compared with the true values. Unfortunately, a first-principles model is not available to generate such synthetic data, because an important goal of endocrine pulse analysis is to find such models. Setting up a representative synthetic secretion profile is complicated because it is not well understood which of the estimated pulses are true pulses. Some pulses are without a clear biological effect and occur at unknown frequencies and amplitudes, complicating their biological validation.

Expert knowledge may be used to construct synthetic secretion profiles, but an alternative is to use the secretion pulse profile identified by our method. Simulations based on derived synthetic data, called realistic simulations, are broadly accepted as a means for validation and method comparison (2). However, such simulation studies rely strongly upon the assumptions of the model that estimated the realistic profile. Although the basic model underlying VisPulse and AutoDecon is identical, the additional assumptions and constraints on the secretion terms are not, and these are dominantly present in the estimated secretion profile. Thus, even in the context of a synthetic truth, the assumptions used in deriving this truth influence the model comparison.

Therefore, to eliminate bias, two simulation studies were performed, one that used a synthetic truth estimated by VisPulse and one that was estimated by AutoDecon. The data used in the simulation study were based on the hormone concentrations estimated by both methods. Thus, the secretion pulses and clearance rate were known, and a random model residual was added to induce variation in the concentration profile. The synthetic hormone concentration profiles were based on a measured LH concentration profile and defined as
*x̂*_{syn,j} equals the *x̂*_{j} of *Eq. 2* and *ε*_{j} is a random number obtained from randomized fit residuals. A total of 25 (VisPulse) and 20 (AutoDecon) profiles were made, with independent noise realizations derived from the model fit residual.

The constructed synthetic hormone concentration profiles were subsequently analyzed by VisPulse and AutoDecon (5). The secreted mass, half-life, sensitivity (rate of correctly finding a pulse), selectivity (rate of correctly finding no pulse), and zero secretion rate (*q* value) were calculated, as were the sums of the squared residuals between the modeled and real concentration profiles (model residual) and between the modeled and real pulse profiles (pulse residual).

#### Method validation using releasing hormone infusion.

A limitation of simulation studies is their inherent dependence on the underlying model and intrinsic discrepancy between real and simulated data series. In a physiological context, releasing hormone (RH) infusion data sets are commonly used to assess the secretion response in vivo. Subjects that have an induced RH deficiency are used, and the controlled infusion of RH causes the relevant tissue to secrete the hormone with a reasonably well-known temporal relation. Smaller secretion events may emerge between the infusions, potentially caused by heterogeneity or refraction differences in the secreting tissue. However, the secretion events directly following the infusion are commonly severalfold larger than those between infusions. The emergence of small secretion events also depends on the specific hormone and subject. For example, LH seems to suffer more from this effect than growth hormone (GH) or testosterone (see supplemental documentation; Supplemental Material for this article is available at the *AJP-Endocrinology and Metabolism* web site).

Three RH infusion experiments were analyzed by VisPulse and AutoDecon. The experiments consisted of a measured LH, GH, or testosterone response to the RH infusion (9). The LH experiment consisted of six consecutive pulses (10-min infusion) of LH followed by a single gonadotropin-releasing hormine (GnRH) pulse, allowing a method comparison on fully known pulses. The GH and testosterone method comparison was only qualitative, because the methods had different assumptions about the underlying processes.

## RESULTS

The modeling strategy applied to the earlier-mentioned LH data is illustrated in Fig. 2 and revealed that the modeled LH concentrations matched the measured hormone concentrations well. The estimated secretion pattern was sparse, consisting of long silent periods (especially during the night) and few secretion time points. The larger peaks were evenly distributed over the day, whereas the smaller events tended to cluster together. The residuals were small compared with the measured hormone profile and were without obvious structure. The estimated decay rate, or metabolic clearance rate, of LH in this subject was 108 min, slightly longer than the previously reported 75 min (13).

More complex models were also fitted to the hormone concentration data. In the biexponential model, a faster decay rate was prescribed by our understanding of dilution in the transport (blood, serum, and interstitial) compartments. In most instances, this rate exceeded the sampling frequency of six measurements per hour, which may explain why the biexponential model specification was not selected in most of our data by the AICu. An exception to this was the infusion experiment of testosterone. Model fitting is detailed further in the appendix.

### Comparing VisPulse and AutoDecon

We used VisPulse and AutoDecon, a commonly used method for hormone secretion estimation, to model a cortisol concentration series (10, 11), and representative results are shown in Fig. 3. AutoDecon underestimated the high cortisol peaks in the middle of the time profile, whereas a large overestimation occurred at the end of the time series. AutoDecon usually identified fewer secretion events; therefore, the VisPulse fit was somewhat better because it led to lower residual errors.

To compare the results between VisPulse and AutoDecon in a more quantitative manner, various estimated parameters were compared. The estimated clearance rate was 76 min for VisPulse and 81 min for AutoDecon; both were slower than the previously reported 66 ± 18 min (16, 17). The percentage of zero secretion points (*q*) was 75% for VisPulse and 60% for AutoDecon (36 and 58 secreting times, respectively). The secretion mass was estimated at 2,222 U·l^{−1}·day^{−1} for VisPulse and 2,113 U·l^{−1}·day^{−1} for AutoDecon. Finally, the sum of squared residuals was slightly higher for AutoDecon. Table 1 shows the average results in detail and gives standard deviations calculated over the 36 individuals in our study. The most important difference was the larger number of silent secretion periods of VisPulse.

### Method Validation

#### Method validation using synthetic data.

As mentioned above in materials and methods, the synthetic data depend upon the assumptions used to create the data, and the assumptions and constraints on the secretion terms were dissimilar between VisPulse and AutoDecon. This is visualized in Fig. 4, where a small part of the simulation study is shown to illustrate the problem. To eliminate bias, we performed two simulation studies using a synthetic truth estimated by VisPulse or by AutoDecon. Tables 2 and 3 shows the results of the VisPulse-derived (Table 2) and AutoDecon-derived (Table 2) synthetic profiles.

Although the simulated true values of the secreted mass and half-life of the simulated approaches were rather different (VisPulse had a much lower true mass and higher half-life than AutoDecon), both methods estimated these values well. The major differences between the methods were found in the selectivity and sensitivity results. The sensitivities of the methods were similar in VisPulse-based simulations, whereas the selectivities were similar in AutoDecon-based simulations. However, the selectivity of VisPulse was much better in VisPulse-based simulations, whereas the sensitivity of AutoDecon was much better in AutoDecon-based simulations. The lower sensitivity of VisPulse in the AutoDecon-based simulation can be explained by the assumed Gaussian pulse shape, which is recovered better in a method that constraints the pulses to be Gaussian (see also Fig. 4). The *q* value was higher in VisPulse than in AutoDecon, indicating that VisPulse estimated less secreting occasions than AutoDecon. These major differences between the methods are shown in detail in Fig. 4.

The VisPulse secretion profile contained more unit width and small amplitude secretion events than the AutoDecon profile. The combination of the model and pulse residual indicated that the large amplitude pulses were correctly identified, although AutoDecon estimated these events to be wider than the synthetic profile. In the AutoDecon scenario a different situation was observed, since VisPulse estimated the events too narrowly given the synthetic profile, although the total mass estimate was similar in accuracy to AutoDecon. The true negative ratio was identical in the AutoDecon scenario, but VisPulse clearly outperformed AutoDecon in the VisPulse scenario due to pulses being too wide for the synthetic profile.

#### Method validation using infusion experiments.

Figure 5 displays the results of an experiment measuring LH. This experiment consisted of six LH infusions and one GnRH infusion. The LH was infused over 10 min; since this interval was equal to the sampling interval, it was safe to assume a block pulse of 10 min to be wide.

Figure 5*A* displays the LH concentration profiles of the first five subjects. Most profiles clearly showed the infused LH pulse with the correct width at the expected position. Figure 5*B* shows for each time index whether there was an estimated secretion. For clarity, the color scheme was inverted for VisPulse vs. AutoDecon. The results suggested that VisPulse had a higher temporal resolution than AutoDecon, since the pulses appeared to be more strictly defined to a specific time index. Figure 5*C* is the average secretion over all subjects per method, which indicated a well-performed experiment and differences in secretion event length between VisPulse and AutoDecon. The sensitivity and selectivity for VisPulse were 0.88 and 0.95, respectively, and for AutoDecon were 0.69 and 0.91, respectively.

Using either method, we also noted very small secretion events between the infusions, which we speculate resulted from the refractionary principles detailed above. Finally, in these infusion experiments the same hormone peak identification differences that were discussed above and shown in Fig. 4 were found again. AutoDecon typically identified the peak at an earlier time point than VisPulse and estimated the peaks to be broader.

Figure 6 displays the recovered GH pulse episodes induced by GHRH infusions. The VisPulse-recovered pulse episodes are generally narrower and variable in width when compared with AutoDecon; the latter also misses more induced pulse episodes.

Figure 7 displays the recovered testosterone pulse episodes induced by LH infusions. Again, the VisPulse-recovered pulse episodes are narrower while being confined to strictly defined time slots over subjects, whereas the AutoDecon-recovered pulse episodes are much wider.

## DISCUSSION

Applying our modeling strategy to previous LH data revealed a good correlation between the modeled and measured LH hormone concentrations, and the estimated secretion pattern was sparse, in agreement with biological expectations. The small and structureless residuals indicated a proper model with a good fit. VisPulse appeared to generate a larger number of silent secretion periods, which was verified by method validation with synthetic data. Method validation using (RH) infusion data revealed that the sensitivity and selectivity for VisPulse were 0.88 and 0.95, respectively, and for AutoDecon were 0.69 and 0.91, respectively.

Pulsatility analysis is complicated, since the true mechanistic details of secretion and clearance are essentially unknown. There is a large body of evidence suggesting sparse or episodic secretion and exponential clearance. However, models that combine secretion and clearance do not result in unique secretion profiles due to mathematical problems. Different solution paths have been tried. For instance, AutoDecon assumes parameterized (Gaussian) secretion pulses. This assumption works well in fitting the model parameters but may not closely match the underlying secretion mechanism. In contrast, our VisPulse approach did not assume any shape for the hormone pulse. This idea of sparsity or simplicity is well rooted in science and goes back to William of Ockham (Ockham's razor). Implementing this idea using a penalized method gave very good results and thus seems legitimate.

In conclusion, the sensitivity and selectivity statistics suggest that both methods did reasonably well in recovering the true secretion masses and metabolic clearance rates and in identifying pulses. The largest difference between the two methods was the significantly higher VisPulse selectivity, suggesting that this method correctly identified more time periods of secretory silence.

## DISCLOSURES

No conflicts of interest are declared by the author(s).

## ACKNOWLEDGMENTS

We graciously thank Johannes D. Veldhuis, Mayo Clinic (Rochester, MN), for providing us with the infusion experiment data, which were very helpful in strengthening our argument.

- Copyright © 2010 the American Physiological Society

## APPENDIX

This section explains the mathematics used in this paper in more detail. A basic understanding of calculus, linear algebra, and the use of nonlinear fitting methods is assumed throughout this appendix. We performed the entire procedure in a MathWorks Matlab version 2007-A. The notations throughout this appendix follow the following markup definitions: the vectors are in lower case boldface italics (**x**), and scalars are in lowercase lightface italics (*x*). The experimental time series used here were comprised of 145 equidistant samples over 24 h. The first measurement point (*j* = 1) was taken at 9:00 AM, and the samples were 10 min apart so that the second sample (*j* = 2) was taken at 9:10 AM. The last sample (*j* = 145) was taken at 9:00 AM on the next day.

#### Model Fitting

##### Monoexponential model fitting.

In the model given by *Eq. 2*, we used *x*_{j} below for the concentration level at time index *j* and *φ*_{j} for the hormone secretion at *j*. We want *u*_{j}, which describes the fit error at *j*, to be small and nonsystematic. The decay was *e*^{−kΔt}, where *Δt* was the difference between successive time points. Since the time difference was 10 min, which matched our sampling frequency, we omitted *Δt* from the equation. Thus, *e*^{−k} depicted the decay ratio for 10 min (1 time unit). Finally, the predicted concentration at time point *j* was introduced as a new variable, *◯*_{j}:

The residual error was defined as the difference between the predicted and measured concentrations
*k* and φ that give the minimal fit error *u*_{j}:

*Equation 6* has several solutions, as can be seen in Fig. 1. In mathematics, this phenomenon is called nonidentifiability. Basically, it is the same problem as solving the equation *y* + *z* = 1 without knowing either *y* or *z*; the equation has an infinite number of solutions and is therefore not identifiable. For our model, it can also be shown that, from a single solution, an infinite number of equivalent solutions can be constructed, given *k*_{α} > *k*_{β} (*Eqs. 7*–*9*):

Also, consider a model with some different decay and secretion parameters (*Eq. 8*):

This solution with *k*_{β} and *φ*_{β,j} is equivalent to the solution with *k*_{α} and *φ*_{α,j} if

Thus, it is impossible to identify only one solution to *Eq. 4*.

A related problem is that the fit of the model is perfect (i.e., the error is zero at all time points). This is a model artifact, since measurement noise is always present, and is called overfit. One way to tackle the problems of identifiability and overfit is to restrict the solution space of the model parameters. The restrictions should come from the biological domain. The restriction that was used here was the number of possible secretion events. Thus, φ should hold many zero-value elements. However, the true number of secretory events is unknown and varies for each hormone and subject. This problem is dealt with using a search strategy for the number of events. An extra parameter is introduced for the number of zero φ elements set by the ratio *q*. Its operator (named *Q*) selects the lowest percentage of the φ elements. For example, if the lowest 40% were chosen, then *q* = 40%. The ratio *q* can vary from 0 to 100%.

To enforce the solution of the search to have a large number of zeros is to put a penalty on the nonzeros in the fraction *q*. To ensure that the penalty induced by the restriction weighs sufficiently in the total minimization criterion, it is multiplied by five times the maximum value of the concentration.

The model is solved by minimizing the sum of the residual and the penalized elements:

The setting for the scalar penalty weight depends on the relative dominance of the model fit vs. the penalty weight. The penalty should dominate the solution in the φ space; thus, multiplication by the data maximum generally suffices. However, if the variation in the data is limited, the maximum data value might not be dominant enough; for that reason, this value was multiplied by five. Proof that this multiplication does not greatly influence the results is given in the supplemental documentation.

##### Biexponential model fitting.

Biexponential model fitting is similar to the monoexponential fitting discussed above. The difference is that a decay function was added, and an extra scalar (α) dictated how much of the total decay went via the fast (1 − *α*) or slow (*α*) decay route. Assuming a nonzero *φ*_{i}, the biexponential decay of this secretion event from *time index i* to *time index j* is defined by

The decay parameters should be chosen such that the *k*1 decay is faster than the *k*2 decay. The concentration at *j* is a summation over all secretion events up to *time index j*, considering the decay of each secretion event. The initial offset or start concentration is modeled using a monoexponential slow decay only, because it is unknown how to model this in a biexponential model. The following equations combine the preceding concepts:

When *α* = 1, these equations reduce to the monoexponential model. The result of the biexponential fitting is a parametrization of the equation

##### Monoexponential model with basal secretion.

The monoexponential model with basal secretion follows the same strategy as the monoexponential fitting. The difference is that a basal secretion is added, effectively increasing the secretion at the φ element with basal secretion. Some of the estimated secretion bursts may be of slightly different masses. Therefore, the solver is given the secretion positions and estimates the secretion masses with the basal secretion:

Fitting is helped by seeding with the secreting positions of the best monoexponential φ; then the *basal* secretion and decay rate are estimated.

#### Entropy-Based Model Selection

The episodic secretion of most hormones (26) translates to a sparse φ. Thus, we preferred a model with the least parameters and an adequate model fit error. AIC is a measurement that allows for this model selection (1, 4). The AIC may be minimized through the number of parameters to form a tradeoff between the model fit (i.e., sum of the squared residuals) and the complexity (i.e., number of parameters). In this formula, *n* is the sample count and *K* is the number of parameters, which is the number of nonzero φ elements plus 1 for the exponential decay. In the biexponential case, two extra parameters are added to *K*: one for the extra decay term and one for the ratio. For problems with few observed samples relative to the number of parameters, the use of AICu is recommended by McQuarrie and Tsai (18):

#### Sensitivity, Selectivity, Secretion Mass, and Clearance Rate

The sensitivity and selectivity statistics indicate the ratio of the correct classification, and are defined as
*T* and *F* stand for true and false, respectively, and *P* and *N* for positive and negative, respectively. When both of these statistics are 1, the classification is perfect, such as when every healthy patient is classified as healthy and every diseased patient as diseased.

The total secretion mass is an indicator of the total hormone secretion in 1 day and is simply the sum of the φ:

The metabolic clearance rate is expressed in minutes using the following formula in which *si* is the sampling interval, *si* = 10 (min), in our data:

#### Method Comparison

In this article, we compared our method to a different method, AutoDecon (version 2008.09.28). The basic model for both methods is identical, but the definition of the pulses is handled in different ways. In VisPulse, the pulses are without shape, whereas in AutoDecon, Gaussians form the pulses. AutoDecon provides output that is easy to compare with VisPulse. We included only pulses >0.01 [arbitrary units (AU)/l]. This restriction is a technicality, because the numerical techniques used will not always give zero amplitude pulses but minutely small (e.g., ≤10^{−6}) pulses that are not realistic and are thus purged.

#### Nonlinear Model Solver Details

The model was implemented in MathWorks MATLAB version 2007-A on openSuSE 11.0 64bit, using a nonlinear least-squares minimizing solver from the specialized Optimization toolbox. MathWorks documents that this solver uses a subspace trust region method and is based on the interior-reflective Newton method. We set the iteration ceiling to 15,000 and the iteration tolerance change to 10^{−3}, which were sufficient to solve the problem. Since nonlinear solvers are sensitive to starting conditions, we used a rational start based on the five highest first derivative values of the concentration profile. The software will be made available to the general public upon publication as a standalone program on our website (www.bdagroup.nl).

#### The Uniqueness and Optimality of the Solutions

The use of nonlinear solvers for fitting nonlinear equations frequently suffers from what is known as a local minimum. Such a solution is not optimal, but the search strategy of the solver is unable to find a better solution. Frequently, this is related to the initialization or starting values of the fitting procedure. However, using different starting values increases the likelihood of finding a global minimum, thereby corroborating that the solution is a global minimum, without a formal mathematical proof. For similar reasons, we cannot prove that the solution found is unique.