A mathematical model of luteinizing hormone release from ovine pituitary cells in perifusion

K. Heinze1,2, R. W. Keener3, and A. R. Midgley Jr.2

Departments of 1 Biomedical Engineering and 3 Statistics, Reproductive Sciences Program, and 2 National Center for Infertility Research at Michigan, University of Michigan, Ann Arbor, Michigan 48109-0404

    ABSTRACT
Top
Abstract
Introduction
Methods
Results
Discussion
References

We model the effect of gonadotropin-releasing hormone (GnRH) on the production of luteinizing hormone (LH) by the ovine pituitary. GnRH, released by the hypothalamus, stimulates the secretion of LH from the pituitary. If stimulus pulses are regular, LH response will follow a similar pattern. However, during application of GnRH at high frequencies or concentrations or with continuous application, the pituitary delivers a decreased release of LH (termed desensitization). The proposed mathematical model consists of a system of nonlinear differential equations and incorporates two possible mechanisms to account for this observed behavior: desensitized receptor and limited, available LH. Desensitization was provoked experimentally in vitro by using ovine pituitary cells in a perifusion system. The model was fit to resulting experimental data by using maximum-likelihood estimation. Consideration of smaller models revealed that the desensitized receptor is significant. Limited, available LH was significant in three of four chambers. Throughout, the proposed model was in excellent agreement with experimental data.

desensitization; gonadotropin-releasing hormone; secretion

    INTRODUCTION
Top
Abstract
Introduction
Methods
Results
Discussion
References

ONE OF THE HALLMARKS of reproductive hormones is their pulsatility. Quasi-regular pulses of gonadotropin-releasing hormone (GnRH) are released from the hypothalamus at variable frequency, often approximately once an hour, and travel through the portal circulation to the pituitary (20, 21). There they stimulate the coincident release of luteinizing hormone (LH). Normally, LH release is also quasi-regular and pulsatile, with an amplitude that varies with the input. However, high-frequency or continuous administration of GnRH disrupts the LH-GnRH relationship, and LH production is reduced, a process termed desensitization (18). Reduced LH retards follicular development, steroid production, and ovulation. GnRH and LH also display different pulse shapes. Typically, in vivo, GnRH input pulses closely resemble square waves; the LH secretion profile exhibits an initial peak followed by an exponential decay (20).

Given the complexity of the GnRH-LH relationship, it is useful to augment experimental study with a theoretical approach, formulating and fitting a model. The model developed here organizes the relationships among variables, describing system behavior concisely and quantitatively. The model provides estimates for unobserved system variables over time and can be used to predict future system response. By formulating the model in a quantitative framework, statistical procedures can be used to assess agreement between model and experimental data, testing specific physiological hypotheses. Here these tests provide evidence of receptor desensitization.

Our goal is to model the observed behavior of ovine pituitary cells in vitro under the pulsatile administration of GnRH, including the nonlinear phenomenon of desensitization. The true mechanisms of LH desensitization are not clear. Our intent was not to discover the details and causes of gonadotrope desensitization but, rather, to determine whether proposed mechanisms can account for observed behavior. The model is based on underlying physiological concepts and experimental observations. Observed relations between diminishing LH response pulse area and increasing input frequency suggest that the pituitary requires recovery time between application of GnRH pulses to produce regular LH pulses. Lessened LH response may result from insufficient time between pulses to reactivate receptors and/or replenish LH pools. Our in vitro model postulates two schemes to account for desensitization: a desensitized receptor and limited LH availability. Because cells cannot replenish LH immediately, the second mechanism must be present, at least to some extent. Receptor desensitization is of more interest, because it may account for how cells adapt to receive and amplify signal information from a pulse in the presence of background GnRH. Only by considering a model with both mechanisms can their effects be separated.

Experiments involving variable-frequency application of GnRH to sheep cells in perifusion provided data for the model fit. In an intact animal, although LH responses are governed by GnRH, they are modulated by ambient levels of gonadal steroids. An advantage of perifusion systems is that they allow isolation of desired components without confounding from feedback present in vivo. In addition, these systems allow control of input for the detailed study of the dynamics of secretory patterns over time with controlled signal input. Finally, minimization of input and output signal path (length) reduces distortion of these signals. Our perifusion procedure differs from most conventional approaches in that, immediately after dissociation, sheep pituitary cells are combined with beads, loaded in a cell chamber, and perifused with media. With this regimen, regular GnRH pulses can be initiated immediately, thereby avoiding the usual prolonged incubation period without GnRH (12, 15, 26, 27). Utilizing these conditions, we have been able to maintain cells and perform experiments over the course of several days.

    MATERIAL AND METHODS
Top
Abstract
Introduction
Methods
Results
Discussion
References

Experimental Protocol

Microperifusion system. The perifusion apparatus was largely as described by Midgley et al. (19) with the following modifications.

GnRH injections. GnRH was administered as 4-min injections. Desired concentrations were generated by dilution of GnRH in experimental media. GnRH was introduced into the injection loop to warm >= 2 min, before use. Once the air valve was switched, GnRH was forced into the cell chamber at the media flow rate in a pattern approaching a square wave. At system flow rates, residence time in the system was short (~10 s) relative to the response time of LH. Thus distortion of the input signal was limited, and the input GnRH pattern was preserved.

Cleaning. At the conclusion of the experiment, cell chambers were removed from the system, disassembled, and rinsed in water. Medium lines and injection loops were flushed first with water and then with 70% ethanol before air drying.

Assay. For the described experiment, LH concentration was determined by RIA using highly purified ovine LH (oLH; NIADDK-oLH-25, AFP-5551B) for iodination and polyclonal antiserum to oLH described by Niswender et al. (23). An internal serum pool, B1371, calibrated against NIH-LH-S12, was the reference standard in all assays.

The length of the experiment (33 h), use of four chambers, and sampling at high frequencies (2.5 or 10 min) generated large numbers of samples and required two RIAs. The number of samples collected and assayed was 227, 225, 225, and 226 for chambers A, B, C, and D, respectively. Samples from chambers A and B and from chambers C and D were assayed together, so that intrachamber results could be examined within an assay. In addition, for each assay, data were collected on 75 tubes for the purposes of calibration.

Experimental Design

The experiment was conducted to observe the effect on LH response when pulse pairs were separated by different time intervals.

Treatments consisted of a pair of 4-min, 0.5 nM GnRH pulses. The first pulse was separated from the second by intervals of 5, 10, 20, or 40 min for a total of four different treatments. Each chamber received each of the four possible treatments once. Each treatment (pulse pair) was separated from the next by 2 h. During this time, no GnRH was administered.

By submitting each chamber to each of the treatments, chamber became a blocking variable, and one source of variability (chamber differences) was mitigated. Another large source of variability is elapsed time, and treatment order was chosen as the second blocking variable in an attempt to control for the effect of elapsed time. The final design was a 4 × 4 Latin square design (7), with chamber and treatment order as blocking variables and interpulse interval as the treatment. Chamber, treatment order, and treatment level were randomly assigned to variable sequence. Each level of treatment (interpulse interval) was applied to each chamber and treatment order. After randomization, the Latin square design was as follows: 20, 5, 10, and 40 min for chamber A; 10, 40, 5, and 20 min for chamber B; 5, 20, 40, and 10 min for chamber C; and 40, 10, 20, and 5 min for chamber D. The time line in Fig. 1 presents this experimental design.


View larger version (12K):
[in this window]
[in a new window]
 
Fig. 1.   Time line for each chamber illustrating timing and sequence of gonadotropin-releasing hormone (GnRH) injections. Each injection is represented by a vertical dash. Gaps between pairs of dashes represent intervals between pulse pairs (numbers on right are exact intervals in minutes); larger distances separating pairs of dashes represent 2-h intervals between pulse pairs.

Model Formulation: Physiological Considerations

The relation between decreasing GnRH pulse frequency and increasing LH pulse amplitudes was noted in ovariectomized ewes and monkeys in which endogenous GnRH was ablated and replaced with pulsatile GnRH (5, 30). Similar results were noted when rat cells were examined in perifusion systems (29). In experiments utilizing dispersed sheep cells, McIntosh and McIntosh (17) concluded that the LH release mechanism needs time to recover sensitivity between stimulations. Because these experiments were done in the absence of gonadal steroids, they suggested that the frequency of GnRH administration could modulate pituitary gonadotropin secretion independently of changes in gonadal steroid feedback.

The recovery time required between GnRH stimulations for regular, uniform LH release might be the time required for LH synthesis or diffusion to the cell membrane. An observed biphasic LH response led Bremner and Paulsen (2) to speculate on compartmentalization of LH into releasable and nonreleasable pools within each gonadotrope. They contend that the LH peak is determined by the size of the readily available, releasable LH pool, whereas the exponential decay phase of the response requires mobilizing LH from nonreleasable to releasable compartments or synthesis of new hormone. Clarke and Cummins (4) theorize that limitations in the releasable pool are responsible for the drop in LH peak height with associated increases in GnRH pulse frequency.

However, insufficient synthesis or depletion of LH stores cannot completely account for desensitization. Pituitary cells unresponsive to further stimulation by GnRH secrete LH when a 50 mM K+ stimulus (22) or phorbol myristate acetate (a potent secretagogue in pituitary cell cultures) or A-23187 (a Ca2+ ionophore, another effective secretagogue) (26) is applied. Furthermore, in perifused rat cells desensitized after 12 h of continuous stimulation, 43% of the original LH content remained (26), and >45% of the initial LH remained after desensitization of sheep cells (16). This indicates that cells desensitized to GnRH are not necessarily depleted of their gonadotrophin stores.

Knox et al. (13) formulated a theoretical model of cellular response for Dictyostelium based on the distribution of the receptor population among active and desensitized states. Receptors are desensitized during stimulation, and resensitization (recovery from desensitization) begins as soon as stimulus is removed or decreased. With lengthening pulse interval, more time exists for resensitization between pulses, a larger proportion of receptors are in the active state, and the response of the receptor system progressively rises. In the limit, the receptor is fully resensitized between pulses, no further increases in the fraction of active receptor can occur, and further increases in pulse interval have no effect (10).

With these theories on desensitizing systems taken into account, our model incorporates availability and activity of GnRH receptors and LH availability as possible causes for desensitization. A loss of receptor number, modification of the receptor to a less active state, and/or uncoupling of receptor from the effector system are processes that appear capable of leading to reduced responses in systems in which they occur. Reduced response may also result from insufficient time to replenish and/or distribute LH through the cell. In vivo, with machinery for LH synthesis intact, synthesis may compensate to some degree for LH depletion, and the magnitude of desensitization may be diminished. Depletion or inadequate synthesis is a larger concern in vitro. Thus, in this in vitro model, limited LH release and limited availability due to diffusion constraints within the cell and lack of LH synthesis must be taken into consideration.

Specifically, we model this system with three receptor states and two compartments of LH within the cell. Active, bound, and desensitized receptor states are incorporated in this model. A simpler model without the desensitized state is also considered but does not fit the observed data as well, and the difference seems statistically significant (P < 0.05). Time constants for desensitization (conversion of bound to desensitized receptor) and the reverse process, resensitization (desensitized to active receptor), were derived. In developing the model, it was presumed that all LH within the cell was contained within two pools: releasable and interior (i.e., not readily available). In the model, LH from the interior feeds releasable LH at a constant rate s. This should be appropriate if the interior pool is large enough that it remains substantially intact over the course of the experiment. This may oversimplify reality, but more elaborate models with two flow rates between interior and releasable LH (this would give an equilibrium between these pools in the absence of secretion) have been considered but do not fit appreciably better. The simpler model presented here incorporates the relevant physiology with fewer parameters and an equivalent statistical fit.

Release of LH was held to have constitutive and regulated components, with constitutive release being responsible for the baseline levels observed in perifusion. Binding of GnRH to receptor, generating bound receptor, was assumed to result in regulated release. McIntosh and McIntosh (16) did not observe LH synthesis in vitro. Synthesis was not measured in our system, and no single term is devoted to its representation; i.e., it is not explicitly represented in the model.

Model Construction and Differential Equations

This model is composed of two parts. The first part (Fig. 2, top) is concerned with exchange among three proposed GnRH receptor states: F, B, and D, representing free, bound, and desensitized receptor, respectively. In this model, all flows are irreversible. Free receptor is available to bind GnRH. Binding of GnRH to free receptor generates the bound receptor state. Eventually, bound receptor loses GnRH and returns to the active state. At some point in this latter pathway, unbound receptor is desensitized or not readily available and does not promote LH release. When the receptor recovers from the desensitized state, it is again available and free. Each reaction rate in this model is governed by a single, unknown kinetic constant. The first kinetic constant, kdf, describes conversion of desensitized to free receptor, the second, kfb, governs free to bound receptor, and the third, kbd, describes bound to desensitized receptor. This portion of the model accounts for loss of receptor activity with the desensitized receptor state. The processes governing alteration of receptor states from bound to desensitized and desensitized to free are not well understood, and the model does not attempt to distinguish between possible receptor outcomes leading to this desensitized receptor. (All possibilities are lumped into one receptor state in the interest of parsimony.)


View larger version (13K):
[in this window]
[in a new window]
 
Fig. 2.   Components and reactions involved in initial model of luteinizing hormone (LH) release (model 1). kfb, kdf, and kbd, kinetic constants (free to bound, desensitized to free, and bound to desensitized); F, B, and D, free, bound, and desensitized receptors; s, rate of reaction; R, readily releasable LH; B, bound LH.

The second part of the model describes the distribution of LH within the cell at time t. R(t) is the readily releasable portion of LH, and B(t) is the proportion of receptor in the bound state. Releasable or available LH is presumed to reside in granules near the cell plasma membrane and can exit these granules easily during exocytosis. These granules are replenished by migration of more interior granules as required. This component of the model utilizes limitations in the mass transfer of LH to account for desensitization.

The model specifies that LH release consists of two components: constitutive and regulated. The rate of regulated LH secretion from the cell at any time is a function of the current, bound receptor population and the releasable LH content of the cell. Specifically, the rate of LH release is [a1 + a2B(t)]R(t). With this form for the rate of LH release, in the absence of GnRH and bound receptor, there will be LH release at a declining baseline level, governed by the parameter a1. Then, superimposed on this basal secretion are pulses of LH released whenever GnRH binds to receptor, with this regulated component of release arising from a2 × B in Fig. 2.

Obviously, this model is simpler than physical reality. All the LH in the cell is limited to one of two locations, and at each location LH concentration is assumed to be uniform. This is surely an idealization. A more sophisticated description of the internal dynamics might include the partial differential equations common in the study of diffusion (i.e., Fick's law). Even this approach may be a simplication; the cell is likely compartmentalized with many organelles, and modeling diffusion of granules would involve multiple parameters, significantly increasing model complexity. Our interpretation of the model does not contain LH synthesis, in accordance with what has been observed in vitro (16) (but the rate s could be taken to be a rate of synthesis rather than a rate of replenishment, so the same model might be appropriate for a system with synthesis). Nor does the model incorporate a term for cell death; any gradual temporal decline in baseline must be accounted for by reduction in LH stores. Furthermore, in the model, all receptor transitions are unidirectional and irreversible. Chemical reactions are generally reversible, with rate constants describing forward and reverse reactions. However, using a single constant to describe reaction kinetics reduces the number of unknown parameters with little loss of realism, since the reverse rates are presumably very slow. To minimize dispersion, the perifusion system was designed with a short distance for the GnRH input and LH response signals to travel. To avoid additional parameters, any remaining distortion of either signal is not accounted for in the model.

The model was developed to describe aggregate behavior in a parsimonious fashion. LH compartments were assumed to be homogeneous (uniform concentration). The assumption of a totally homogeneous cell may lead to an appropriate model for average behavior, even if there is substantial cell-to-cell variation. Furthermore, kinetic and flow rate constants were held to be unknown but constant in time. The pituitary contains multiple cell types, e.g., gonadotropes, somatotropes, and lactotropes, each with variable behavior and characteristics. Although GnRH receptors have been described in multiple cell types (3), we have assumed that only gonadotropes have a significant number of GnRH receptors. Furthermore, although receptor number, total LH content, and GnRH threshold required for LH release vary among gonadotropes, generalizing the behavior of the cell in the interest of simplicity enhances our ability to fit the model without loss of critical information.

The schematic in Fig. 2 is a pictorial representation of a system of differential equations or a compartmental model. In this model the state of the pituitary at time t is described by four state variables. Three of these are receptor states, F, B, and D, and the final state variable is the amount of releasable LH, R. In receptor transitions, transfer between compartments occurs by diffusion and biochemical reactions involving GnRH. Transfer of material between LH compartments occurs by physical transport. The LH content is expressed in absolute terms, whereas each receptor state is expressed as a fraction of total amount of receptor. For convenience, we also introduce an additional state variable, C(t), which represents the total amount of extracellular LH secreted by time t. The change in this variable between two times, t1 and t2, represents the amount of LH released during the time interval (t1t2). The differential equations that govern the evolution of these state variables over time are
<FR><NU>dF(<IT>t</IT>)</NU><DE>d<IT>t</IT></DE></FR> = <FR><NU>d(free)</NU><DE>d<IT>t</IT></DE></FR> = <IT>k</IT><SUB>df</SUB>D(<IT>t</IT>) − <IT>k</IT><SUB>fb</SUB>F(<IT>t</IT>)GnRH(<IT>t</IT>)
<FR><NU>dD(<IT>t</IT>)</NU><DE>d<IT>t</IT></DE></FR> = <FR><NU>d(desensitized)</NU><DE>d<IT>t</IT></DE></FR> = <IT>k</IT><SUB>bd</SUB>B(<IT>t</IT>) − <IT>k</IT><SUB>df</SUB>D(<IT>t</IT>)
<FR><NU>dB(<IT>t</IT>)</NU><DE>d<IT>t</IT></DE></FR> = <FR><NU>d(bound)</NU><DE>d<IT>t</IT></DE></FR> = <IT>k</IT><SUB>fb</SUB>F(<IT>t</IT>)GnRH(<IT>t</IT>) − <IT>k</IT><SUB>bd</SUB>B(<IT>t</IT>)
<FR><NU>dR(<IT>t</IT>)</NU><DE>d<IT>t</IT></DE></FR> = <FR><NU>d(releasable LH)</NU><DE>d<IT>t</IT></DE></FR> = <IT>s</IT> − [<IT>a</IT><SUB>1</SUB> + <IT>a</IT><SUB>2</SUB>B(<IT>t</IT>)]R(<IT>t</IT>)
<FR><NU>dC(<IT>t</IT>)</NU><DE>d<IT>t</IT></DE></FR> = <FR><NU>d(cumulative LH)</NU><DE>d<IT>t</IT></DE></FR> = [<IT>a</IT><SUB>1</SUB> + <IT>a</IT><SUB>2</SUB>B(<IT>t</IT>)]R(<IT>t</IT>)
These equations contain three unknown kinetic constants and the unknown parameters a1, a2, and s. Parameter a1 determines the rate of basal LH secretion, LH release under conditions of no bound receptor [B(t) = 0]. Parameter a2 gives the rate gain for LH secretion in the presence of bound receptor. The constant supplement to the releasable LH pool is designated parameter s.

Although LH-containing effluent from the perifusion system is collected at a point downstream of the cells, the model predicts LH concentration just outside the cell. Predicted and observed measurements will be offset with respect to time. (Observed measurements are also subjected to dispersion, which is not incorporated in the model.) If rho  denotes the delay, i.e., the time needed for effluent to flow from the pituitary to the point of collection, LH collected in a time interval (t1t2) should be the LH secreted by the pituitary during the interval (t1 - rho , t2 - rho ). According to the model, this would be C(t2 - rho - C(t1 - rho ). Because of the slow flow rate of media, effluent is collected over time intervals of 2.5 min. Because this interval is large relative to the dynamics of the system, it is inappropriate to view LH measurements as instantaneous, and this differencing is necessary to predict LH collection from the model. Note that rho  depends on medium flow rate, tubing diameter, and length of tubing through which it travels. Because each chamber has its own experimental setup, i.e., syringe pump and assorted tubing, differences among these setups result in some variation in medium flow rate with chamber, thus necessitating incorporation of different delay parameters for each chamber.

The system of differential equations is nonlinear (due to the product of B × R in the equation for dR/dt) and cannot be solved analytically. Explicit results can only be obtained with numerical methods.

To solve these equations, a set of initial conditions is required. If the differential equations are started hours before any GnRH input, it is expected that almost all receptors will be in the free state. Therefore, initial proportions of free, desensitized, and bound receptor are set to 1.0, 0.0, and 0.0, respectively. By convention, initial cumulative LH concentration, C(0), is taken to be zero. It is difficult to make assumptions or measurements of the starting amounts of releasable LH, which can be expected to vary across chambers and experiments. So these starting values are viewed as unknown parameters.

The model just developed has eight physiological parameters (Table 1). It will be designated model 1. Model 1 contains two mechanisms to account for desensitization, a desensitized receptor and a limited, yet replenishable pool of releasable LH. To determine which mechanisms were significant, two smaller models were fit to the experimental data. The model in which s, the constant addition to the releasable pool, was removed is designated model 2 and is depicted in Fig. 3. In model 2 the releasable LH pool receives no contributions from LH transport or synthesis. In model 3 (Fig. 4) the desensitized receptor was eliminated.

                              
View this table:
[in this window]
[in a new window]
 
Table 1.   Physiological parameters


View larger version (12K):
[in this window]
[in a new window]
 
Fig. 3.   Model 2, with no replenishment of releasable LH; a1 determines rate of basal LH secretion and a2 rate gains for LH secretion in presence of bound receptor.


View larger version (12K):
[in this window]
[in a new window]
 
Fig. 4.   Model 3, containing no desensitized receptor.

Assay Parameters

The oLH concentration is measured in an RIA. The assay is dependent on a competition between radiolabeled (tagged with isotope 125I) oLH and unlabeled, or sample, LH for antibody specific for LH (23). The amount of bound, tagged LH present in the sample is inversely related to the amount of sample LH in the tube and is quantified by the number of gamma counts emitted in 1 min. For calibration, counts are observed for tubes with known concentrations of LH. Maximum binding of radioactive material occurs when a tube contains only buffer; i.e., no LH is present. These tubes are responsible for the highest radioactive counts and determination of the 100% binding point, designated b0. Nonspecific binding (NSB) tubes represent background counts that must be subtracted from each tube. The radioactive (gamma) counts emitted from samples of known oLH concentration are utilized to generate a standard calibration curve. To generate this curve, logarithm of concentration [log (concn)] is plotted vs. logit of percent binding, logit (Y/b0), in which Y is the radioactive count of an assay tube before background subtraction, b0 is the average count of the buffer control tubes before background subtraction, NSB is the average background count, and logit is the logistic function, logit (x) = log [x/(1 - x)]. The usual equation that describes the standard curve is
logit <FR><NU><IT>Y</IT></NU><DE>b<SUB>0</SUB></DE></FR> = <IT>c</IT><SUB>1</SUB> + <IT>c</IT><SUB>2</SUB> log (concn) (1)
An RIA analysis program, RIANAL (6), uses Eq. 1 to estimate LH concentration. Calibration of the constants c1 and c2 is based on counts generated from tubes of LH standards.

Diagnostic plots suggest that the standard curve (Eq. 1) is not completely linear; thus a quadratic term was introduced to account for the slight curvature. The quadratic coefficient is c3.

For this analysis, we extend this approach using the observed radioactive counts to estimate assay parameters (b0, NSB, c1, c2, and c3) and model parameters simultaneously. This approach uses all information in the assay, counts from sample and calibration tubes, to estimate all the unknown parameters. An alternative approach, treating estimates from RIANAL as the raw data, may introduce some systematic error into the concentration calculations, which the model will then attempt to predict. (We do, however, use this approach to obtain starting values for the physiological parameters in the search algorithm.) By modeling assay output (raw count data) directly, we avoid this potential bias. Starting values for the assay parameters are obtained from RIANAL. Rearranging the standard curve to solve for expected counts with the quadratic term included yields the following equation
count = <FR><NU>b<SUB>0</SUB> exp {<IT>c</IT><SUB>1</SUB> + <IT>c</IT><SUB>2</SUB> log (concn) + <IT>c</IT><SUB>3</SUB>[log (concn)]<SUP>2</SUP>} + NSB</NU><DE>1 + exp {<IT>c</IT><SUB>1</SUB> + <IT>c</IT><SUB>2</SUB> log (concn) + <IT>c</IT><SUB>3</SUB>[log (concn)]<SUP>2</SUP>}</DE></FR>
This formula relating the expected value of the count to LH concentration has five parameters: b0, NSB, c1, c2, and c3, given for reference in Table 2.

                              
View this table:
[in this window]
[in a new window]
 
Table 2.   Assay parameters

Model Evaluation

Assay parameters presents a model consisting of a set of nonlinear ordinary differential equations with a number of unknown parameters that describe LH response in an in vitro system. These parameters are estimated by maximum likelihood (28). This approach seeks estimates for which the observed data and predictions from the model agree well in a probabilistic sense.

To be more specific, the observed counts, Y1, Y2, ... , Yn, are viewed as independent normal variables. The variance of count Yi is assumed proportional to the mean, with proportionality constant gamma . Let theta  be a vector with entries listing the eight physiological parameters, five assay parameters, and the variance coefficient gamma  just introduced. The mean for any given count Yi is determined by theta  and will be denoted fi(theta ). These functions are complicated but can be calculated by numerical quadrature of the system of differential equations.

The maximum-likelihood approach seeks a value theta as an estimator of theta  that maximizes the log-likelihood function
<IT>l</IT>(&thgr;) = − <LIM><OP>∑</OP><LL><IT>i</IT>=1</LL><UL><IT>n</IT></UL></LIM> <FR><NU>[<IT>Y</IT><SUB><IT>i</IT></SUB> − <IT>f</IT><SUB><IT>i</IT></SUB>(&thgr;)]<SUP>2</SUP></NU><DE>2&ggr;<IT>f</IT><SUB><IT>i</IT></SUB>(&thgr;)</DE></FR> − <FR><NU>1</NU><DE>2</DE></FR> <LIM><OP>∑</OP><LL><IT>i</IT>=1</LL><UL><IT>n</IT></UL></LIM> log [2&pgr;&ggr;<IT>f</IT><SUB><IT>i</IT></SUB>(&thgr;)] (2)

Diagnostic work suggests that this stochastic model, in which the counts are normal with the variance proportional to the mean, is adequate. Plots of Yi - fi(theta ) vs. fi(theta ) indicate some heteroskedasticity, and the direct proportionality assumed is suggested by the Poisson distribution often used for count data (although there is extra variation not explained by the Poisson distribution). A normal probability plot of [<IT>Y</IT><SUB><IT>i</IT></SUB> − <IT>f</IT><SUB><IT>i</IT></SUB>(<A><AC>&thgr;</AC><AC>ˆ</AC></A>)]/<RAD><RCD><IT>f</IT><SUB><IT>i</IT></SUB>(<A><AC>&thgr;</AC><AC>ˆ</AC></A>)</RCD></RAD> has been used to check the assumption regarding Gaussian errors. The true distribution seems fairly symmetrical, with somewhat heavier tails (starting around the 5th and 95th percentiles) than the presumed normal distribution. So the Gaussian assumption seems tenable. Only a single data point was declared an outlier and removed: a buffer control tube that was more than three standard deviations away from the mean. An additive error structure seems natural, viewing the errors as largely due to measurement. There are other natural possibilities, but these have not been investigated. Finally, sample autocorrelations for the four chambers are 0.37, 0.30, 0.39, and 0.23. These values are small, so little gain in efficiency from a two-stage estimation scheme is envisioned. Although the model violations just noted seem fairly minor, the attained significance levels reported should be interpreted with caution.

The numerical work estimating model parameters was implemented using code from the Numerical Algorithms Group (Downers Grove, IL) library. Numerical solution of the system of differential equations was based on the Runge-Kutta-Merson method (11). Two algorithms have been utilized to maximize the likelihood function. The first, in the Numerical Algorithms Group library, is a quasi-Newton algorithm developed by Gill and Murray (9). The other is based on Powell's method (24). To try to ensure convergence to a global minimum, the optimization routine was restarted several times with random starting values. Multiple minima appear not to be a problem (restarts generally converge to near-identical values).

In the work reported here, data from runs for the four chambers were analyzed separately. This may be a bit artificial, since certain parameters should remain constant across chambers. For instance, the kinetic constants, kfb, kbd, and kdf, should be the same for the four chambers, and some chambers have assay parameters in common. There are credible reasons why the other parameters should vary from chamber to chamber; different chambers will have different pituitary cell densities and different media flow dynamics. A combined analysis of the four chambers together might be preferable but did not seem practical because of the substantial increase in the number of parameters involved.

    RESULTS
Top
Abstract
Introduction
Methods
Results
Discussion
References

Overview

LH concentrations, estimated from RIANAL, for all four chambers of the experiment are graphed against time in Fig. 5. When the shortest interval between pulses, i.e., 5 min, was applied, only one pulse was discernible. As the interval separating pulses was lengthened, the two pulses became more apparent. LH response level returned to baseline or nearly so between pulses, and two distinct pulses were observed when interpulse intervals were 20 and 40 min.


View larger version (60K):
[in this window]
[in a new window]
 
Fig. 5.   LH concentration (from RIANAL) vs. time. Vertical dashed lines indicate times of GnRH injections. All GnRH injections were 0.5 nM for 4 min. Injection timing (after hour 24) is given in Fig. 1. Top left: chamber A; top right: chamber B; bottom left: chamber C; bottom right: chamber D.

Parameter estimates under model 1 and their standard errors are listed in Table 3. The covariance matrix of the estimator theta is estimated as minus the inverse of the Hessian of the log-likelihood function, l(theta ), at theta  = theta (8). There is substantial uncertainty estimating several of the parameters. This is due in part to strong dependence between the estimators, indicated by an estimate for the correlation matrix given in Table 4 for chamber D. High correlations arise because certain linear combinations of the estimators are difficult to determine accurately. This, in turn, degrades the accuracy of parameters with significant weights in these linear combinations. The most evident example of this occurs in chamber A, where standard errors for ŝ, â1, and &Rcirc;(0) are considerable.

                              
View this table:
[in this window]
[in a new window]
 
Table 3.   Parameter estimates: model 1 

                              
View this table:
[in this window]
[in a new window]
 
Table 4.   Correlation of estimates: chamber D, model 1 

Comparison of Competing Models

Plots of the data and model estimates are given for chamber D in Fig. 6. Only observed data and the model 1 estimates are shown, since distinctions between the three model versions are not apparent visually.


View larger version (33K):
[in this window]
[in a new window]
 
Fig. 6.   Actual data and model estimates of chamber D for model 1.

To consider whether either or both of the two mechanisms used to model desensitization are significant, we examine the fit of the smaller models to determine whether they describe the data adequately. Formally, this is accomplished by computing a likelihood-ratio test statistic, determined for each model comparison and for each chamber. These test statistics are defined as
2 log &lgr; = 2<IT>l</IT>(<A><AC>&thgr;</AC><AC>ˆ</AC></A>) − 2<IT>l</IT>(<A><AC>&thgr;</AC><AC>˜</AC></A>)
where theta and <A><AC>&thgr;</AC><AC>˜</AC></A> are maximum-likelihood estimates under the full and smaller models, respectively, and l( · ) is the log-likelihood function defined in Eq. 2.

If the smaller model is adequate, 2 log lambda should have a chi 2 distribution on one degree of freedom (approximately) (28). Values of these test statistics for the comparisons of interest are given in Table 5, along with attained significance levels or P values. (As a technical aside, the P values for the model 1 vs. model 3 comparison may be too conservative: since negative values for s are untenable, a one-sided test would seem more natural, which would halve the P values reported.) Model comparison could also be considered using model selection criteria, such as Akaike's (1) information criteria. Here the change in the Akaike information criteria is (2 log lambda  - 2)/2. This approach would lead to similar conclusions.

                              
View this table:
[in this window]
[in a new window]
 
Table 5.   Testing hypothesis of significantly better fit with smaller models

The desensitized receptor is significant in all chambers (P = 0.0038, <0.0001, <0.0001, and <0.0001 for chambers A, B, C, and D, respectively). However, s, the replenishment rate for releasable LH, is not significant for chamber A but is significant for chambers B, C, and D (P = 0.8839, 0.0013, <0.0001, and 0.0002 for chambers A, B, C, and D, respectively). The insignificance of s in chamber A arises, because the model can compensate for a small value of s by adjusting other parameters.

Evidence of the significance of a term for desensitized receptor is found in graphs of bound receptor and LH release with and without desensitized receptor, shown for chamber D in Figs. 7 and 8. In Fig. 7, pulse heights for bound receptor are higher for model 3. Differences in the two plots of LH release shown in Fig. 8 are apparent at early pulse troughs and at pulse peaks when pulses are close together, i.e., when the interval between pulses is 5, 10, and 20 min.


View larger version (31K):
[in this window]
[in a new window]
 
Fig. 7.   Chamber D: proportion of bound receptor vs. time for model 1 (with desensitized receptor) and model 3 (without desensitized receptor).


View larger version (26K):
[in this window]
[in a new window]
 
Fig. 8.   Chamber D: LH release vs. time for model 1 (with desensitized receptor) and model 3 (without desensitized receptor).

In all four chambers a similar pattern was noted in regard to the kinetic parameter estimates. The relative ordering of processes from fastest to slowest is as follows: GnRH binding, desensitization, and recovery. Kinetic rates may be interpreted through their associated half-lives. The half-life for binding GnRH at 0.5 nM is log (2)/(kfb/2). Half-lives for desensitization and recovery are log (2)/kbd and log (2)/kdf, respectively. The predicted rate of binding was fast, with half-lives varying between 0.5 and 4 min. The half-lives of desensitization were approximately equal, i.e., 4 min, in all chambers. Within the four chambers the combined half-life for both stages of the recovery process (from bound to free via desensitized) was variable: 11-23 min.

    DISCUSSION AND CONCLUSIONS
Top
Abstract
Introduction
Methods
Results
Discussion
References

A model, based on theorized mechanisms of LH release from ovine pituitary cells, particularly as they relate to desensitization, was fit to experimental data generated from a perifusion system. The data consisted of radioactive counts, used to measure LH concentrations, from an experiment in which desensitization of the cells was provoked in four cell chambers. Smaller, competing models were also fit to the data. All fits utilized maximum-likelihood estimation. Likelihood-ratio tests, performed to study the relevance of the omitted parameter, revealed that the desensitized receptor was significant in all chambers. In the context of the models studied here, receptor desensitization is statistically supported by the current data. These results suggest that there is some physical mechanism that allows for alteration in the availability or activity of receptors. The LH replenishment mechanism, represented by the parameter s, was statistically significant in all but chamber A.

Although parsimony has been a major factor in our modeling process, the end result is still complicated with a total of 14 unknown parameters for model 1. Perhaps inevitably this leads to a certain amount of redundancy, in which different collections of these parameters may lead to near-identical predicted counts. As a consequence, certain parameters are difficult to estimate accurately. This is reflected in the correlation among the estimators, given for model 1, Chamber D in Table 4. For instance, a model with a small value for s and large R(0) (initial value for releasable LH) may be quite similar to a model with larger s if R(0) is reduced. Hence, it is difficult to estimate s and difficult to know whether s is necessary for an adequate model.

Estimates for the delay parameter rho  were fairly consistent among the chambers. In all four chambers, &kcirc;fb > &kcirc;bd > &kcirc;df; this ordering also defined the relative rates of the respective processes: GnRH binding, desensitization, and resensitization. Although the ordering is consistent across chambers, actual estimates for these parameters show considerable variation: for kfb there is more than an eightfold difference between the smallest and largest estimates. The half-lives of binding and desensitization were approximately 0.5-4, and 4 min, respectively, whereas the recovery half-life varied from 11 to 23 min among the chambers. Because samples were collected every 2.5 min, most of the estimated binding half-lives are below the resolution of the system. Binding is rapid, but the estimate may have substantial relative error. Together these rate estimates suggest an interpulse interval of two half-lives or 22-45 min to avoid desensitization. Examination of the data shows this to be appropriate: pulses separated by 5 and 10 min were not consistent, whereas pulses separated by 40-min intervals were distinct and fairly regular. In general, for pairs of pulses separated by 20 min, a small amount of desensitization can be detected (Fig. 5).

Possible sources of error are the assumptions and simplifications made in developing a model with a reasonably small number of parameters: idealizing the cells, receptors, reactions, and locations of LH and neglecting signal dispersion. Additional sources of error are the cells and system hardware: e.g., possible contamination, cell damage from the dissociation, cell chamber packing, and pump and/or flow rate variation. Differences among the pituitaries used in the experiment may add additional error.

Segel et al. (25) modeled systems in which communication occurs by periodic signals, such as aggregation in Dictyostelium discoideum, to demonstrate the efficiency of pulsatile signaling. Their model consisted of multiple receptor states in which conversions among receptor states were reversible, kinetic parameters were measured, and no component for diffusion within the cell existed. Li and Goldbeter (14) did a similar analysis for the pulsatile secretion of GnRH. However, for this analysis Li and Goldbeter lacked experimental values for the parameters characterizing desensitization. Instead, parameter values were taken so as to yield good agreement with the optimal pattern of GnRH stimulation in vivo as determined by Wildt et al. (30). This pattern consists of a 6-min pulse every hour. Their analysis suggests that, with these input characteristics, optimal stimulation is achieved with a relatively rapid process of desensitization followed by a slower resensitization process. Their models suggest a half-life of desensitization on the order of 1-2 min while the half-life of resensitization is ~30 min. Our half-lives for desensitization and recovery differ: 4 and 11-23 min, respectively.

The equations and parameters used by Li and Goldbeter (14) differ substantially from those presented in this formulation. The different results are likely due to the varying GnRH input patterns utilized for the fits: our model fit data with input pulses of 4 min administered at 5- to 40-min intervals to sheep in vitro, whereas Li and Goldbeter assumed a 6-min pulse applied every hour to rhesus monkeys in vivo.

To summarize our findings, conclusions regarding limited LH were difficult to draw, inasmuch as this mechanism involved a parameter that was difficult to estimate separately. In most of the chambers the parameter representing LH replenishment was significant. However, conclusions regarding the desensitized receptor were clear; it was significant in each case.

    ACKNOWLEDGEMENTS

The authors thank our colleagues Vasantha Padmanabhan, William Lemon, and Paul Favreau for making the experiment possible and the Assay and Reagents Core of the Reproductive Sciences Program for supplying assay standards and reagents.

    FOOTNOTES

This work was performed as part of the National Institute of Child Health and Human Development National Cooperative Program on Infertility Research and was supported by Grants U54 HD-29184 and R01 HD-18018.

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. §1734 solely to indicate this fact.

Address for reprint requests: K. Heinze, Pharmacokinetics, Dynamics, and Metabolism, Parke-Davis Pharmaceutical Research/Warner Lambert, 2800 Plymouth Rd., Ann Arbor, MI 48105.

Received 9 February 1998; accepted in final form 28 August 1998.

    REFERENCES
Top
Abstract
Introduction
Methods
Results
Discussion
References

1.   Akaike, H. A new look at statistical model discrimination. IEEE Trans. Automat. Control 19: 716-723, 1974.

2.   Bremner, W. J., and C. A. Paulsen. Two pools of luteinizing hormone in evidence from constant administration of the human pituitary: luteinizing hormone releasing hormone. J. Clin. Endocrinol. Metab. 39: 811-815, 1974[Medline].

3.   Childs, G. V. Division of labor among gonadotropes. Vitam. Horm. 50: 215-286, 1995[Medline].

4.   Clarke, I. J., and J. T. Cummins. GnRH pulse frequency determines LH pulse amplitude by altering the amount of releasable LH in the pituitary glands of ewes. J. Reprod. Fertil. 73: 425-431, 1985[Abstract].

5.   Clarke, I. J., J. T. Cummins, J. K. Findlay, K. J. Burman, and B. W. Doughton. Effects on plasma luteinizing hormone and follicle-stimulating hormone of varying the frequency and amplitude of gonadotropin-releasing hormone pulses in ovariectomized ewes with hypothalamo-pituitary disconnection. Neuroendocrinology 39: 214-221, 1984[Medline].

6.   Duddleson, W. G., A. R. Midgley, Jr., and G. D. Niswender. Computer program sequence for analysis and summary of radioimmunoassay data. Comp. Biomed. Res. 5: 205-217, 1972[Medline].

7.   Fisher, R. A. The arrangement of field experiments. J. Ministry Agric. 33: 503-513, 1950.

8.   Gay, D. M., and R. E. Welsh. Maximum likelihood and quasi-likelihood for nonlinear exponential family regression models. J. Am. Stat. Assoc. 83: 990-998, 1988.

9.  Gill, P. E., and W. Murray. Minimization Subject to Bounds on the Variables. National Physical Laboratory, 1976. (Tech. Rep. 72)

10.   Goldbeter, A., and Y.-X. Li. Frequency coding in intercellular communication. In: Cell to Cell Signalling: From Experiments to Theoretical Models, edited by A. Goldbeter. New York: Academic, 1989, p. 415.

11.   Hall, G., and J. M. Watt (Editors). Modern Numerical Methods for Ordinary Differential Equations. New York: Clarendon, 1976, chapt. 3, p. 59.

12.   Kamel, F., J. A. Balz, C. L. Kubajak, and V. A. Schneider. Effects of luteinizing hormone (LH)-releasing hormone pulse amplitude and frequency on LH secretion by perifused rat anterior pituitary cells. Endocrinology 120: 1644-1650, 1987[Abstract].

13.   Knox, B. E., P. N. Devreotes, A. Goldbeter, and L. A. Segel. A molecular mechanism for sensory adaptation based on ligand-induced receptor modification. Proc. Natl. Acad. Sci. USA 83: 2345-2349, 1986[Abstract].

14.   Li, Y.-X., and A. Goldbeter. Frequency specificity in intercellular communication: influence of patterns of periodic signaling on target cell response. Biophys. J. 55: 125-145, 1989[Abstract].

15.   Liu, T., and G. L. Jackson. Long-term superfusion of rat anterior pituitary cells: effects of repeated pulses of gonadotropin-releasing hormone at different doses, durations and frequencies. Endocrinology 115: 605-613, 1984[Abstract].

16.   McIntosh, J. E. A., and R. P. McIntosh. Varying the patterns and concentrations of gonadotrophin-releasing hormone stimulation does not alter the ratio of LH and FSH released from perifused sheep pituitary cells. J. Endocr. 109: 155-161, 1986[Abstract].

17.   McIntosh, R. P., and J. E. A. McIntosh. Influence of the characteristics of pulses of gonadotrophin releasing hormone on the dynamics of luteinizing hormone release from perifused sheep pituitary cells. J. Endocr. 98: 411-421, 1983[Abstract].

18.   McIntosh, R. P., and J. E. A. McIntosh. Dynamic characteristics of luteinizing hormone release from perifused sheep anterior pituitary stimulated by combined pulsatile and continuous gonadotropin-releasing hormone. Endocrinology 117: 169-179, 1985[Abstract].

19.   Midgley, A. R., Jr., R. M. Brand, P. A. Favreau, B. G. Boving, M. N. Ghazzi, V. Padmanabhan, E. Y. Young, and H. C. Cantor. Monitoring dynamic responses of perifused neuroendocrine tissues to stimuli in real time. In: Quantitative Neuroendocrinology, edited by M. L. Johnson, and J. D. Veldhuis. New York: Academic, 1994, p. 188-219.

20.   Moenter, S. M., R. M. Brand, and F. J. Karsch. Dynamics of gonadotropin-releasing hormone GnRH secretion during the GnRH surge: insights into the mechanism of GnRH surge induction. Endocrinology 130: 2978-2984, 1992[Abstract].

21.   Moenter, S. M., A. Caraty, and F. J. Karsch. The estradiol-induced surge of gonadotropin-releasing hormone in the ewe. Endocrinology 127: 1375-1384, 1990[Abstract].

22.   Naor, Z., M. Katikineni, E. Loumaye, A. G. Vela, M. L. Dufau, and K. J. Catt. Compartmentalization of luteinizing hormone pools: dynamics of gonadotropin releasing hormone action in superfused pituitary cells. Mol. Cell. Endocrinol. 27: 213-220, 1982[Medline].

23.   Niswender, G. D., A. R. Midgley, Jr., S. E. Monroe, and L. E. Reichert, Jr. Radioimmunoassay for rat luteinizing hormone using anti-ovine LH serum and radioiodinated ovine luteinizing hormone. Proc. Soc. Exp. Biol. Med. 128: 807-811, 1968.

24.   Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes: The Art of Scientific Computing. Cambridge, UK: Cambridge University Press, 1986, p. 294-301.

25.   Segel, L. A., A. Goldbeter, P. N. Devreotes, and B. E. Knox. A mechanism for exact sensory adaptation based on receptor modification. J. Theor. Biol. 120: 151-179, 1986[Medline].

26.   Smith, M. A., and W. W. Vale. Desensitization to gonadotropin-releasing hormone observed in superfused pituitary cells on Cytodex beads. Endocrinology 108: 752-759, 1981[Abstract].

27.   Smith, W. A., and M. Conn. Microaggregation of the gonadotropin-releasing hormone-receptor: relation to gonadotrope desensitization. Endocrinology 114: 553-559, 1984[Abstract].

28.   Vonesh, E. F., and V. M. Chinchilli. Linear and Nonlinear Models for the Analysis of Repeated Measures. Statistics: Textbooks and Monographs. New York: Dekker, 1997, vol. 154, p. 262-263 and 446-448.

29.   Weiss, J., and J. L. Jameson. Perifused pituitary cells as a model for studies of gonadotropin biosynthesis and secretion. Trends Endocrinol. Metab. 4: 265-270, 1993.

30.   Wildt, L., A. Häusler, T. M. P. G. Marshall, J. S. Hutchison, P. E. Belchetz, and E. Knobil. Frequency and amplitude of gonadotropin-releasing hormone stimulation and gonadotropin secretion in the rhesus monkey. Endocrinology 109: 376-385, 1981[Abstract].


Am J Physiol Endocrinol Metab 275(6):E1061-E1071
0002-9513/98 $5.00 Copyright © 1998 the American Physiological Society




This Article
Abstract
Full Text (PDF)
Alert me when this article is cited
Alert me if a correction is posted
Citation Map
Services
Email this article to a friend
Similar articles in this journal
Similar articles in PubMed
Alert me to new issues of the journal
Download to citation manager
Google Scholar
Articles by Heinze, K.
Articles by Midgley, A. R.
Articles citing this Article
PubMed
PubMed Citation
Articles by Heinze, K.
Articles by Midgley, A. R., Jr.


HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
Visit Other APS Journals Online