Modeling and high-throughput experimental data uncover the mechanisms underlying Fshb gene sensitivity to gonadotropin-releasing hormone pulse frequency

Neuroendocrine control of reproduction by brain-secreted pulses of gonadotropin-releasing hormone (GnRH) represents a longstanding puzzle about extracellular signal decoding mechanisms. GnRH regulates the pituitary gonadotropin's follicle-stimulating hormone (FSH) and luteinizing hormone (LH), both of which are heterodimers specified by unique β subunits (FSHβ/LHβ). Contrary to Lhb, Fshb gene induction has a preference for low-frequency GnRH pulses. To clarify the underlying regulatory mechanisms, we developed three biologically anchored mathematical models: 1) parallel activation of Fshb inhibitory factors (e.g. inhibin α and VGF nerve growth factor-inducible), 2) activation of a signaling component with a refractory period (e.g. G protein), and 3) inactivation of a factor needed for Fshb induction (e.g. growth differentiation factor 9). Simulations with all three models recapitulated the Fshb expression levels obtained in pituitary gonadotrope cells perifused with varying GnRH pulse frequencies. Notably, simulations altering average concentration, pulse duration, and pulse frequency revealed that the apparent frequency-dependent pattern of Fshb expression in model 1 actually resulted from variations in average GnRH concentration. In contrast, models 2 and 3 showed “true” pulse frequency sensing. To resolve which components of this GnRH signal induce Fshb, we developed a high-throughput parallel experimental system. We analyzed over 4,000 samples in experiments with varying near-physiological GnRH concentrations and pulse patterns. Whereas Egr1 and Fos genes responded only to variations in average GnRH concentration, Fshb levels were sensitive to both average concentration and true pulse frequency. These results provide a foundation for understanding the role of multiple regulatory factors in modulating Fshb gene activity.

The evolution and function of the multicellular organism depends on intracellular and intercellular communication mediated by biological chemicals. One form of biological information transfer relies on sensing and responding to signaling chemical concentration. The mammalian endocrine system predominantly encodes intercellular communication as chemical concentrations (e.g. in the control of metabolic rate by thyroid hormone). Frequency-and pattern-encoded communication supports much higher-fidelity and -capacity information transfer. The specialized, high-capacity information transfer of the nervous system is associated with stimulus pattern-based electrochemical biological information transfer protocols, such as neurochemically regulated synaptic potentials and the all-ornone action potential. Although the mechanisms of information coding and decoding of many concentration-based biological chemical signals are well understood, the frequency and pulse pattern-signal information transfer protocols are largely unknown. At one key interface between the nervous system and the endocrine system, the hypothalamus controls the reproductive system largely via the release of gonadotropin-release hormone (GnRH) 5 in short discrete pulses that vary in frequency during the reproductive cycle (1). Since the discovery of this apparently frequency-sensitive, brain-endocrine information transfer system in the 1970s, the mechanisms underlying the GnRH signal information decoding by the recipient pituitary gonadotrope cell have been an area of widespread research interest.
Pulsatile release of GnRH from specialized hypothalamic neurons orchestrates the biosynthesis and secretion of the pituitary gonadotropin hormones follicle-stimulating hormone (FSH) and luteinizing hormone (LH), which in turn regulate gonadal development and steroidogenesis (1). Reproduc-tive disorders such as hypogonadotropic hypogonadism and anosmia (Kallmann syndrome) show impaired GnRH pulse secretion and subsequent abnormal FSH and LH levels; treatment with either pulsatile GnRH or gonadotropins restores fertility in those patients (for a review, see Ref. 2).
Differential regulation of FSH and LH secretion by GnRH during puberty and throughout the female menstrual cycle is characteristic of the hypothalamic-pituitary-gonadal axis (3). Higher GnRH pulse frequencies result in preferential LH secretion, whereas relatively lower GnRH pulse frequencies favor FSH production. The control of FSH and LH synthesis is linked to the transcription of their specific ␤ subunits, Fshb and Lhb. Moreover, FSH secretion is mainly constitutive, and its synthesis and release are highly correlated (4 -7). Thus, elucidating how the pituitary gonadotrope decodes the GnRH pulse signal into an Fshb transcriptional response is important due to its central role in reproductive function and dysfunction and as a paradigm for pulse-encoded signaling.
One barrier to achieving a mechanistic understanding of this regulatory system is that the specific features of the GnRH stimulus responsible for differential gene induction have not been fully considered. The importance of this aspect in studying frequency-dependent phenomena has recently been emphasized (8). In a typical experiment or simulation, changing the GnRH pulse frequency with other stimulus parameters held constant also alters the average GnRH concentration. As illustrated in Fig. 1A, high-frequency pulses (30-min cycle period) represent 4 times the average GnRH concentration per cycle period relative to lower-frequency pulses (120-min cycle period) of the same pulse duration (or width) and pulse amplitude (or height; concentration). Thus, whether the gonadotrope is responding to a specific pulse frequency or to a specific average GnRH concentration cannot be determined by only varying pulse frequency. Previous mathematical studies of GnRH signaling have largely neglected the effects of changes in average concentration delivered via altered pulse frequency on differential induction of Fshb. The various characteristics of the signal (pulse width, shape, height, and frequency) responsible for downstream responses may be complex and have non-linear relationships (9). The limited published experimental data relevant to this question about GnRH signal features (10,11) are also inconclusive in establishing true frequency sensing due to the limited sampling possible in perifusion systems and other technical considerations, such as non-linearity of receptor stimulus-response systems (12). To test the effects of altered pulse frequency at constant average concentration, the pulse duration or the pulse amplitude can be altered (Fig. 1, B and C). However, to compensate by concentration, the concentration Average GnRH concentration, which is calculated as A/T over a cycle period, is plotted for four different cycle periods (T ϭ 30, 60, 120, and 240 min; bottom). B, in the pulse duration-compensated experiment, GnRH pulses of 5-min duration are applied every 30 min (top), whereas pulses of 20-min duration are applied every 120 min (middle); pulse amplitude is invariable (A ϭ 1 nM). Subsequent to this pulse duration alteration, average GnRH concentration remains constant across the four different cycle periods (bottom). C, in the pulse amplitude-compensated experiment, GnRH pulses of 1 nM amplitude are applied every 30 min (top), whereas pulses of 4 nM amplitude are applied every 120 min (middle); pulse duration is invariable ( ϭ 5 min). Subsequent to this pulse amplitude change, average GnRH concentration remains constant across the four different cycle periods (bottom). GnRH concentration and average concentration values were expressed in arbitrary units (AU).
of the stimulus has to be in the near-linear regime of the concentration-response curve. Compensating by altering pulse duration requires testing stimulus patterns unrelated to normal physiology that may elicit complex short-term signaling regulation responses. Furthermore, the gonadotropin gene response may vary over time, and sampling that is beyond the capacity of the experimental protocols used to study these responses is needed to obtain reliable mapping of this complex parameter space. We have addressed the question of which components of GnRH pulse stimulation are responsible for Fshb induction using simple, biologically inspired mathematical models of the gonadotrope and the development of a new high-throughput experimental system.

Construction of model 1
We previously reported that G␣ s activation by GnRH promoted Lhb but suppressed Fshb gene expression in L␤T2 cells and found that this differential effect was mediated, at least in part, via the secretion of autocrine factors including inhibin ␣ (13) and VGF nerve growth factor-inducible (VGF)/neuroendocrine regulatory peptide-1 (NERP-1) (14). In contrast with Fshb gene expression, Inha and Vgf were preferentially induced by high GnRH pulse frequency. Thus, we hypothesized that secreted inhibitory factors might contribute to the lower GnRH pulse frequency preference of Fshb induction. To explore this possibility, we constructed a simple mathematical model that recreated the dependence of Fshb expression on GnRH pulsatile signal and its negative regulation by inhibitory factor(s). Model 1 comprised the GnRH stimulus, Fshb, and inhibitory factor(s) (IF), with IF exerting a negative feedback on Fshb synthesis rate (Fig. 2i); see model 1 equations under "Experimental procedures"). Although this reduced model was not a complete description of the gonadotrope network, it isolated the role of inhibitory autocrine or intracellular factors to investigate minimal requirements for the Fshb gene response to low GnRH pulse frequency.

Model 1 is consistent with experimental data
To model the conditions of the standard experiment protocol (Fig. 1A), we set pulse duration at 5 min and varied the pulse period without compensating for average concentration changes, which occur when altering frequency in this protocol. Cooperativity among the inhibitory factor(s) was required for the GnRH low-frequency (long cycle period) preference of the Fshb gene (n Ͼ 1; see Equation 4). As shown in Fig. 3, when n ϭ 1, Fshb expression increases monotonically with decreasing GnRH period. By contrast, when n ϭ 2 or n ϭ 3, Fshb gene expression shows a U-shaped GnRH frequency-response relationship with a maximum at a period of roughly 120 min. Both the GnRH stimulus for Fshb production and the inhibitory factor level increase with increasing frequency. Without cooperativity, the effect of this inhibitory factor remains proportional to the Fshb-stimulating GnRH stimulus, and a monotonic curve results. With cooperativity, the inhibitory effect increases more than the stimulatory effect at increasing frequency, generating a U-shaped response curve. n ϭ 3 was used as the cooperativity parameter in the fitting of model 1 to exper-imental data. Consistent with experimental data on the inhibitory factor VGF (14), in model simulations, the inhibitory factor showed increasing expression at decreasing GnRH pulse period (Fig. 4A), whereas Fshb showed maximal expression at about a 120-min cycle period (Fig. 4B). The model simulations showed good agreement with experimental data for Vgf and Fshb obtained using the same protocol (14) (Fig. 4).

Fshb gene induction in model 1 responds to average GnRH concentration
As noted above, changing GnRH pulse frequency in the standard protocol also alters the average GnRH concentration stimulus. To determine whether model 1 was a true frequency decoder, we carried out simulations where cycle period was varied while keeping average GnRH concentration constant by adjusting pulse duration. This compensated simulation completely eliminated the apparent frequency sensitivity seen in the standard protocol (Fig. 4). Conversely, when average GnRH concentration was varied while holding cycle period constant, a concentration sensitivity was observed (data not shown) that paralleled the results obtained with the standard protocol ( Fig.  4B (ii)). Thus, in model 1, Fshb gene activity is actually dependent on the average GnRH concentration, with the standard uncompensated protocol using GnRH pulse frequency changes as a mechanism for delivering different average GnRH concentrations. These effects can be seen more clearly in the contour plots shown in Fig. 4, A and B (iv), which summarize the results of many simulations performed at finely varied pulse durations and cycle periods. For instance, in the Fshb contour plot (Fig. 4B (iv)) the color, from dark blue to dark red, represents the level of Fshb expression. The result of simulation under any conditions can be determined from this plot. For example, the dotted horizontal line reflects a standard experiment (shown in Fig. 4B (ii)) with pulse duration held constant. The solid black lines represent coordinate variation of pulse duration and period, such that average GnRH concentration is held constant (constant /T) while the pulse frequency is changed. The level of Fshb can be seen to depend entirely on average GnRH concentration, because Fshb is unchanging along any constant average GnRH concentration (constant /T). Fshb induction has an optimum average concentration that corresponds to the concentration delivered by the 120-min GnRH pulse period in the standard 5-min pulse duration protocol.

Construction of model 2
A pulse frequency-sensitive model has been proposed by the Goldbeter group (15) to explain the effect of periodic stimuli on cellular response. The model included a receptor system capable of desensitization toward its ligand. According to this desensitization paradigm, a fast process of desensitization and a slower resensitization process were required for optimal stimulus frequency. Desensitization schemes were suggested to play a role in the GnRH system as well as in other systems (16 -18). Moreover, recent experimental studies showed that the G␣ q pathway, unlike G␣ s , is subject to rapid desensitization in response to pulse stimulation (19,20) and slowly resensitizes (21). Thus, we evaluated a model in which, in response to GnRH stimulus, a desensitizing stimulatory factor (DF) rapidly moves from a resting state (DF r ) to an active state (DF a ) and then rapidly inactivates into an inactive state (DF i ); in the absence of GnRH, DF i slowly resensitizes to the resting state. As illustrated in Fig. 2 (ii), Fshb production depends on active DF (DF a ). Model 2 equations are provided under "Experimental procedures."

Model 2 is consistent with experimental data
Computer simulations of this model were carried out at two different cycle periods of GnRH stimulus: 30 and 120 min. Pulse duration (5 min) was unchanged from model 1. As in model 1, simulations showed that Fshb expression was optimal at the long cycle period (120 min). DF a also showed maximal expression at the 120-min cycle period (Fig. 5, A and B (i)); at the long cycle period, DF a first desensitized and then was able to resensitize between pulses, whereas at the short cycle period (30 min), DF a desensitized quickly and then did not have time to resensitize. Fshb reached steady state levels after ϳ13 h. We fitted the model to previously reported experimental data (14). As shown in Fig. 5B (ii), this produced a cycle period-Fshb response curve with a maximum at the 120-min cycle period, which was consistent with the data. In contrast, the cycle period-response plot for mean DF i decreased monotonically as cycle period increased ( Fig. 5A (ii)).

Fshb gene response in model 2 is dependent on GnRH pulse frequency
To determine whether model 2 was a true frequency decoder, we carried out further simulations where cycle period was varied while keeping average GnRH concentration constant, by adjusting pulse duration. As shown in Fig. 5, A and B (iii), the plots for DF i and Fshb exhibited similar shapes as in Fig. 5A (ii), with the characteristic cycle period-Fshb response curve peaking at the 120-min cycle period. This indicated that model 2 was a true frequency decoder.
As for model 1, we examined the effects of variations in the stimulus parameters, namely cycle period and pulse duration, Mean Fshb gene expression calculated over one cycle at steady state is plotted over a range of GnRH periods for different values of n. When n ϭ 1, Fshb increases to a saturated level with decreased GnRH period. When n ϭ 2, Fshb is expressed maximally when GnRH period is 120 min. When n ϭ 3, Fshb still has the same frequency preference (120-min period), but the curve is more sharply tuned to this frequency.
on Fshb and DF i steady state values. As illustrated in Fig. 5, A and B (iv), mean DF i decreased monotonically as pulse duration and cycle period increased. By contrast, mean Fshb was maximal in a specific region of cycle period values of 100 -200 min and low pulse duration and decreased as both cycle period and pulse duration decreased. These results corroborated that Fshb gene response in model 2 was contingent upon frequency of the GnRH stimulus.

Construction of model 3
We previously identified GDF9 as a secreted autocrine factor, which induces the Fshb gene but is preferentially decreased by high-frequency GnRH pulses (22). Knockdown experiments also suggested that GDF9 contributes to the regulation of GnRH pulse frequency preference of the Fshb gene. We studied a model where GnRH initiates rapid inactivation of a DF (e.g. GDF9); in the absence of GnRH stimulation, the inactive form of DF slowly resensitizes, and DF then activates Fshb gene expression. As illustrated in Fig. 2 (iii), Fshb production depends on both DF and the GnRH stimulus. Model 3 equations are provided under "Experimental procedures."

Model 3 is consistent with experimental data
Computer simulations of this model were carried out at two different cycle periods of GnRH stimulus: 30 and 120 min. Pulse duration (5 min) was unchanged from model 1. As in model 1, simulations showed that Fshb expression was optimal at the long cycle period (120 min). DF also showed maximal expression at the 120-min cycle period (Fig. 6, A and B (i)); at the long cycle period, DF was able to resensitize between pulses, whereas at the short cycle period (30 min), DF desensitized quickly and then did not have time to resensitize. Fshb reached steady state levels after ϳ13 h. This model also was consistent with experimental data, as represented in Fig. 6B (ii). The cycle periodresponse plot for mean DF increased monotonically as cycle period increased ( Fig. 6A (ii)).

Fshb gene response in model 3 is dependent on GnRH pulse frequency
To determine whether model 3 was a true frequency decoder, we carried out further simulations where cycle period was varied while keeping average GnRH concentration constant, by adjusting pulse duration. As shown in Fig. 6, A and B (iii), the

GnRH pulse frequency decoding
plots for DF and Fshb exhibited similar shapes as in Fig. 6A (ii), with the characteristic cycle period-Fshb response curve peaking at the 120-min cycle period. This demonstrated that model 3 showed true frequency decoding.
As in model 1, we examined the effects of variations in the stimulus parameters, namely cycle period and pulse duration, on Fshb and DF steady-state values. As illustrated in Fig. 6, A and B (iv), mean DF increased monotonically as pulse duration and cycle period increased. By contrast, mean Fshb was maximal in a specific region of cycle period values of 100 -200 min and low pulse duration and decreased as both cycle period and pulse duration increased. These results corroborated that Fshb gene response in model 3 was contingent upon frequency of the GnRH stimulus.

Experimental investigation of the signal components mediating Fshb induction
The model simulations revealed that experimental data obtained using the standard protocol, which does not compensate for average concentration changes that occur when GnRH pulse frequency is altered, were consistent with mechanistically distinct models. In the case of the inhibitory factor model 1, changes in GnRH pulse frequency alter Fshb levels only because they alter the average GnRH concentration to which the cells are exposed. The research community has accepted for decades the view that Fshb is GnRH frequency-controlled. Nonetheless, the modeling results led us to the question of what components of the GnRH signal were in fact being sensed by the gonadotrope and the gene control mechanisms underlying Fshb induction. Previous attempts to evaluate GnRH frequency sensing by compensating for changes in average concentration are inconclusive due to the assumption of linearity of compensation effects, limited assays, and lack of non-gonadotropin gene response controls (11,23). The gonadotrope cell shows nonlinear GnRH concentration stimulus-response curves for both proximal signaling responses and early gene induction that have narrow linear-response ranges (12, 24 -28). Moreover, the effect of average concentration changes via pulse duration alterations on Fshb response cannot be assumed to be linear, due to the potential for rapid desensitization. Therefore, extensive study of variation of signal parameters in conjunction with comparing the responses of Fshb versus other control genes is required to definitively address this question. Due to the number of confounding variables to evaluate, conclusive pulse duration-or pulse amplitude-compensated GnRH stimulation experiments require analysis of a large number of samples. Such studies are not feasible using perifusion systems due to their limited sample throughput. To overcome this limitation, we developed a high-sample throughput cell-on-coverslip transfer protocol that can be operated either manually or using a programmable articulated joint robot arm, allowing large numbers of samples to be assayed under varying treatment protocols and time points (see "Experimental procedures"). Mechanical stimulation (such as removing media) can alter gene expression in gonadotrope cells (data not shown). We tested whether the transfer protocol limited mechanical effects by comparing the expression of the early genes Egr1 (early growth response 1) and Fos (FBJ osteosarcoma oncogene) in cells that either (i) remained in medium (no pulse, no handling control), (ii) were exposed to vehicle pulses (vehicle pulse, handling control), or (iii) were exposed to GnRH pulses. As shown in Fig. 7A, no differences in Egr1 or Fos mRNA levels were seen between the no-pulse controls (i) and the vehiclepulse samples (ii), whereas robust phasic gene responses were observed after each GnRH pulse (iii). Furthermore, Fshb induction was detected only with GnRH exposure (Fig. 7B). Similar to Fig. 7A, Egr1 and Fos exhibited robust phasic gene responses after each GnRH pulse, both at the 60-min (Fig. 7C (ii)) and 240-min cycle periods (Fig. 7C (iii)). In contrast, at the 30-min cycle period, the early gene responses remained analogous throughout the GnRH pulses ( Fig. 7C (i)).
To guide time point selection and data summarization, we evaluated the detailed temporal responses of gene changes in cells exposed to multiple pulses of GnRH. Shown in Fig. 8 is an example of the Fshb, Egr1, and Fos time trajectories obtained with low GnRH pulse frequency (120-min cycle period). Fshb levels increased steadily from pulse to pulse (Fig. 8A). In contrast, the early genes Egr1 (Fig. 8B) and Fos (Fig. 8C) showed a high level of expression 30 min after the pulse and then decreased to basal level within 60 min following the pulse.
Based on the results of a series of parameter exploration experiments that were deposited in GEO (GSE85179), we quantified gene expression levels by averaging gene copy numbers of all of the samples collected between 4 and 8 h after the initial pulse (i.e. nine time points; Fig. 9). We examined the effects of altered GnRH pulse frequency on the induction of Fshb, Egr1, and Fos by performing either pulse duration-or pulse amplitude-compensated experiments; gene responses to GnRH pulses administered every 60 min versus 120 min were com-

GnRH pulse frequency decoding
pared at constant average GnRH concentration by varying either the duration () or the amplitude (A) of the individual pulses. In pulse duration-compensated experiments, Fshb levels were reproducibly higher at the 120-min cycle period (T) than at the 60-min period for the same average GnRH concentration (Fig. 9, A and B). In contrast, Egr1 and Fos levels were insensitive to changes in GnRH pulse frequency when average GnRH concentration was unchanged. For the frequency pat-terns tested, these early genes responded only to the average GnRH concentration with no actual sensitivity to GnRH pulse frequency other than the effect of frequency changes on average GnRH concentration. These effects can be seen more clearly in the contour plots shown in Fig. 9B. At a given average GnRH concentration (where pulse duration/pulse period (/T) is constant), there was more Fshb at the 120-min cycle period than at the 60-min cycle period, whereas there was no such effect for Egr1 and Fos. Comparable results for Fshb and early gene induction were obtained in pulse amplitude-compensated experiments (Fig. 9C). Notably, whereas Fshb is responsive to GnRH pulse frequency, it is also responsive to changes in average concentration. At a given frequency, Fshb levels increase with increasing average GnRH concentration. These results suggest that both the concentration-sensitive (model 1) and frequencysensitive (models 2 and 3) circuits may be relevant to aspects of Fshb gene control by GnRH.

Discussion
The regulation of Fshb by pulsatile GnRH involves the interaction of multiple factors (for reviews, see Refs. 7 and 29). Our analysis of simple mathematical models derived from various molecular mechanisms and modeling efforts demonstrated that, although a model may be consistent with experimental GnRH pulse response data, it may not actually represent a frequency-sensing regulatory mechanism. In particular, the Fshb response seen in model 1, which is based on GnRH stimulating both Fshb and an Fshb-suppressing factor, was unaffected by compensated changes in GnRH pulse frequency that held average GnRH concentration constant and responded only to changes in the average GnRH concentration.
These simulation results made it crucial to determine experimentally whether Fshb is responding to pulse frequency itself or to average concentration changes caused by variations in pulse frequency. Studies of this fundamental question to date have been limited, in part by the potential confounds of compensating concentration experimentally. Furthermore, only a small number of samples can be feasibly studied in perifusion systems, which does not provide the resolution needed to accurately determine the effects of different GnRH stimulus features. As shown by Fletcher et al. (8), for output frequency decoding, stimulus amplitude cannot compensate for variations in stimulus frequency when a non-linearity in stimulus amplitude is present in the signaling from input to output. Additionally, we studied responses of the intact Fshb gene, thereby avoiding the potential problems of non-physiological (error bars) of three biological replicates, with markers (ϩ) denoting the individual data points. Of note, the nominal copy numbers for Fos differ from those in Fig. 7, due to the use of a different PCR primer set, as indicated under "Experimental procedures." levels of expression and responses introduced by transfected receptors or reporter constructs. We assayed early genes as controls. Our results, replicated in thousands of separate samples, demonstrated conclusively that Egr1 and Fos levels change only in response to GnRH average concentration, whereas Fshb levels are sensitive both to changes in average GnRH concentration and to the actual frequency of GnRH pulses.
Desensitization schemes were previously suggested to have a critical role in the frequency sensitivity of various systems (16 -18). Li and Goldbeter (15) discussed the application of their desensitization model to cyclic AMP signaling in Dictyostelium cells and to GnRH pulse stimulation of pituitary cells. We constructed model 2 based on a modified Goldbeter scheme and validated it experimentally. In model 2, the slow resensitization of a DF results in a refractory-like period (in a manner similar to an action potential), which suppresses further Fshb stimulation in response to high-frequency GnRH pulses. Several mechanisms of desensitization of the response to GnRH stimulation have been described, including desensitization of the G q signal (20), induction of regulator of G protein signaling protein (RGS) (43,44), induction of MAPK phosphatases (31), and down-regulation of inositol 1,4,5-trisphosphate receptors (45).
The mature murine gonadotrope L␤T2 cell line that we have studied (46,47) may not fully recapitulate the physiology of the gonadotrope cell in vivo. Pituitary cell heterogeneity and the small proportion of gonadotropes in vivo have made primary culture and in vivo gonadotrope studies challenging. The recent advent of high-throughput single-cell transcriptomic experimental and computational approaches (48 -50) may facilitate testing and refining hypotheses developed in cell line models using data obtained from gonadotropes.
Which molecular mechanisms are most important for sensing different components of the GnRH stimulus remains to be determined. Experimentally, using a standard GnRH pulse protocol that does not distinguish frequency from concentration effects, we found that knockdown of the model 1-like factor Inha changed Fshb amplitude but not Fshb response to frequency changes (22). In contrast, knockdown of Gdf9, which forms a circuit topology corresponding to model 3 that shows "true" frequency sensitivity in simulations, shifted the frequency preference of the Fshb response (22). However, further study will be required to determine which among modulators we have studied and those investigated by other groups contribute to sensing components of the GnRH signal in vitro and in vivo. Our results clarify the questions about GnRH pattern that must be resolved to understand the control mechanisms of Fshb and provide prototype models for the underlying mechanisms. Developing a predictive gonadotropin expression model, which incorporates intracellular and extracellular regulatory loci and is validated in vivo, has the potential to improve our understanding of frequency-encoded signaling and to advance the rational treatment of many diseases that are sensitive to the function of the hypothalamic-pituitary-gonadal axis.

Mathematical modeling
Definitions and concepts used in the models-For simplicity, we consider only the case of an idealized square-wave GnRH pulse having constant concentration during the pulse as well as any test-pulse exposure occurring at a single frequency having regular (constant) pulse intervals. Within these constraints, several parameters of the GnRH stimulus can be varied experimentally and in simulations: pulse amplitude, which is the concentration of GnRH during each GnRH pulse (A), pulse duration or width (), and cycle period or the interval between the onset of each pulse (T). As mentioned above, most experiments and models addressing frequency decoding have thus far considered only variations in GnRH cycle period. Notably, changes to any of these three parameters also cause changes to the average GnRH concentration (A/T).
To help separate the effects of frequency differences from the effects of average concentration differences, we use a working definition of frequency decoding. Our goal is to use the definition to assess whether a model or experimental system is actually sensitive to the frequency of the incoming stimulus or whether downstream effects depend only on average stimulus concentration.
Consider an arbitrary biological system or mathematical model that takes a square pulse stimulus, S(t), as input and leads to a response, R(t). The square pulse stimulus function, S(t), has pulse amplitude (A), pulse duration (), and cycle period (T), such that the following is true, for n ϭ 0, 1, 2, 3… The average value of S(t) over one cycle period, at steady state, is given by the following.
When GnRH pulse frequency ( ϭ 1/T) is increased (i.e. cycle period is decreased) without changing pulse amplitude (A) or pulse duration (), there is a consequential increase in average stimulus concentration. Thus, it is impossible to tell whether changes in R(t) are due to changes in the stimulus concentration or in the frequency of pulses.
To try to distinguish the effects of pulse frequency from that of concentration, we can use changes to pulse duration or pulse amplitude to compensate for the changes in average stimulus concentration as cycle period is varied. In other words, as we vary T, we also vary either A or so that A/T remains constant (see Fig. 1).
If the model is frequency-insensitive, then we expect that average response, ϽRϾ, depends only upon the average concentration of the stimulus, ͗S͘, and because ͗S͘ is constant, ͗R͘ will be constant for all T. If the model is frequency-sensitive, then we expect that when T varies while the average concentration of S(t) is constant, the output ͗R͘ is sensitive to T.
Furthermore, we differentiate here between frequency sensitivity and frequency decoding. Frequency sensitivity refers to changes in ϽRϾ as frequency varies. Frequency decoding specifically refers to the case where ϽRϾ has a maximum for a particular frequency of stimulation. At frequencies higher or lower than this preferred frequency, ϽRϾ is smaller than the maximum. A model may be frequency-sensitive, but not a frequency decoder. However, all frequency decoders are also frequency-sensitive.
However, compensating the changes in frequency by either changes in pulse amplitude or changes in pulse duration is only partly satisfactory. For each pulse, the system may have a threshold concentration, a maximal activating concentration, and a highly non-linear relationship between the two (20). Hyperbolic stimulus-response relationships are known for GnRH receptor activation (12), and inverted U concentrationresponse curves are not uncommon. However, achieving equivalent average concentration at different frequency patterns by changes of pulse duration also has the limitation that complex activating or inhibitory regulatory mechanisms may be activated at different pulse durations. Thus, concentration compensation studies of frequency sensitivity need to be interpreted cautiously, and parameter space needs thorough exploration due to non-linearity of response curves and other confounding effects. For evaluating the three prototype models, we focused on pulse width compensation at constant pulse concentration. The actual biology of the system is more complex than the models, and our subsequent detailed experiments compared both pulse concentration and pulse width compensation effects.
Model 1-Model 1 is based on GnRH activation of a factor that suppresses FSH levels (see Fig. 2), such as inhibin (13), and consists of two equations where IF represents an inhibitory factor, ␣1 is the inhibition constant, and n reflects cooperativity of the inhibitory mechanism. Because we focused on pulse duration-compensated experiments, we have set the pulse amplitude to 1. The model parameters are found in Table 1. Both IF and Fshb depend directly on S(t) and follow first-order decay. There is a single negative-feedback term, which multiplies the expression rate parameter of Fshb.
Using Equations 3 and 4, we consider the case where IF and FSH have reached a periodic steady state, with IF(0) ϭ IF(T) and FSH(0) ϭ FSH(T). Based on our definition of frequency decoding above, we want to evaluate the following quantities. Decay rates for IF and Fshb mRNAs (0.35/h for both) were derived from the measured half-lives of GnRH-induced inhibitory factor Vgf (14) and Fshb, respectively. Fshb half-life was 2.68 Ϯ 0.20 h (40) and approximately the same for Vgf. 6 A decay rate () is calculated as ϭ t1 ⁄ 2 /ln 2, where t1 ⁄ 2 is the half-life.
Model 2-Model 2, shown schematically in Fig. 2, is based on GnRH activating an intermediate that leads to Fshb induction and that spontaneously desensitizes and cannot be immediately reactivated. In response to the GnRH stimulus, DF moves rapidly from a resting state (DF r ) to a transiently active state (DF a ) and then rapidly into an inactive state (DF i   Similar to model 1, we consider the case where DF r , DF a , DF i , and FSH have reached a periodic steady state, with DF r (0) ϭ DF r (T), DF a (0) ϭ DF a (T), DF i (0) ϭ DF i (T), and FSH(0) ϭ FSH(T). As above, to check for true frequency decoding, we evaluate the average values of each component at steady state (see Equations 5 and 6). Model parameters are listed in Table 1.

Model 3 includes an intermediate (DF) that stimulates
Fshb and that is inactivated by GnRH stimulation (see Fig. 2). This corresponds to the behavior of GDF9 (22). The GnRH stimulus (S(t) as in Equation 1) causes a rapid deactivation of a desensitizing factor (DF) to a desensitized state (DFЈ), whereas in the absence of stimulus, DFЈ slowly resensitizes and returns to the active state (DF). DF has some maximal level in the absence of GnRH stimulation. Model  The model parameters are listed in Table 1. In the model, Fshb concentration depends on both DF 2 (t) and S(t) and follows first-order decay.
Again, we consider the case where all model components have reached a periodic steady state and evaluate their average values over a cycle period to check for frequency decoding (see models 1 and 2 and Equations 5 and 6).
Simulations-We used the ODE45 solver in MATLAB to simulate the models while varying patterns of the GnRH stimulus (S(t)).

Materials and cell culture
GnRH was purchased from Bachem (Torrance, CA). L␤T2 cells were obtained from Dr. Pamela Mellon (University of California San Diego, La Jolla, CA). Cells were cultured at 37°C in DMEM (Mediatech, Herndon, VA) supplemented with 10% fetal bovine serum (FBS; Gemini, Calabasas, CA) in a humidified air atmosphere of 5% CO 2 . Our L␤T2 cell line was regularly tested (every 3-6 months) for mycoplasma and interspecies contamination and authenticated by analysis of short tandem repeat DNA profiling using 27 mouse-specific microsatellite markers (Mouse Cell Check Plus, IDEXX BioResearch, Columbia, MO). Cell line authentication was achieved by comparing our cells with an original aliquot of L␤T2 cells provided by Dr. P. Mellon and used as a standard reference. Our results confirmed that our L␤T2 cells were mycoplasma-free, were of mouse origin, and had similar markers as the original cell line aliquot.

High-throughput parallel pulse experiments
750,000 authenticated L␤T2 cells were seeded on each tissue culture-treated coverslip (Thermanox, Thermo Fisher Scientific, Waltham, MA) and were grown for 2 days in DMEM supplemented with 10% FBS. Cells were incubated in DMEM supplemented with 2% charcoal-treated FBS (Gibco) and 20 mM HEPES (Mediatech, Herndon, VA) overnight before the pulse experiment. Coverslips were placed in inert coverslip racks, and pulse patterns were achieved by moving racks among GnRH, wash, and resting solutions in chambers maintained at 37°C in a water bath. The chamber solution in the water bath was DMEM supplemented with 2% charcoal-treated FBS and 20 mM HEPES. For each condition/time point, a minimum of 3 biological replicates were collected.

Quantitative real-time PCR
Coverslips were collected with forceps directly in the chamber and placed immediately in 360 l of guanidium thiocyanate RNA lysis buffer (4 M guanidium thiocyanate, 25 mM sodium citrate, pH 7, 0.5% sarcosyl (N-lauroyl sarcosine), and 0.1 M 2-mercaptoethanol). Total RNA was isolated with the "Absolutely RNA 96 microprep" kit (Agilent, Santa Clara, CA) according to the manufacturer's protocol, subjected to an ethanol precipitation to remove salts, and resuspended in elution buffer (Agilent). RNA concentrations were determined with Quant-iT RiboGreen RNA reagent (Invitrogen) using a fluorescence microplate reader (SpectraMax M3, Molecular Devices, Sunnyvale, CA). After reverse-transcription of 1 g of RNA with the Affinity Script reverse transcriptase (Agilent), samples were diluted 1:20 in molecular biology grade H 2 O (Cellgro, Manassas, VA). Later, SYBR Green qPCR assays were performed (40 cycles) in an ABI Prism 7900HT thermal cycler (Applied Biosystems, Foster City, CA) using 5 l of cDNA template and 5 l of master mix containing the specific primers for the targeted gene, Platinum TaqDNA polymerase, and the required qPCR buffer, according to manufacturer's recommendations. Three technical qPCR replicates were run for each biological replicate. Results were exported as cycle threshold (Ct) values, and Ct values of target genes were normalized to that of RPS11 in subsequent analysis. Data were expressed as arbitrary units by using the formula, E ϭ 2,500 ϫ 1.93 (rps11 Ct value Ϫ gene of interest Ct value) , where E is the expression level in arbitrary units. Primer sequences (5Ј-3Ј) were as follows: Fshb-sense, ACGAGACCGTAAGATTGCCT; Fshb-antisense, CGGCAATGTCCATCGTCGTT; Egr1-sense, ACGTCTTGGTGCCTTTTGTG; Egr1-antisense, ACATTC-TGGAGACCGAAAGC; Fos-sense, GTGTTCCTGGCAATA-GCGTG; Fos-antisense, GCAAGAAGGTGGTCGCATTC (used in Fig. 8); Fos-sense2, TTCCTGGCAATAGCGTGTTC; Fos-antisense2, TTCAGACCACCTCGACAATG (used in Fig.  7); Rps11-sense, CGTGACGAAGATGAAGATGC; Rps11antisense, GCACATTGAATCGCACAGTC. All qPCR data, which represented 36 independent experiments (with over 4,000 samples), were deposited in Gene Expression Omnibus (GEO; GSE85179).

Statistical analysis
Statistical calculations were performed using the GraphPad Prism statistical software package version 5 (GraphPad Inc., La Jolla, CA). Data were analyzed for normality followed by calculation of analysis of variance. Statistical significance was set as indicated in the figure legends with at least p Ͻ 0.05.