Firing patterns of gonadotropin-releasing hormone neurons are sculpted by their biologic state

Gonadotropin-releasing hormone (GnRH) neurons form the final pathway for the central neuronal control of fertility. GnRH is released in pulses that vary in frequency in females, helping drive hormonal changes of the reproductive cycle. In the common fertility disorder polycystic ovary syndrome (PCOS), persistent high-frequency hormone release is associated with disrupted cycles. We investigated long- and short-term action potential patterns of GnRH neurons in brain slices before and after puberty in female control and prenatally androgenized (PNA) mice, which mimic aspects of PCOS. A Monte Carlo (MC) approach was used to randomize action potential interval order. Dataset distributions were analysed to assess (i) if organization persists in GnRH neuron activity in vitro, and (ii) to determine if any organization changes with development and/or PNA treatment. GnRH neurons in adult control, but not PNA, mice produce long-term patterns different from MC distributions. Short-term patterns differ from MC distributions before puberty but become absorbed into the distributions with maturation, and the distributions narrow. These maturational changes are blunted by PNA treatment. Firing patterns of GnRH neurons in brain slices thus maintain organization dictated at least in part by the biologic status of the source and are disrupted in models of disease.


Introduction
Reproduction is controlled by interactions among the brain, anterior pituitary and gonads. Gonadotropinreleasing hormone (GnRH) neurons in the ventral diencephalon secrete GnRH near the hypothalamopituitary portal vasculature [1][2][3]. GnRH induces the anterior pituitary to synthesize and secrete luteinizing hormone (LH) and follicle-stimulating hormone (FSH) [4,5]. These hormones activate the gonadal functions of gametogenesis and steroidogenesis. Steroids provide feedback regulation to modulate the GnRH release pattern.
Inappropriate patterns of GnRH release can cause infertility [6,7]. One example is polycystic ovary syndrome (PCOS), the leading cause of infertility in women. Diagnosis requires two of the following three symptoms: absent/infrequent ovulation, elevated androgens and/or polycystic ovarian morphology [8]. Despite its high incidence (approx. 20% of women [9]), the root causes of PCOS are unknown. The present work is guided by findings that women with hyperandrogenemic PCOS, about 50% of affected women [7,10], have persistently high-frequency LH release, indicative of high-frequency GnRH release [11][12][13]. Action potential firing in neuroendocrine neurons like GnRH neurons is correlated with hormone release [14,15], thus a greater understanding of the activity of GnRH neurons in healthy versus PCOS states may identify mechanisms underlying increased GnRH/LH release frequency in PCOS.
Investigations of GnRH neuron physiology are not possible in humans. Prenatal exposure to androgens is a commonly used animal model to study this disorder. Prenatally androgenized (PNA) mice, rats, sheep and primates exhibit phenotypes that are similar to symptoms of women with PCOS, including disrupted reproductive cycles, increased androgens, and high LH pulse frequency [16][17][18][19]. Further, prenatal exposure to anti-Müllerian hormone (AMH) in mice elevates maternal neuroendocrine drive, resulting in increased prenatal androgen exposure and similar outcomes; AMH is elevated during gestation in women with PCOS [20]. While ovarian morphology aspects of PCOS in the rodent species have limitations [21], probably attributable to the polyovulatory nature of these species, the neuroendocrine phenotypes studied here are strikingly similar among species [22] and to the elevated LH pulse frequency observed in women [11][12][13]. In brain slices from adult PNA female mice [16], overall GnRH neuron firing rate is increased compared with controls, consistent with elevated LH pulse frequency observed in vivo [23]. Before puberty, however, firing rate is reduced in cells from PNA mice [23][24][25]. These observations suggest the postulate that PNA programmes a different developmental trajectory for GnRH neurons that results in different action potential firing output.
This research had two main goals. First, to address the overarching question of the validity of using GnRH neuron firing activity data to understand biological mechanisms, specifically if short-and/or long-term GnRH neuron firing activity is organized in a non-random manner. While GnRH release in vivo is clearly organized into discrete pulses [1][2][3], these occur in the context of the whole animal's physiology. Making brain slices removes both peripheral and central inputs to GnRH neurons that may contribute to this organization. Work has examined if physiologic state of the originating animal alters the fairly simple measure of mean firing rate of these cells in brain slices [20,[26][27][28][29][30]. More formal investigations of pattern organization are limited [31][32][33], however, and whether or not these patterns differ from a distribution of possible datasets with the same inter-event intervals has not been examined. The second goal was to determine if elements of pattern organization differ with reproductive state. To achieve these goals, we used Monte Carlo (MC) randomization [34] to perform additional analyses on a subset of the data from Dulka and Moenter [25]. MC is a quantitative method that generates multiple randomizations of an aspect of the original dataset (here, intervals) to create a distribution of possibilities to which the original dataset can be statistically compared. Randomization tests like Monte Carlo analyses provide power by reducing the constraints of small samples, distributions and unequal variance. These methods do not make assumptions about data distribution and can be used to address questions such as how likely it is to achieve a specific outcome.

Data used
Data were from [25]; collection of those data was approved by the Institutional Animal Care and Use Committee of the University of Michigan ( protocol 6816). We focused on female control and PNA mice at three weeks of age (3wk) and in adulthood (adult, 17-38 wks). These groups were chosen as they exhibited the greatest difference in mean firing rate, interspike intervals and burst patterning and thus made interesting points to examine if GnRH neuron firing activity is organized, and if any organization changes with development or disease model. Action currents, the currents associated with action potentials, were recorded from green-fluorescentprotein-identified neurons in acutely prepared brain slices [25]. Action currents (events) were detected and confirmed, and inter-event intervals obtained. These intervals, and the MC randomizations thereof (below), were analysed in two ways. The Cluster algorithm examines what we define as long-term patterns. Cluster was originally designed to find patterns in the release of hormones such as LH [35], which is typically sampled at 5 to 10 min intervals from the blood of experimental subjects. Electrophysiological data are collected at sub-millisecond intervals and are thus oversampled for this analysis aimed to detect long-term patterns that may be associated with hormone release. To format our data more suitably for use with the Cluster algorithm, data were divided into 120 s bins. Multiple sequential bins are compared with one another to detect peaks and nadirs in firing rate. The vary burst window (VBW) algorithm examines what we define as short-term patterns. This algorithm works on raw intervals, determining if an event can be grouped with the previous event based on a user-determined inter-event interval [25,[36][37][38]. VBW automates iterative changes in the user-defined interval, and groups of events are referred to as bursts.

Monte Carlo approach
To address whether or not short-and/or long-term GnRH neuron firing activity is organized in a nonrandom manner, we used a Monte Carlo approach in which we generate random permutations of the original data to create surrogate datasets (referred to as MC datasets below). Because we were interested in the organization of action potentials, we randomized the order of event intervals. This approach has the advantage that we do not need to have a priori knowledge about the underlying distribution of GnRH neuron firing activity with these interval characteristics. Using the surrogate data, we are then able to test the null hypothesis that the firing intervals and the resulting groupings of events observed (individual spikes, bursts and clusters) are generated by chance.

Permutation generation
To generate the permuted data for the MC approach, the ordering of the intervals between events of each original recording was randomized using a version of Durstenfeld's shuffle algorithm [39] implemented in IgorPro. Random values were obtained through the language's standard library functions using the ran2 algorithm as the underlying pseudo-random number generator [40]. This process was repeated 1000 times for each cell [41,42]. Original data and the 1000 randomized MC datasets were subjected to the two analyses described above. Limitations to MC approaches include: (i) readily available pseudorandom number generators cannot generate all possible permutations for all of the cells studied, (ii) 1000 is lower than the maximum permutations possible in the data being analysed, and (iii) the number of permutations is different among the biological groups examined because action potential frequency is altered, leading to different numbers of intervals in similar length recordings.

Long-term patterns
In vivo, GnRH pulse frequency is modulated by steroids and varies from once every several minutes to once every few hours during the female reproductive cycle [43,44]. Frequency can also change in response to experimental manipulation [45][46][47], to disease states such as PCOS [11][12][13] or to natural changes in fertility such as in seasonal breeders [48]. Peaks and nadirs within the firing rate data are of interest as they are hypothesized to be associated with neurosecretion, based on their interval being similar to that of LH pulses in vivo [31,32]. Peaks were identified using a version of the Cluster algorithm [35] implemented in IgorPro [37]. Based on previous studies of LH and GnRH pulses and GnRH neuron activity [31,32,44,49], Clusters of 2 × 2 bins (for peaks and nadirs, respectively) with a t-score of 2 were used to identify increases and decreases, and the local standard deviation was used to estimate error. Proximal variations in these parameters (cluster sizes of 1-5, t-scores of 1-3, local versus global errors) were tested but did not alter the outcomes. Cluster output parameters analysed include number of peaks in firing rate, peak frequency, amplitude and duration. These outputs and the methodology are well described and illustrated in the following reviews [50,51].

Short-term patterns
Whereas Cluster analysis gives insight into the long-term organization of the cell's activity, burst properties characterize comparatively short-term patterns within the data. Bursts are groups of action potentials separated by short (millisecond to second) intervals and are often characteristic of neurons royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 201040 3 [33,[52][53][54][55][56][57][58]. Bursts are of interest as they are correlated with neurosecretion [14,15]. Bursts were detected using burst windows from 0 to 2000 ms in 10 ms increments. Based on the distribution of number of bursts detected (see below), most analyses were confined to burst windows at 150 ms intervals from 60 ms to 810 ms as these span the burst windows at which the highest number of bursts were detected, as well as intervals used for previous studies of GnRH pattern analysis (210 ms and 360 ms, [25,53]). Burst parameters analysed include burst frequency, spikes/burst, burst duration, inter-event interval (bursts and single spikes included as events), intraburst interval and single-spike frequency.

Analyses
Cluster and VBW algorithms were run as described on each of the 1000 MC datasets, yielding frequency distributions for the output parameters (Cluster: frequency amplitude and duration of elevated firing; VBW: burst frequency, inter-event and intra-event intervals, burst duration, spikes/burst, single-spike frequency). Given the large number of datasets and their randomized nature, we assume these distributions are representative of the underlying probability distribution of data with these interval characteristics. Because these analyses account for the original data as well as the MC datasets, we refer to these as total datasets (1000 MC runs plus original data). Initially, the proportions of total datasets with a parameter greater than or equal to versus less than or equal to the corresponding parameter of the original data were calculated to test if the original falls towards the upper or lower tails of the distribution, respectively. The inclusion of values that were equal to the result from the original data in these calculations revealed a high percentage of overlap for some parameters, particularly in the number of peaks detected with Cluster (attributable to their whole number nature and limited number of possible values), as well as burst parameters at burst windows that were shorter than the majority of the observed intervals within the original data. Datasets were thus divided into three additional classifications: those with the parameter being less than that of the original, those with parameters equal to that of the original and those with parameters greater than that of the original. Because these analyses include original data plus the MC datasets, the label for the ordinate was changed from %MC datasets to %total datasets.
Additional analyses were then performed to quantify the relationship between the original data and the generated frequency distributions, and to determine if there were differences among treatment groups.

Permutation tests
Permutation tests estimate the likelihood that the order of inter-event intervals in the original data is arbitrary by comparing it to a distribution estimated by the MC datasets [59][60][61][62]. This approach is non-parametric and does not require assumptions about whether or not the distributions are normal or had equal variance. Permutation tests were thus used to determine if two sets of data have different distributions for any of the parameters quantified by Cluster or VBW analysis. Random values used in the permutation tests were generated using Python's standard library functions for generating random numbers, using Mersenne twister [63] as the underlying pseudo-random number generator. These tests were independently run using medians and means. Results were similar and we chose to use medians because not all data were normally distributed, indicating non-parametric methods are more appropriate.
In brief, to generate permutations, the difference of medians for two sets A and B was calculated and the observations from each set combined. Two new sets A' and B' (the size of A and B, respectively) were then constructed by randomly assigning each observation to one of the sets. The difference of medians for A' and B' were recorded, and this process repeated 1000 times. The number of times the difference of medians for the generated sets A' and B' were greater than or equal to the difference of medians for A and B was divided by 1000 to obtain a true two-tailed p-value.

Cell versus itself permutation test
To determine if the original data differ from the distribution of MC datasets generated from it, we used the above permutation test method, taking the original data and 1000 MC runs from each cell as the two sets to be compared. In this design, we test the null hypothesis that intervals in the original data are random.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 201040 2.9. Pairwise group permutation test To generate p-values for comparing effects of age and treatment, the above permutation method was performed between the development and treatment groups (i.e. all 3wk versus all adult; all control versus all PNA). Logical individual pairwise comparisons within this two-by-two design (e.g. adult control versus adult PNA) were then compared (i.e. <, >, =, ≤, ≥ corresponding original data) when justified. In this design, we test the null hypotheses that age and/or treatment does not affect the proportion of MC datasets that are <, >, =, ≤, ≥ the corresponding original data for any of the parameters quantified by Cluster or VBW analysis.

Binomial exact tests
To identify trends within a group for the number of cells with differences to the MC distribution as determined by the cell versus itself permutation tests, binomial exact tests were used to compare the proportions of cells in each group with p-values less than 0.05 and with p-values ≥ 0.05 to the proportion of values in groups of this size that would be expected to fall into these categories if distribution were random. Expected proportion for p < 0.05 was defined as 1/n for each group (range 7.7-14%), with the assumption being that typically one value in groups of these sizes would fall towards a tail of the distribution. This expectation was validated by performing the same cell versus itself analysis (100×) on a set of randomly generated datasets of similar value ranges; p-values were distributed fairly evenly across the possible range. Outcomes did not change when the expected proportion off cells with p-values less than 0.05 was set to 10% for all groups, thus the variation between 7.7 and 14% did not affect outcomes. Binomial exact tests are suitable for this comparison as they allow for the testing of deviations from an expected distribution of observations into two categories and are not subject to the same limitations on sample sizes as other related tests (Chi-squared test).

Do long-term firing patterns of gonadotropin-releasing hormone neurons differ from the Monte Carlo distribution?
To examine if the firing patterns generated by GnRH neurons differ from random, the cell versus itself analysis was used to compare the original data from each individual cell to its own MC datasets. Cell versus itself analyses also indicated differences between the original data and distributions for amplitude and duration of these peaks in firing rate (figure 1e). In three-week-old and adult vehicle (VEH) control mice, a higher than expected proportion of cells had original data that were different from its MC distribution for both frequency and duration of peaks in firing rate (binomial exact test, table 1). Amplitude was also different from the MC distribution in all groups except adult vehicle. The relationship of the MC datasets to the original data for Cluster-detected peak frequency, amplitude and duration is shown in figure 2. Each dot shows the percentage of total datasets from a cell that is less than, equal to or greater than the original data's value for that particular parameter. Unlike the cell versus itself permutation test, this analysis reveals not only if the original data are different from the distribution of MC values for a particular parameter, but also the direction of change. This analysis supports the above conclusion that the long-term patterns arising from the original data exhibit fewer, longer duration peaks than the MC distributions. For amplitude, however, differences revealed by binomial exact tests may be attributable to the spread of data, rather than a specific directional shift as was observed for peak frequency and duration. These observations suggest long-term GnRH neuron activity is organized to generate fewer, royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 201040 longer peaks than random activity, and that PNA treatment disrupts this organization both before and after puberty; we thus reject the null hypothesis that intervals in the original data are random for long-term patterns.

Burst window analysis of short-term firing patterns
MC datasets were next examined using the VBW algorithm. How VBW identifies bursts is shown in figure 3a and the relationship between interval duration and number of occurrences in the original datasets is in figure 3b. As burst window increases in duration, the number of bursts detected increases to a peak, then declines as bursts are merged with one another (figure 3c). The median per cent of total datasets with different relationships (<, =, >, ≤, ≥) to the original data for burst frequency is shown for all burst windows examined up to 1 s in figure 4. Original data from adult vehicle controls did not exhibit any differences from the corresponding total MC dataset distributions. This is in contrast with the other three groups, in which burst frequencies of most MC datasets were greater than in the original. This suggested the postulates that (i) burst frequency changes with typical development, and (ii) development of adult burst patterns is disrupted by PNA treatment, as adult PNA mice more closely resemble three-week-old mice from either group than adult vehicle controls.
3.4. Do short-term firing patterns of gonadotropin-releasing hormone neurons differ from the Monte Carlo distribution?
To examine these postulates in a more rigorous manner, selected burst windows (every 150 ms from 60 to 810 ms, figure 3c) were examined as above. VBW parameters of a cell's original data were first   royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 201040 observations suggest the developmental trajectory of short-term firing is different in VEH and PNA mice. This results in a split decision on the null hypothesis that intervals in the original data are random for short-term patterns, with acceptance for adult controls but rejection for other groups.

Do age and/or treatment affect short-term firing patterns of gonadotropin-releasing hormone neurons?
Comparing the effects of age and treatment on the MC distributions of these VBW parameters using pairwise group permutation tests revealed an effect of age on burst frequency, duration, spikes/burst and single-spike frequency but no effect on inter-event or intraburst interval (table 4 and figures 6-8).
These differences were observed in the physiologic range of 210-360 ms burst windows. The 60 ms burst window failed to detect any bursts, not surprising as this is shorter than most intervals between action potentials in GnRH neurons. Differences with age were largely lost at the 510 ms burst window, reappearing at 660 and 810 ms. These three burst windows are past the peak of interval durations (figure 3a) and past the peak of bursts/group as a function of burst window (figure 3c). Pairwise analysis was thus confined to 210, 360 and 510 ms burst windows. These analyses revealed differences between 3wk and adult VEH mice for burst frequency, duration, spikes per burst and single-spike frequency (supporting the first postulate), and between adult VEH and adult PNA for burst frequency and duration (supporting the second postulate, table 5). We thus reject the null hypothesis that age and/or treatment does not affect the proportion of MC datasets that are <, >, =, ≤, ≥ the corresponding original data for long-term patterns.

Discussion
Episodic, frequency-modulated GnRH release is critical for fertility. In vivo, this output is typically monitored as LH release frequency, as direct monitoring of GnRH requires access to the vasculature connecting the brain to the pituitary. In reduced preparations, such as brain slices used for electrophysiology, other parameters such as action potential firing rate become available for study. Here, we used MC analysis to examine both short-and long-term patterns in GnRH neuron firing rate. The present studies revealed that GnRH neurons have lower frequency, longer duration bouts of firing activity for both short-and long-term firing patterns than would be expected if activity was random. They further revealed that elements of these patterns change with age, and that this maturation is incomplete in PNA mice. The Cluster algorithm was used to study long-term firing patterns. Cluster was originally designed to analyse the patterns in serum hormone levels and has been used to examine LH [35] and GnRH [44]    royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 201040 Table 4. release. Cluster does not rely on hormone half-life to establish peaks, thus is appropriate for time-series data, such as firing rate, that do not have a half-life component. When the original data from individual cells were compared with that cell's total data distribution, the contribution of GnRH neuron physiology to long-term firing patterns was evident. Specifically, original data from both the 3wk and adult vehicle groups featured fewer peaks of longer duration than their respective MC datasets. Such a shift in vivo would be expected to facilitate an episodic pattern versus more continuous release. Treatment with a continuous GnRH regimen downregulates pituitary response, essentially shutting off the downstream reproductive system [64]. This phenomenon has been used to develop long-acting GnRH analogues for treatment of conditions such as precocious puberty [65] and illustrates the critical nature of the episodic GnRH release pattern. Interestingly, the long-term patterns of individual cells are not different from their MC distribution in PNA mice. This may help explain the high-frequency LH release in these mice and potentially the same phenomenon in women with PCOS. Cluster has been used to examine patterns in firing rate of GnRH neurons [29,31,32] and arcuate kisspeptin neurons [37], which are putative upstream drivers of episodic GnRH release [66,67]. In those studies, Cluster revealed increased frequency of firing peaks in castrated mice versus  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 201040 castrated mice in which homeostatic feedback was replaced via steroid implants, mirroring in vivo changes in LH pulses and validating this approach. Cluster analysis failed to reveal differences among groups in the present study when only the original data were examined. This may be attributable to the recordings being of shorter duration for the present study, precluding a rigorous estimation of long-term firing patterns from raw data alone and re-emphasizing the utility of the MC analysis.
The contribution of GnRH neuron physiology was also evident in analyses of short-term patterns, or bursting. Burst firing is thought to facilitate neurosecretion as repeated arrival of depolarizing action potentials in nerve terminals prolongs and enhances the rise in intracellular calcium required for vesicle fusion [38,[68][69][70][71]. As with long-term patterns, GnRH neurons organize their bursts to be longer and less frequent than the MC dataset distributions. The exception was adult vehicle controls, for which burst frequency did not differ from the MC dataset distributions. This is also the only group examined that exhibits typical reproductive cycles, suggesting the hypothesis that maturation of GnRH neurons, and subsequent successful reproduction, involves a shift in short-term burst organization. The median percentage of MC datasets with parameter values less than or greater than the  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 201040 original data for adult controls is repeatedly different from that of three-week-old control mice. By contrast, values for adult PNA mice are intermediate between three-week-old and adult controls. Consistent with this, pairwise comparisons of groups at the physiologic burst windows revealed differences with maturation between three-week-old and adult controls, and differences with treatment in adults between control and PNA mice. The failure of adult PNA mice to undergo a similar transition suggests PNA treatment programmes a failure of maturation of the reproductive neuroendocrine system. These observations are of interest with regard to LH release relative to sleep stage and how this relationship changes in women with PCOS and with puberty. Women with PCOS [72] and normal (male and female) pubertal subjects [73,74] exhibit similar relationships between sleep stages and LH pulse initiation. Specifically, LH pulses are more typically associated with slow-wave sleep patterns. By contrast, in normal mature mid-late follicular phase women, LH pulse initiation is more often preceded by increased wake episodes and fewer REM epochs compared with randomly selected time points. Wake and REM were not associated with LH pulse initiation in women with PCOS [72]. LH pulse initiation in the early follicular phase of the cycle is also preceded by brief awakenings [75].   royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 201040  LH pulse initiation in women with PCOS is thus more comparable to the immature state, as are GnRH neuron burst parameters in adult PNA mice in the present study. Differences in burst parameters among groups were dependent upon the burst window chosen to define bursts. Given there are more bursts at the 360 ms burst window, intentionally set near the peak, increased observation of differences could reflect an increase in observations, thus statistical power. The MC approach should minimize this caveat by repeatedly permuting the intervals. Further arguing against the number of bursts contributing to the findings, low p-values were also observed across many categories for the 210 ms burst window, in which fewer bursts are detected. Of note, the increase in number of bursts detected when moving from the 210 to 360 ms burst window was similar to the decrease in the number of bursts detected when moving from the 360 to 510 ms burst window; despite a similar number of bursts, differences were typically not observed at 510 ms. This supports the concept that physiologic ranges for GnRH burst generation are more appropriate for these analyses.
The phenomena that generate episodic GnRH release, the 'GnRH pulse generator' are still not understood, although evidence is mounting for a possible location within the arcuate nucleus kisspeptin neurons [37,66,76]. The present work confirms action potential firing patterns of GnRH neurons in coronal brain slices, which would not be subject to ongoing input from this region, exhibit both long-and short-term patterns that change with age and a disease model. These observations suggest the biologic state of these cells contributes considerably to their output patterns. Understanding the intrinsic and synaptic properties that underlie this biology and how important inputs such as a pulse generator sculpt these properties are topics of interest for future studies.
Ethics. The work in this manuscript did not require approval by any ethics board as it was entirely in silico. We did include an approval statement for the collection of the original data.