Mathematical modeling of gonadotropin-releasing hormone signaling

Abstract

Gonadotropin-releasing hormone (GnRH) acts via G-protein coupled receptors on pituitary gonadotropes to control reproduction. These are G_q-coupled receptors that mediate acute effects of GnRH on the exocytotic secretion of luteinizing hormone (LH) and follicle-stimulating hormone (FSH), as well as the chronic regulation of their synthesis. GnRH is secreted in short pulses and GnRH effects on its target cells are dependent upon the dynamics of these pulses. Here we overview GnRH receptors and their signaling network, placing emphasis on pulsatile signaling, and how mechanistic mathematical models and an information theoretic approach have helped further this field.

Keywords: ERK; GPCR; GnRH; Mathematical modeling; Mutual information; NFAT.

Fig. 1

A simplified GnRHR signaling network. Panel A: GnRH activates GnRHR causing a Gq/11-mediated activation of phospholipase C (PLC). This generates IP₃ which drives IP₃ receptor (IP₃R)-mediated mobilization of Ca²⁺ from intracellular stores, and diacylglycerol (DAG) which (with Ca²⁺) activates conventional PKC isozymes. GnRH increases cytoplasmic Ca²⁺ and this drives the regulated exocytotic secretion of LH and FSH from within secretory vesicles. Ca²⁺ also activates calmodulin (CaM), which activates CaM-dependent protein kinases (CaMK) and the phosphatase calcineurin (Cn), which activates the Ca²⁺-dependent transcription factor NFAT (nuclear factor of activated T-cells). GnRH also activates MAPK cascades, including the (largely PKC-mediated) activation of the Raf/MEK/ERK cascade shown. NFAT and ERK-activated transcription factors (amongst others) then act in combination to control gene expression. GnRH target genes include the gonadotropin subunits; GnRH acutely regulates the rate of vesicle fusion with the plasma membrane, and chronically regulates the gonadotropin content of these vesicles. Panels B and C: data from HeLa cells transduced to express GnRHR and also ERK2-GFP (B) or NFAT-EFP (C) that translocate from the cytoplasm to the nucleus on activation, providing live cell readouts for the Raf/MEK/ERK and CaM/Cn/NFAT activation, respectively. The data shown are the nuclear:cytoplasmic ratios (N:C) and are from an experiment in which cells received 5 min pulses of 10⁻⁷ M GnRH at 30, 60 or 120 min intervals. Note that each GnRH pulse causes nuclear translocation of each reporter and the ERK2-GFP translocation responses have more rapid on-set and off-set than the NFAT-EFP responses. Note also that with the highest pulse frequency there is insufficient time for the NFAT-EFP to return to the pre-stimulation value. Similar experiments (and experimental details) are published elsewhere (Armstrong et al., 2009a, Armstrong et al., 2009b, Armstrong et al., 2010).

Fig. 2

Simulating GnRH signaling. The GnRH signaling network has been simulated with a series of thirty-four ordinary differential equations and parameters trained on ERK2-GFP and NFAT-EFP translocation data from HeLa cells transduced with GnRHR (Tsaneva-Atanasova et al., 2012). This model was modified to add agonist-induced GnRHR internalization (and recycling), trained against data from GnRH time-course and concentration-dependence experiments in LβT2 cells (see Supplemental Data) and then used to simulate responses to GnRH pulses. The figure shows system input (square wave pulses of 10⁻⁷ M GnRH with 5 min width and 60 min period) as well as model-predicted concentrations of hormone-occupied GnRHR (HR), active PLC, cytoplasmic Ca²⁺, nuclear ppERK, nuclear Egr1 (all μM) and the nuclear fraction of NFAT (NFAT-NF). Note that the simulated upstream signals are rapid in onset and offset whereas the downstream responses (NFAT translocation and Egr1 levels) are much slower.

Fig. 3

Increasing efficiency and specificity of signaling with pulses: simulations with a minimal model. We modelled activation of an effector E1, that in turn activates two downstream effectors, E2 and E3. The traces show active effector (E1*, E2* and E3* in arbitrary units) from simulations with square wave input pulse. Activation follows Michaelis-Menten type kinetics and parameters are set for rapid activation and inactivation of E1 and E3 and for slower activation and inactivation of E2 (see parameters in Supplemental Data). Fig. 3A shows simulations with a fixed pulse width of 4 min and varied pulse period (including constant stimulation with width and period both 4 min in the top row). In addition to the time-courses (top 5 rows) we show integrated outputs as area under the curve (AUC) for the activated effectors plotted against pulse frequency (bottom row). These are non-compensated frequency-response relationships where the input integral increases in direct proportion to the frequency. For comparison, Fig. 3B shows compensated pulsatile-stimulation where any increase in frequency is offset by a reduction in pulse width so that the input integral is identical for all frequencies. Note that for the non-compensated scenario, E1* and E3* AUCs are almost directly proportional to pulse frequency because responses are rapid in onset and reversal, but slower activation and inactivation causes a non-linear relationship between pulse frequency and E2* AUC. This effect is more obvious for the compensated scenario (Fig. 3B) where the rapid E1* and E3* responses again mirror the input integral and are therefore similar at all pulse frequencies, whereas for the slower E2* responses AUC increases with pulse frequency in spite of the fact that the integrated input is identical at all frequencies (i.e. the E1* and E3* plots are effectively flat lines whereas there is an increasing monotonic relationship for E2*). Fig. 3B therefore provides a simple illustration of an integrative tracking system with rapid outputs closely mirroring the integrated input and slower responses leading to a non-linear input-output relationship. This increases efficiency (multiple brief pulses cause greater output than single long pulses) and specificity (because the same change biases signaling toward E2* as compared to E3*).

Fig. 4

Avoiding desensitization with pulses: simulations with an LβT2 cell-trained model. The data shown are concentrations of active PLC, ppERK and Egr1 from simulations of responses to 10⁻⁷ M GnRH as a constant stimulus or as 5 min pulses at 30 or 120 min intervals as indicated. The model incorporates agonist-induced receptor internalization at a rate derived from fitting wet-lab data (1×) as well as at an extremely low rate (0.001×) and with an 8-fold increased rate (8×), as indicated. The data predict receptor internalization to have a pronounced effect with constant stimulation (compare grey and blue traces in column 1) but that its effect becomes increasingly negligible with pulsatile stimulation as period increases (compare grey and blue traces in columns 2 and 3).

Fig. 5

Maximal output with sub-maximal inputs: simulations with varied feedback strength. The upper three rows show simulated Ca²⁺ responses (μM cytoplasmic Ca²⁺ concentration) for the LβT2 cell-trained model using 5 min pulses of 10⁻⁷ M GnRH at 60, 30 or 15 min and incorporating upstream negative feedback as agonist-induced receptor internalization at a rate derived from fitting wet-lab data (1×) as well as at two increased rates (8× and 16×), as indicated. The AUC of the Ca²⁺ transients is calculated (for 960 min simulations) and for each GnRHR internalization rate the condition giving the highest Ca²⁺ AUC is shown in red. Note that as internalization rate is increased, pulse-frequency-dependent desensitization becomes more evident and, as a consequence of this the greatest output is achieved with sub-maximal GnRH pulse-frequency when GnRHR internalization is set at 8× or 16×. The bottom row shows GnRH pulse frequency-response relationships from a more extensive series of simulations with GnRHR internalization varied from 0.03125× to 32× and output AUCs shown for both active PLC and Ca²⁺. Note that maximal Ca²⁺ responses only occur at sub-maximal pulse frequency when GnRHR internalization rate is 4× or greater (i.e. where pronounced desensitization of Ca²⁺ responses occurs) and that the PLC responses are maximal with constant stimulation (i.e. 12 pulses of 5 min width per hour) for all GnRHR internalization rates.

Fig. 6

Maximal output with sub-maximal input: simulations with co-operative convergent regulation of gene expression. The LβT2 cell-trained model was used to simulate GnRH signaling at various levels in the GnRHR network (PLC activity, nuclear ppERK, cytoplasmic Ca²⁺, nuclear NFAT) and also for ERK-driven transcription (ERK-DT), NFAT-driven transcription (NFAT-DT) and the situation where ERK and NFAT converge and act co-operatively to drive transcription (ERK- & NFAT-DT) as described (Tsaneva-Atanasova et al., 2012). Panel A shows output AUCs for 960 min simulations with 5 min pulses of 10⁻⁷ M GnRH at varied frequency (including constant stimulation with 12 pulses/hr) and with GnRHR internalization at a rate derived from fitting wet-lab data (1×) as well as at negligible or low rates (0.001× and 0.5×). Note that for all conditions increasing monotonic frequency-response curves are obtained except for the ERK- & NFAT-DT, for which bell-shaped frequency-response relationships are seen, even with negligible negative feedback (Fig. 6A, lower right). Panel B shows data from simulations with constant stimulation at varied GnRH concentration. As shown, increasing monotonic concentration-response curves are obtained for all outputs except for ERK- & NFAT-DT where maximal responses are predicted for sub-maximal GnRH concentration when GnRHR internalization is at 0.5× or 0.001×.

Fig. 7

Cell-cell variability and information transfer. The solid sigmoid curves in the upper cartoons illustrate population averaged responses, with individual dots representing single cell responses from which the population averages are derived. For panels A and B the population averaged data are identical but there is higher cell-cell variability in A. Consequently, frequency distribution plots shown on the left (for the stimulus concentrations indicated by the dotted lines) overlap for panel A. This creates a region of uncertainty, in that any individual cell in the area of overlap cannot “know” which stimulus concentration it has been exposed to. For panel B, cell-cell variability is much lower so the frequency-distributions do not overlap and there is no area of uncertainty. Mutual information is a statistical measure of inference quality (how reliably the system input can be inferred from the output). It is measured in Bits (with an MI of 1 indicating a system that can unambiguously distinguish two equally probable states of the environment) and would be higher in B than in A. We also illustrate the situation where the cells adapt to their environment such that the population averaged response is reduced either with a proportional reduction in cell-cell variability (A→A′) or with no change in cell-cell variability (B→B′). Note that the frequency-distributions overlap in A′ just as they do in A, and in B′ whereas they don't in B. Accordingly, the B→B′ adaptive response reduces information transfer whereas the A→A′ adaptation does not. In this scenario, consideration of population averaged responses alone can clearly deliver the wrong conclusion; if this were a hormone pre-treatment protocol one would conclude that the system has desensitized from A to A′ in spite of the fact that the quality of hormone sensing has not altered.

Fig. 8

MI as an information theoretic measure of GnRH sensing. Panels A and B show concentration and time-dependent effects of GnRH and PDBu on ERK activity in LβT2 cells, with nuclear ppERK values measured by automated fluorescence microscopy and reported in arbitrary fluorescence units (AFU, mean ± SEM, n = 3–4). The single cell measures underlying these plots were also used to calculate MI between ppERK and each of these stimuli and these values are plotted (I(ppERK; stimulus) in Bits) against time in panel C. These cells were also transduced with recombinant adenovirus for expression of an ERK-driven transcription reporter (Egr1-zsGREEN). Panel D shows the concentration-dependence of GnRH and PDBu on zsGREEN expression (in AFU, mean ± SEM, n = 3) after 360 min stimulation and the MI between zsGREEN and each of these stimuli is also shown for this time. Adapted from Garner et al., (2016).