Neuronal Network Dynamics in the Posterodorsal Amygdala: Shaping Reproductive Hormone Pulsatility

neuromodulation, amygdala, mathematical model Author for correspondence: Kateryna Nechyporenko e-mail: [email protected] Krasimira Tsaneva-Atanasova e-mail: [email protected] Neuronal Network Dynamics in the Posterodorsal Amygdala: Shaping Reproductive Hormone Pulsatility Kateryna Nechyporenko 1,2, Margaritis Voliotis1,2, Xiao Feng Li3, Owen Hollings3, Deyana Ivanova4, Jamie J Walker1, Kevin T O’Byrne3, Krasimira Tsaneva-Atanasova1,2,5 1 Department of Mathematics and Statistics, University of Exeter, Stocker Rd, Exeter, EX4 4PY , United Kingdom 2 Living Systems Institute, University of Exeter, Stocker Rd, Exeter, Ex4 4PY , United Kingdom 3 Department of Women and Children’s Health, School of Life Course and Population Sciences,King’s College London, Guy’s Campus, London, SE1 1UL, United Kingdom 4 Division of Endocrinology, Diabetes, and Hypertension, Brigham and Women’s Hospital, Harvard Medical School, Boston, MA, USA 5 EPSRC Hub for Quantitative Modelling in Healthcare, University of Exeter, Stocker Rd,Exeter, Ex4 4PY , United Kingdom Normal reproductive function and fertility rely on the rhythmic secretion of gonadotropin-releasing hormone (GnRH), which is driven by the hypothalamic GnRH pulse generator. A key regulator of the GnRH pulse generator is the posterodorsal subnucleus of the medial amygdala (MePD), a brain region which is involved in processing external environmental cues, including the effect of stress. However, the neuronal pathways enabling the dynamic, stress-triggered modulation of GnRH secretion remain largely unknown. Here, we employ in-silico modelling in order to explore the impact of dynamic inputs on GnRH pulse generator activity. We introduce and analyse a mathematical model representing MePD neuronal circuits composed of GABAergic and glutamatergic neuronal populations, integrating it with our GnRH pulse generator model. Our analysis dissects the influence of excitatory and inhibitory MePD projections’ outputs on the GnRH pulse generator’s activity and reveals a functionally relevant MePD glutamatergic projection to the GnRH pulse generator, which we probe with in vivo optogenetics. Our study sheds light on how MePD neuronal dynamics affect the GnRH pulse generator activity, and offers insights into stress-related dysregulation. © The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 2ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1. Introduction The rhythmic secretion of gonadotropin-releasing hormone (GnRH) from the hypothalamus into the portal circulation is crucial in triggering the pulsatile release of gonadotropin hormones (luteinizing hormone (LH) and follicle- stimulating hormone (FSH)) from the pituitary gland [1, 2]. This dynamic process significantly contributes to the initiation of puberty and plays a pivotal role in ensuring fertility [2, 3]. The pulsatile release of GnRH is controlled by an upstream brain network known as the "GnRH pulse generator". This network is composed of KNDy neurons found in the arcuate nucleus (ARC) of the hypothalamus, which co-express kisspeptin, neurokinin B and dynorphin A [4, 5]. We have previously used a mathematical model of the KNDy network to show that periodic pulsatile activity emerges as the basal activity or external activation of the network is increased, and confirmed this modelling prediction in vivo using optogenetic stimulation of the KNDy network [6, 7]. The GnRH pulse generator is sensitive to inputs from neural centres relaying information about the psychological and physiological state of the organism. In particular, the posterodorsal subnuclei of the medial amygdala (MePD), a limbic structure responsible for emotional processing of complex external cues, is responsive to stress [8] and regulates pubertal onset [9–11] as well as LH pulsatility [12, 13]. Whilst there is limited understanding of how the MePD processes stress-related information and relays it to the GnRH pulse generator, it has been shown that optogenetic stimulation of kisspeptin neurons in the MePD increases LH pulse frequency [12]. This effect is mediated by both GABAergic and glutamatergic signalling within the MePD, since pharmacological antagonism of GABA receptors within the MePD prevents the increase in LH pulse frequency, while pharmacological antagonism of glutamate receptors within the MePD terminates the LH pulses altogether [13]. It is likely that MePD regulation of LH pulsatility is mediated, at least in part, through direct projections from the MePD to KNDy neurons within the ARC. Indeed, viral-based monosynaptic tract-tracing in mice has shown that the amygdala provides inputs to ARC KNDy neurons [14, 15]. Furthermore, although glutamatergic projections from the MePD have not yet been explored experimentally, it has been shown that stimulation of MePD GABAergic projection terminals in the ARC causes a suppression of LH pulses [16]. Taken together, the available data suggest that the MePD can modulate GnRH pulse generator activity through local kisspeptin signalling. However, how the GABA-glutamate neuronal network within the MePD integrates kisspeptin activity and relays this to the GnRH pulse generator is less clear. Here, we investigate this using a mathematical mean-field model of the MePD GABA-glutamate neuronal circuit (see Figure 1(a)). We calibrate the model using available experimental observations and analyse it to understand how changes in functional connectivity within the MePD’s neuronal network affect the system’s output. We then couple this model of the MePD with our earlier mathematical model of the GnRH pulse generator [6, 7] to study the effects of manipulating both GABAergic and glutamatergic signalling on the activity of the GnRH pulse generator. 2. Methods (a) Animals Adult Vglut-flp mice heterozygous for the allele Slc17a6 (Strain #:030212, B6;129S- Slc17a6tm1.1(flpo)Hze/J; Jackson Laboratory, Bar Harbor, ME, USA) were bred in house. Mice were genotyped by PCR using the following primers: mutant reverse, 13007- ACA CCG GCC TTA TTC CAA G; common, 34763- GAA ACG GGG GAC ATC ACT C; and wildtype reverse, 34764- GGA ATC TCA TGG TCT GTT TTG. Mice, aged7-9 weeks at time of initial surgery, were group housed unless chronically implanted with fiber optic cannulae and kept at25°C ± 1°C, 12:12 hr light/dark cycle (lights on 0700 h), with ad libitum access to food and water. All procedures were performed in accordance with UK home office regulation and approved by The King’s College London Animal Welfare and Ethical Review Body. (b) Stereotaxic injection of viral constructs and implantation of fibre optic cannula All surgical procedures were carried out under general anaesthesia using ketamine (Vetalar,100 mg/kg, i.p.; Pfizer, Sandwich, UK) and xylazine (Rompun, 10 mg/kg, i.p.; Bayer, Leverkusen, Germany) under aseptic conditions. Mice were secured in a David Kopf stereotaxic frame (Model 900, Kopf Instruments) and bilaterally ovariectomised (OVX). A midline incision of the scalp was used to expose the skull. The periosteum was removed and two small bone screws were inserted into the skull. Using a robot stereotaxic system (Neurostar, Tubingen, Germany) two windows were drilled intracranially directly above target coordinates for the posterodorsal subnucleus of the medial amygdala (MePD) (2.35 mm lateral, 1.45 mm posterior to bregma, at a depth of 5.49 mm below the skull surface) and arcuate nucleus (ARC) (0.24 mm lateral, 1.51 mm posterior to bregma, at a depth of 5.85 mm below the skull surface) obtained from the mouse brain atlas of [17]. Either AAV- CAG-FLEXFRT-ChR2(H134R)-mCherry (n = 7, 200 nL, 3 × 1012 GC/mL, Serotype: 9; Addgene, MA, USA) to express channelrhodopsin (ChR2) or control virus (n= 4, AAV-Ef1a- fDIO-mCherry, 200 mL, 3 × 1012 GC/mL; Serotype: 9; Addgene) was unilaterally injected over 10 minutes into the MePD using a 2-µL Hamilton microsyringe (Esslab, Essex, UK). The needle was left in place for 5 minutes, before withdrawing 0.2 mm, waiting another minute, and then fully withdrawing over the course of 1 minute. A fiber optic cannula (200 µm,0.39 NA, 1.25 mm ceramic ferrule, Doric Lenses, Quebec, Canada) was inserted into the brain .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 3ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . targeting the ARC. This was fixed to the skull surface and bone screws using dental cement (Super-Bond Universal Kit, Prestige Dental, UK) before suturing closed the incision. A one week recovery period was given post-surgery. After this period, the mice were handles daily to acclimatize them to the tail-top blood sampling procedure [18]. Mice were left for 4 weeks to achieve effective opsin expression before experimentation. (c) Experimental design and blood sampling for LH measurement For measurement of LH pulsatility during optogenetic stimulation, the tip of the mouse’s tail was removed with a sterile scalpel for tail-tip blood sampling [19]. The chronically implanted fiber-optic cannula was attached to a multimode fiber-optic rotary joint patch cables (Thorlabs, Ltd, Ely, UK) via a ceramic mating sleeve which allows mice to freely move while receiving blue light (473 nm wavelength). Laser (DPSS laser, Laserglow Technologies, Toronto, Canada) intensity was set to10 mW at the tip of the fibre optic patch cable. The frequency and pattern of optical stimulation was controlled by software designed in house. After 1 h acclimatization, blood samples (5µl) were collected every 5 min for 2 h. After 1 h controlled blood sampling without optical stimulation to determine baseline LH pulse frequency, mice received patterned optical stimulation (5 s on 5s off, 10-ms pulse width) at2, 5, 10 or 20 Hz for 1 h while blood sampling continued. Non-stimulation controls were performed in the same manner, with no stimulation during the second hour. The control virus injected animals received 5 Hz optical stimulation. Mice received all treatments in a random order, with at least 3 but typically 5 days between experiments. The blood samples were snap-frozen on dry ice and storing at −80°C until processed. In-house LH enzyme- linked immunosorbent assay (LH ELISA) similar to that described by [18] was used for processing of the mouse blood samples. The mouse LH standard (AFP- 5306A; NIDDK-NHPP) was purchased from Harbor-UCLA along with the primary antibody (polyclonal antibody, rabbit LH antiserum, AFP240580Rb; NIDDK-NHPP). The secondary antibody (donkey anti-rabbit IgG polyclonal antibody [horseradish peroxidase]; NA934) was from VWR International. Validation of the LH ELISA was done in accordance with the procedure described in [18] derived from protocols defined by the International Union of Pure and Applied Chemistry. Serially diluted mLH standard replicates were used to determine the linear detection range. Nonlinear regression analysis was performed using serially diluted mLH standards of known concentration to create a standard curve for interpolating the LH concentration in whole blood samples, as described previously [6]. The assay sensitivity was 0.031 ng/mL, with intra- and inter-assay coefficients of variation of 4.6% and 10.2% respectively. (d) LH pulse detection and statistical analysis LH pulses were determined using the DynPeak algorithm [20], with settings adjusted to accommodate for the high LH pulse frequency in OVX mice as previously outlined by Breen and colleagues [21]. These include using the programmes’ default parameters except the global threshold was increased to35%, the nominal peak threshold was reduced to 20 min and the 3-point peak threshold was removed. Average LH inter-pulse interval (IPI) (the period of time between two LH pulse peaks) was calculated for the 1 h control period and 1 h optogenetic stimulation period or equivalent non-stimulation period in control animals. Statistical significance was tested using a two-way repeated measures ANOVA and post-hoc Tukey test. Data was represented as mean ± SEM and p < 0.05 was considered significant. (e) Mean-Field Model of the MePD Given the established presence of GABA and glutamate neuronal populations in the MePD [22, 23], we model their interplay employing Wilson-Cowan framework [24, 25]. The framework allows us to take a system- level approach to describe the dynamic evolution of excitatory/inhibitory activity in neuronal populations due to functional interactions within a synaptically- coupled neuronal network, incorporating both cooperation and competition mechanisms. Rather than considering individual neurons within the populations, the framework gives a coarse-grained account of the mean activity in the network, which allows the investigation of putative functional interactions between the various populations as well as the overall network output, enabling coupling of our MePD network model to other neuronal networks, such as the GnRH pulse generator as represented by the KNDy network in the ARC[6, 7]. The Wilson-Cowan framework includes an inhibitory and excitatory dependant variables, both receiving an excitatory input. However, straightforward application of this framework is not sufficient to describe the MePD GABA-glutamate neuronal network to represent differential dynamics of GABA-mediated disinhibitory mechanism [26, 27]. Therefore, we have extended the original Wilson- Cowan model by incorporating an additional inhibitory population (GABA), that does not receive excitatory input (see Figure 1(b)). A key component of the Wilson-Cowan modelling framework is a sigmoid stimulus-response function ϕ(a, F, θ), which controls the mean level of activity generated in the populations at a time t [24]: ϕ(a, F, θ) = 1 1 + exp(−a(F − θ)) − 1 1 + exp(aθ) , (2.1) where F indicates the input to a given population and parameters a and θ define the value of maximum slope and half-maximum firing threshold, respectively. The use of the sigmoidal function is motivated by the fact that the majority .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 4ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . of neurons have fluctuating membrane potential near an excitability threshold, such that the probability of firing grows exponentially upon depolarisation [28]. Additionally, the constant 1/(1 + exp(aθ)) ensures that under the absence of stimulatory input to the population the firing ceases, i.e. ϕ(a, 0, θ) = 0 [24]. The input F is given by the linear sum of excitatory and inhibitory contributions, as follows: Fl(Gl, Gi, Ge) =(1 − β2)cllGl − (1 − β1)cilGi + αKp, (2.2) Fi(Gl, Gi, Ge) =(1 − β2)cliGl + (1 − α)Kp, (2.3) Fe(Gl, Gi, Ge) =(1 − β2)cleGl − (1 − β1)cieGi, (2.4) where the parameters c represent the strength of interaction from one population to another, as shown in Figure 1(b). The parameter Kp depicts the overall kisspeptin level of excitatory input to the system, which is then distributed to the populations of glutamatergic neurons and GABA interneurons in accordance with the relative glutamatergic excitation ratio parameter α ∈ [0, 1] representing the proportion of input directed to glutamatergic neuronal population. We minimised the number of inhibitory coupling strength parameters by setting self-inhibition in the GABAergic populations and functional interactions between GABAergic efferent neurons and the other two populations to zero. In the absence of data that specifically supports the inclusion of such inhibitory interactions, a model with fewer parameters is justified and easier to interpret. As one of our aims is to investigate effects of GABA and glutamate receptor antagonism following pharamcological interventions, we incorporate the terms (1 − β1) and (1 − β2), where β1 and β2 represent the proportion of suppressed functional interaction between GABA and glutamate neuronal populations, respectively. Using the stimulus-response function and the proposed interactions between the populations (Figure 1(b)) the averaged activity in the populations is governed by the following functions: fl(Gl, Gi, Ge) = − Gl + (1 − Gl)ϕ(al, Fl(Gl, Gi, Ge), θl), (2.5) fi(Gl, Gi, Ge) = − Gi + (1 − Gi)ϕ(ai, Fi(Gl, Gi, Ge), θi), (2.6) fe(Gl, Gi, Ge) = − Ge + (1 − Ge)ϕ(ae, Fe(Gl, Gi, Ge), θe), (2.7) where the dependent variables Gl, Gi, Ge represent the mean activity in the populations of glutamatergic neurons, GABA interneurons and GABAergic efferent neurons at time t, respectively. The model also includes refractory dynamics via the term (1 − G), which controls the time period during which the populations are unable to produce a signal following an activation, and its primary effect is decreasing the maximum firing rate [29]. The MePD activity is then governed by the following ordinary differential equations (ODEs): dGl dt =fl(Gl, Gi, Ge), (2.8) dGi dt =fi(Gl, Gi, Ge), (2.9) dGe dt =fe(Gl, Gi, Ge). (2.10) The presented MePD system is non-dimensional w.r.t. time [30]. To effectively couple the system, we introduce the time scaling factor d that relates arbitrary time in equations 2.8-2.10 to time in minutes: t = T · δ, (2.11) where t is the original arbitrary time, T is the new time measured in minutes, and δ is the scaling factor (min −1). The time-converted version of the model is as follows: dGl dT =δfl(Gl, Gi, Ge), (2.12) dGi dT =δfi(Gl, Gi, Ge), (2.13) dGe dT =δfe(Gl, Gi, Ge). (2.14) The parameter values for the MePD model can be found in Table 1. (f) Calculating MePD output in the mean-field model The magnitude of the mean glutamatergic and GABAergic MePD projections’ output is found in the same way, i.e. as the integral of mean activity in the respective population over the integration time period T : Mean glutamate output = 1 T Z T 0 (Gl(t)dt, (2.15) Mean GABAergic output = 1 T Z T 0 (Ge(t)dt. (2.16) To quantify the periodic output of the MePD GABA- glutamate circuit we compute the integral of the difference of mean activity in the populations of glutamatergic neurons (Gl) and GABAergic efferent neurons (Ge) over the integration time period T : Mean MePD output = 1 T Z T 0 (Gl(t) − Ge(t))dt, (2.17) where the term 1 T is the reciprocal of the time duration, allowing to normalise the output. (g) Coarse-grained Model of ARC KNDy population with MePD input Based on experimental evidence regarding MePD projections to other brain regions including the ARC [14, 15] we couple our MePD model’s output with our ARC KNDy network model [6, 7] aiming to explore the effects ofperturbations to the MePD circuit on GnRH .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 5ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Table 1. MePD model parameter values Parameter Definition Value Reference .......................................................................................................................................................................................................................................................................................... δ Temporal scaling factor (min−1) 3 [31] α relative glutamatergic excitation ratio (a.u.) 0.9 Derived Kp Excitatory input to the MePD circuit (a.u.) [0, 14] cll glutamatergic self-excitation strength (Glu → Glu) (a.u.) 18 Derived cli interaction strength Glu → GABAint (a.u.) 16 Derived cil interaction strength GABAint → Glu (a.u.) 35 Derived cle glutamatergic excitation of GABA efferents (Glu → GABAeff) (a.u.) 40 Derived cie GABAergic inhibition of GABA efferent neurons (GABAint → GABAeff) (a.u.) 25 Derived al maximum slope of Glu (a.u.) 1.3 [24] ai maximum slope of GABAint (a.u.) 2 [24] ae maximum slope of GABAeff (a.u.) 2 [24] θl half-maximum firing threshold for Glu (a.u.) 4 [24] θi half-maximum firing threshold for GABAint (a.u.) 3.7 [24] θe half-maximum firing threshold for GABAeff (a.u.) 3.7 [24] β1 GABAergic interaction suppression coefficient (a.u.) [0, 1] β2 glutamatergic interaction suppression coefficient (a.u.) [0, 1] pulse generator activity and to validate our model against experimental observations in [12, 13]. The model describing the dynamics in the KNDy neuronal network is given by the following system of ordinary differential equations: dD dT =fD(v) − dDD, (2.18) dN dT =fN (N, v) − dN N, (2.19) dv dT =fv(N, v) − dvv, (2.20) where D and N represent the concentration of Dynorphin and Neurokinin B produced by the population, and v describes the averaged firing activity in the population in spikes/min. Parameters dD,dN , and dv control the linear decay for each variable. Dynorphin and Neurokinin B secretion rates are represented by functions fD and fN , respectively, while fv describes how the firing rate changes in response to the Neurokinin B concentration and current firing rate. The neuropeptides’ secretion rates are given by the following functions: fD(v) =kD v2 v2 + K2 v,1 , (2.21) fN (N, v) =kN v2 v2 + K2 v,2 K2 D K2 D + D2 , (2.22) where kD and kN signify the neuropeptides’ secretion rates; Kv,1 and Kv,2 describe the frequency value for which the rate of Dynorphin and Neurokinin B secretion is half- maximum; and KD describes the Dynorphin concentration that results in half-maximum inhibition of NKB. In the original introduction of the KNDy model [6], the function fv can take both positive and negative values. Here we modified fv by restricting its output to be non-negative via vertical shift of the sigmoid function: fv(N, v) =v0 1 1 + exp (k(−I + m)) (2.23) (2.24) where v0 is the maximum increase of the firing rate in response to synaptic inputs I (Hz). The parameter m signifies the synaptic input level at which the increase in the firing rate becomes half-maximum. The parameter k represents the membrane’s time constant, which determines how quickly the neuron’s membrane potential changes in response to inputs. The synaptic inputs I consist of both external and internal contributions: I = I0 + pv N 2 N 2 + K2 N v + jlGl − jeGe, (2.25) where I0 stands for the basal input in the population. The excitatory effect of Neurokinin B on the firing rate is accounted via a sigmoid function with KN representing Neurokinin B’s half-maximal effect and pv control’s the strength of the connection between the neurons in the KNDy. To account for the effects of the MePD output on the KNDy, we use the terms jlGl and −jeGe, which signify the excitatory glutamatergic and inhibitory GABAergic contribution from the MePD, respectively. The parameters jl and je are presynaptic firing rate conversion parameters for the corresponding populations in the MePD. These parameters allow converting non- dimensional output from the MePD to the input to the KNDy in Hz and assign weight to the contribution and, we alter when simulating the effects of MePD projections’ stimulation. In the case of simulating the effects of neurotransmitter antagonism and kisspeptin stimulation in the MePD we set jl = je = 1. To mimic the effects of stimulating glutamatergic and GABAergic projections, we consider varied levels of MePD network excitation ( Kp) .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 6ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . and alter the presynaptic firing rate conversion parameters to change the weight distribution of the projections. In case of glutamatergic projections stimulation we increase the weight of the glutamatergic contribution and decrease the GABAergic weight contribution (j l = 1.5, je = 0.5), while for GABA projections we increase the weight of the GABAergic contribution and decrease the glutamatergic weight contribution (j l = 0.5, je = 1.5). For further details on the original KNDy model, refer to [6]. The parameters for the KNDy network can be found in Table 2. (h) Numerical simulations and bifurcation analysis. Bifurcation analysis was performed in AUTO 07-p [36], while numerical simulations were carried out in MATLAB using ode45 (Runge-Kutta method) for the MePD system and ode15s (variable-step, variable-order (VSVO) solver) for the coupled MePD-KNDy model. The codes for reproducing the analysis and simulations presented in this manuscript can be found in GitHub repository. 3. Results (a) Modelling the MePD’s GABA-glutamate circuit and its projections to the ARC The diagram depicted in Figure 1(a) provides an overview of our model describing the functional connectivity in the MePD neural circuit along with the MePD’s projections (outputs) to the ARC. In our mode we consider an excitatory population of glutamatergic neurons and an inhibitory population of GABAergic neurons given the experimentally-established presence of glutamatergic and GABAergic neurons in the MePD [22,23]. These populations of glutamatergic and GABAergic neurons interact with each other, and also extend excitatory and inhibitory connections, respectively, to a distinct neuronal population of GABAergic neurons [26, 27], which we refer to here as GABA efferent neurons. The activity in the populations of the GABA efferent neurons and glutamatergic neurons defines the MePD’s output that we consider in the model to be acting on the ARC. It has previously been shown by [13] that a kisspeptin-expressing neuronal population, found in the MePD [37], provides excitatory input to the populations of GABA interneurons and glutamatergic neurons and has been accordingly included in our modelling. Our mathematical model is based on the established Wilson-Cowan modelling framework [24, 25] and hence allows us to simulate the mean activity of the different neuronal populations; namely glutamatergic neurons (Glu), GABA interneurons (GABA int) and GABA efferent neurons (GABAeff). In this study, we couple the MePD neuronal network model to our KNDy neuronal network model [6, 7]. In previous work, we coupled a first-generation model of the MePD neuronal network with our KNDy (pulse generator) model and performed numerical simulations in order to reproduce the results of optogenetic stimulation of kisspeptin and pharmacological antagonism experiments in the MePD in [13]. The MePD neuronal network model used in [13] was based on the same framework as the model in this manuscript, but under the assumption of stationary MePD network activity and hence constant MePD output. In the present study, we investigate the proposed MePD neuronal network in more detail, taking into account the possibility of dynamic (e.g. oscillatory) MePD activity. To this end, we match the temporal activity in the circuit to the time scales of calcium activity recorded in MePD neurons [31]. Such activity is now routinely used as a proxy of mean neuronal activity, and in our case it is mediated by GABA and glutamate neuronal populations in the MePD. Full details of the model are given in Mean-Field Model of the MePD. (b) How does excitatory input decrease inhibitory tone in the MePD circuit? Optogenetic stimulation of MePD kisspeptin neurons has been shown to have a significant effect on LH pulses [12], presumably via exciting GABA-glutamate neuronal circuits and their projections to the ARC. To investigate how this effect could be relayed through the GABA-glutamate MePD neuronal network, we study the model’s behaviour under various levels of excitatory kisspeptin input. Previous analysis of the Wilson-Cowan model [24, 25, 30] has shown that oscillatory dynamics in the model can be induced via glutamatergic self-excitation and a negative-feedback loop between the populations of GABA interneurons and glutamate neurons. Accordingly, in our MePD model, we assume excitatory coupling from the glutamatergic population to the population of GABA interneurons as well as inhibitory coupling from the GABA interneurons to the glutamatergic population. To mimic experimental optogenetic stimulation of the GABA-glutamate network, we perform a bifurcation analysis using the level of MePD kisspeptin ("MePD excitation" in Figure 2) as a free parameter (for further details on the numerical methods, see Numerical simulations and bifurcation analysis). The analysis reveals that at low kisspeptin excitation, the activity of all neuronal populations considered in our model is low and exhibits stationary dynamics (Figure 2( a- c)). This could be explained by the fact that the system does not receive enough excitation to sustain oscillations, and as a result, settles in a state of low mean activity where each population behaves as a pool of independent single neuron oscillators. As excitation increases, the activity in all three populations is amplified. Numerical continuation along the stable equilibrium branch reveals that the system undergoes a change in qualitative dynamics (atKp = 1.4) due to a Hopf bifurcation (HB), giving rise to a branch of stable periodic (limit cycle) solutions. Within the parameter range where the limit cycle solutions exist, at a threshold kisspeptin level of excitatory input (K p = 1.6) the gain in the population of .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 7ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Table 2. KNDy model parameter values Parameter Definition Value Reference .......................................................................................................................................................................................................................................................................................... dD Dyn degradation rate (min−1) 0.2 [6] dN NKB degradation rate (min−1) 1 [6] dv Firing rate reset rate (min−1) 10 [32] kD Dyn signalling strength (nMmin−1) 4 [7] kN NKB signalling strength (nMmin−1) 40 [7] pv Effective strength of synaptic input (a.u.) 0.006 [7] v0 Maximum rate of neuronal activity increase (spikes min−2) 25000 [32] KD Dyn IC50 (nM) 0.3 [33] KN NKB IC50 (nM) 4 [34] Kv,1 Firing rate for half-maximal Dyn secretion (spikes min−1) 600 [35] Kv,2 Firing rate for half-maximal NKB secretion (spikes min−1) 200 [35] k Membrane’s time constant (min) 10 Fixed m Half-maximal firing rate (min−1) 0.5 Fixed I0 Basal activity (min−1) 0.14 Fixed jl presynaptic firing rate conversion parameter for glutamatergic projections (Hz) {1, 1.5, 0.5} je presynaptic firing rate conversion parameter for GABAergic projections (Hz) {1, 0.5, 1.5} GABA efferent neurons switches from positive to negative (Figure 2(c)). As a result, further increase in excitation leads to decrease in the activity of the GABA efferent neuronal population, which approaches zero with further increase in the kisspeptin excitatory input (see Figure 2(c)). This can be explained by the fact that inhibitory input from GABA interneurons to GABA efferents outweighs the excitatory input from the glutamatergic population (Figure 2(a-b)). The presence of a negative feedback loop in the system leads to a higher rate of increase in activity of the GABA interneuron population compared to the population of glutamate neurons. Therefore, by comparing the activity of the populations of glutamatergic neurons and GABA efferents at different excitation levels (Figure 2( d- e)) we observe that, overall, the MePD projections’ output would increase due to the reduction in the inhibitory GABAergic tone. In the parameter range where the model exhibits periodic behaviour, we also find that the oscillatory period decreases as we increase excitation (Figure 2( f)). The numerical range of the oscillation period confirms that the temporal activity in the model aligns well with the experimentally observed average calcium oscillation period reported in [31], where the oscillations in MePD neuronal activity occur, approximately, in the span of a minute. We further extend the bifurcation diagram (see Figure S1), identifying that the further increase in excitation leads to an exponential increase in the oscillation period and subsequent destruction of the limit cycle solution via a global homoclinic bifurcation. Following further increase in the kisspeptin level, the system enters a bistable regime, induced and destroyed via saddle- node bifurcations, followed by constant a high population activity mode. It is unknown whether MePD kisspeptin directly modulates GABA interneurons and/or glutamatergic neurons. To investigate the role of the distribution of excitatory input between the inhibitory (GABA) and excitatory (glutamate) populations in the model, we introduce a parameter that controls the relative kisspeptin excitation ratio (α). We are thus able to continue the loci of the co-dimension one Hopf and saddle-node bifurcations (electronic supplementary material, figure S1) in two- parameter space (namely, the excitatory input Kp and the relative kisspeptin excitation ratio α) (Figure 2(g)). The Hopf bifurcation and saddle-node bifurcation curves coalesce in the two-parameter space, where a co-dimension two Bogdanov-Takens (BT) point emerges. This (BT) point is also related to the appearance of a homoclinic bifurcation curve, representing a situation where the stable and unstable manifolds of a saddle equilibrium intersect, indicating the presence of complex dynamical behaviour in the system. The bifurcation curves representing the Hopf (HB), saddle node (SN) and homoclinic (HC) loci allow us to identify regions in two-parameter space characterised by different qualitative dynamics in the system; and in particular, regions where the system oscillates. We note that for the current choice of parameters (Table 1), the oscillations in the system can be induced only under the condition that the majority of excitation is directed to the excitatory population of glutamatergic neurons. To investigate how the distribution and different levels of excitation affect the output of the system within the oscillatory region, we compute a heat map depicting changes in the mean MePD projections’ output (Figure 2(g)) due to changes in activity of the populations of excitatory (glutamate) neurons and inhibitory (GABA) efferent neurons in our MePD model. For complete details on how the mean MePD projections’ output is defined and computed, see Calculating MePD output in the mean- field model. We find that an increase in the level of MePD kisspeptin excitation leads to a transition from inhibitory to excitatory MePD output. As the proportion of kisspeptin excitation to the glutamatergic population is increased, the inhibitory tone of the MePD circuit can be maintained under .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 8ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . stronger excitation in the model (Figure 2(g)). This suggests that additional (kisspeptin) excitation of the glutamatergic population may lead to an increase in GABAergic tone, depending on the functional interaction strength between the glutamatergic and GABAergic neuronal populations. Taken together, our theoretical findings suggest that the reduction in GABA efferent neuron activity amid increased excitation of the MePD neuronal circuit may be reliant on the intricate interplay between competitive inhibitory and excitatory connections to the population of GABA efferent neurons in the MePD circuit. (c) How does MePD functional network connectivity affect MePD projections’ output? The presence of GABAergic and glutamatergic neuronal populations in the MePD [22,23], as well as their importance in the modulation of GnRH pulse generator activity, has been demonstrated experimentally [13]. However, the role of their functional interactions within the circuit remains unknown. Hence, in this section we investigate how changes in the functional interaction (coupling) strength affect the dynamics of the GABA-glutamate circuit and the MePD output. The aim here is to characterise network interaction patterns associated with oscillatory dynamics and the corresponding periodic MePD excitatory/inhibitory projections’ output. The stimulatory effect of the amygdala on the GnRH pulse generator under the stimulation of kisspeptin neurons has been attributed to the activation of GABAergic interneurones, which in turn inhibit GABAergic efferent neurons [13], forming a GABA-GABA disinhibitory interaction, which is of interest in understanding functional mechanisms and dynamics in the MePD circuit. As competing excitatory and inhibitory signals counterbalance each other, we fix the interaction strength responsible for glutamate input to GABA efferent neurons and investigate how MePD projections’ output changes under the variation of kisspeptin (level) stimulation and the strength of interaction between GABA interneurons and GABA efferents (c ie) Figure 3. The mean glutamatergic projections’ output remains relatively unaffected by the strength of interaction due to the lack of GABA efferent neuronal input to the glutamate neuronal population, but its magnitude moderately increases as the excitation in the circuit increases (Figure 3(a)). Meanwhile, as the system is excited, a low strength of interaction results in an increase in the magnitude of the mean GABAergic projections’ output (Figure 3(b)), as the GABA-GABA interaction is not sufficient to induce a decrease in the activity of the population of GABA efferent neurons under the increased excitation in the MePD circuit. On the other hand, under high interaction strength, the population of GABA interneurons exerts an excessive inhibitory input to the population of GABA efferent neurons, resulting in a reduction in the GABA efferent neuronal population activity to zero (Figure 3(b)). We compute the combined MePD projections’ output as the difference in magnitude between the glutamatergic and the GABAergic projections’ outputs. We show that a low level of functional interaction strength between GABA interneurons and GABA efferents leads to a predominantly inhibitory mean MePD projections’ output, whereas a high interaction strength results in an excitatory mean MePD projection’s output that remains mostly unchanged under further excitation of the MePD circuit (Figure 3(c)). This happens as the change in the activity of the population of GABA efferent neurons is much higher compared to the population of glutamate neurons (Figure 3(a-b)). Hence, the mean MePD projections’ output is heavily dependent on the strength of GABA-GABA disinhibition. Another critical component of the system is the functional interaction strength between the populations of GABA interneurons and glutamate neurons, which facilitates oscillatory behaviour in the model by providing negative feedback between the two populations. Having fixed the kisspeptin excitation level to induce oscillatory dynamics (Kp = 2.3), we perform one parameter bifurcation analysis using the functional interaction strength between GABA interneurons and glutamate neurons (c il) as a bifurcation parameter. Our analysis reveals Hopf ( HB), saddle node (SN) and homoclinic (HC) bifurcations that we then continue in two parameters (using the functional interaction strength between glutamate neurons and GABA interneurons (c li) as a second bifurcation parameter). This enables us to investigate the qualitative dynamics of the model under the variation of GABA-glutamate interaction (Figure 4(a)). The intersections of the Hopf and saddle node curves are associated with the location of Bogdanov- Takens (BT) points, which also gives rise to the homoclinic curve, while the location where two saddle node curves meet tangentially indicates the location of a co-dimension two cusp (CP) point. Our analysis confirms that the existence of oscillatory dynamics in the MePD neuronal network intricately depends on the balance between the inhibitory and excitatory interaction strengths. Specifically, insufficient or excessive strength of functional interaction between glutamatergic neurons and the population of GABA interneurons causes loss of oscillatory dynamics, while high inhibitory strength of interaction does not prevent oscillations, but rather makes the mean MePD projections’ output more inhibitory (Figure 4(a)). Next, we consider the combined effects of the glutamatergic population’s self-excitation and its excitatory input to the population of GABA interneurons. We perform one-parameter bifurcation analysis using the strength of functional interaction between glutamate neurons and GABA interneurons (c li) as a bifurcation parameter, identifying Hopf (HB) bifurcation points (Figure 4(b)). Increasing the strength of functional interaction between glutamatergic population and GABA interneurons causes a decrease in glutamatergic population activity as the effect of the negative feedback loop is amplified, resulting in higher .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 9ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . inhibitory interaction between GABA interneurons and the glutamatergic population. We then continue the detected (HB) points in two parameters (glutamatergic self-excitation strength (cll) and strength of functional interaction between glutamate neurons and GABA interneurons (c li)) to define the oscillatory region in two-parameter space (Figure 4(b)). We find that higher levels of self-excitation require higher strength of functional interaction between glutamate neurons and GABA interneurons in order to give rise to oscillatory dynamics. Increasing the strength of glutamatergic functional interaction switches the mean MePD projections’ output from negative to positive. The above analysis indicates that the intricate excitatory/inhibitory MePD projections’ balance is heavily dependent on the MePD circuit’s functional interactions. It also suggests that coordinated changes in the interaction strengths between the circuit’s neuronal populations may be a critical regulator of the MePD projections’ output, and thus its modulatory effect on the GnRH pulse generator. (d) MePD projections’ dynamic modulation of GnRH pulse generator activity Having analysed the MePD GABA-glutamate neuronal network behaviour and the effects of different model parameters on the mean MePD projections’ output, we now investigate our coupled MePD-KNDy network model (see Coarse-grained Model of ARC KNDy population with MePD input for full model details). The aim here is to characterising the differential effects of MePD dynamic projections’ output on GnRH pulsatility. In essence, coupling the MePD and KNDy models results in feeding external periodic input from the MePD neuronal network with the KNDy neuronal network (a.k.a. GnRH pulse generator). As the timescales of the two network models are significantly different (MePD neuronal network evolves on a timescale of seconds while the KNDy neuronal network operates on a timescale of minutes), in the coupled model there is more than one frequency found in the periodic trajectory which now evolves on a torus, i.e. the limit cycle becomes a limit torus solution (Figure 5(a)). In the case when the extended system has two incommensurate frequencies, the trajectory is no longer closed, leading to quasi-periodic dynamics. This is not surprising as it is well established that in response to periodic input, relaxation oscillators such as the KNDy network model [6] can exhibit complex dynamics, like quasi-periodicity [38, 39]. Next, we validate the coupled model by reproducing in vivo experiments where the effect of optogenetic stimulation of MePD kisspeptin neurons on LH pulse frequency was investigated [12]. To simulate the effects of optogenetic stimulation, we increase the kisspeptin level of excitation within the MePD, which results in activation of the GABA-GABA component of the MePD neuronal circuit, resulting in a decrease in the activity of the population of GABA efferent neurons’(Figure 5(b)). Consequently, the reduction of the inhibitory tone in the MePD projections’ output under increased excitation promotes a decrease in the interpulse interval (IPI) in the KNDy system (Figure 5(b)). To investigate the effects of suppression of GABAergic interaction strength on the system’s dynamics, we compute a two-parameter bifurcation diagram for a range of MePD excitation (K p) and GABA functional interaction strength suppression coefficient β1 (see Figure 1(b)), the latter describing the strength of GABA receptor antagonism (Figure 5(c)). We find that complete suppression of GABAergic interaction leads to the loss of oscillatory dynamics in the MePD circuit via its effect on the negative feedback loop between the populations of GABA interneurons and glutamatergic neurons. It is common that the effect of pharmacological blockers is modelled by complete suppression of functional interactions [13], but in reality only partial suppression may occur. Here, we show that partial blocking of GABAergic interaction is sufficient to decrease the MePD projections’ output under increased excitation (Figure 5(c)). On the other hand, reduction of GABAergic functional interaction strength leads to a decrease in the inhibitory coupling between the population of GABA interneurons and the glutamatergic population and the population of GABA efferent neurons, hence increasing both glutamatergic and GABAergic tone in the MePD (Figure 5(d)). However, the difference between the inhibitory and excitatory tone remains relatively constant as before the suppression, resulting in an unperturbed KNDy network interpulse interval. Combining suppression of GABAergic interactions with increased excitation in the MePD also increases the activity in the populations of glutamatergic neurons and GABAergic efferent neurons, but also amplifies the difference between inhibitory and excitatory mean MePD projections’ output (Figure 5(e)). As a result, the MePD projections’ input to KNDy becomes more inhibitory, leading to an increase in the KNDy network interpulse interval. Similarly, to mimic the effects of a glutamate receptor antagonist, we decrease the strength of glutamatergic interactions in the model and observe that oscillatory dynamics rapidly cease (Figure 5(f)). We also observe that under partial suppression of glutamatergic interaction strength and increased excitation, there is a transition from inhibitory to excitatory tone in the system, while complete abolition of glutamatergic interactions causes the mean MePD output to be exclusively excitatory for all levels of excitation (Figure 5(f)). To preserve inhibitory output from the MePD under low levels of excitation, we consider a partial suppression of the functional interaction strength associated with the population of glutamatergic neurons. For lower levels of excitation (still enabling oscillatory dynamics), moderate suppression of glutamatergic interaction strength results in a loss of oscillatory MePD dynamics, while keeping the mean MePD output relatively constant (and inhibitory); hence the extended system dynamics transitions from evolving on a limit torus to a limit cycle while the IPI of the .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 10ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . KNDy system remains unaffected (Figure 5(g)). Suppression of glutamatergic functional interaction strength combined with an overall increase in excitation of the MePD network causes a switch from inhibitory to excitatory mean MePD output (due to the presence of GABA-GABA interactions). This results in over-stimulation of the KNDy network associated with a transition from a pulsatile to a quiescent mode of operation due to a depolarisation-block-type of phenomenon (Figure 5(h)). (e) How stimulation of MePD projections modulates GnRH pulse generator activity The suppressive effect of MePD GABA projections’ output on the GnRH pulse generator frequency has been recently confirmed experimentally by optogenetically stimulating MePD GABAergic terminals in the ARC [16]. Here, using our coupled model, we interrogate the role of both GABAergic and glutamatergic projections of the MePD circuit in modulating KNDy dynamics. In the previous section the weights of GABA and glutamate projection output were set to be equal (j l = je = 1). Here in order to account for direct optogenetic simulation of these projections as carried out in [16] we increase the weight of the respective projection output in the model (i.e. for simulating the experiments presented in [16] we increase the weight of GABA projections output in the model) and vary the stimulation level in the MePD circuit ( Kp). Specifically, if we set je = 1.5 and jl = 0.5 we are able to show that increasing the stimulation level (K p) in the coupled system initially produces no change in the KNDy IPI, followed by an exponential increase in the KNDy inter pulse interval corresponding to cessation of the LH pulsatile dynamics (Figure 6( a)). These results conform with in vivo stimulation of MePD GABA projections, where a decrease in LH pulsatility (demonstrated by increase in inter pulse interval) eventually leading to loss of pulsatile dynamics have been observed at 10 and 20 Hz stimulation of GABA projections, respectively [16]. Now, we employ the same strategy to simulating the effects of glutamate projections’ stimulation (by fixing jl = 1.5 and je = 0.5). Our model simulations show that under increasing excitation of the MePD network, the inter pulse interval of KNDy population activity initially decreases before increasing and then returning to its initial IPI (Figure 6(a)). Furthermore, bifurcation analysis of the extended model demonstrates that the system undergoes a torus bifurcation (TR), associated with the switch from limit cycle dynamics to dynamics evolving on a limit torus, which we depict in (Figure 5(a)). The observed non-monotonic behaviour of the coupled system in this case is counterintuitive, given the excitatory role of glutamate. Nevertheless, selective in-vivo optogenetic stimulation of the MePD glutamatergic projections in the ARC with increasing levels of stimulation confirms our model predictions as shown in (Figure 6(b)). As expected, given the excitatory nature of the glutamatergic projections, sustained stimulation at 5 Hz results in a significant decrease in LH interpulse interval from 15.42 ± 0.60 min to 11.81 ± 0.81 min. However, further increase in the frequency of stimulation (at 10 Hz and 20 Hz) restores the pre-stimulation IPI levels, as predicted by our modelling. Moreover, our modelling allows us to explore potential mechanisms that govern the non-monotonic response in the GnRH pulse generator found experimentally as described above (Figure 6(c-d)). Continuation analysis of the MePD model dynamics indicates that the decrease in the KNDy network IPI observed in the model is due to the amplification of glutamatergic activity in the MePD network, while GABAergic tone remains very close to zero (see Figure 6(c) at stimulation level ≈ 1.2). Further increase in MePD network excitation, however, switches the balance in excitatory/inhibitory MePD projections’ output (i.e. MePD input to KNDy)(Figure 6(d), which in turn promotes an increase in the KNDy IPI. The MePD network model undergoes a Hopf bifurcation (HB), the location of which is associated with the location of the torus bifurcation (TR) in the extended MePD-KNDy network model (Figure 5( a-b)), demonstrating that the transition from limit cycle to limit torus occurs due to a change in the qualitative dynamics of the MePD network model. 4. Discussion In our study, we have introduced and systematically investigated a model incorporating the interplay between GABA and glutamate neuronal populations within the MePD. This model was coupled to a GnRH pulse generator model [6, 7], allowing us to validate it against experimental findings from [12, 13] as well as offering insights into how perturbations in the MePD could impact the activity of the GnRH pulse generator. Our model could serve as a versatile tool for investigating broader MePD circuit effects, such as for example those stemming from the interactions between urocortin and the GABA/glutamate neuronal populations [40]. The utility of our modelling approach lies in its ability to interrogate the neuronal mechanisms that enable the MePD to modulate the dynamics of the GnRH pulse generator, hence enabling us to better understand the effects of environmental and psychosocial factors on reproductive function, such as pubertal timing [9–11] and modulation of LH secretion [12, 13]. A model of the MePD circuit coupled to the GnRH pulse generator has been previously studied under stationary MePD circuit dynamics [13]. In this mode, the system functions like a collection of independent neuronal oscillators characterised by a constant (averaged) population level of activity. However, it is important to note that while this study has offered valuable insights, the actual patterns of MePD activity are likely to be more complex. Indeed, [31] shows changes in the MePD neuronal network’s oscillatory activity in vivo associated with sex- specific differences in the encoding of social stimuli and sexual experience. Here, we demonstrate that the extended model is able to reproduce experimental findings [12, 13], .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 11ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . suggesting the plausibility of an oscillatory mode of MePD circuit activity. In fact, rhythmicity of neuronal populations is a characteristic feature of neuronal synchronisation, allowing the neuronal networks to manage and process complex stimuli [41]. However, to confirm or reject the hypothesis about the importance of oscillations in the MePD neuronal networks, further experiments involving recordings of calcium activity in individual GABA and/or glutamate neurons in the MePD, and how they synchronise, would be required. Under oscillatory MePD circuit behaviour, our modelling shows that an increase in the excitatory input to the MePD system decreases GABAergic MePD output due to the activation of GABA-GABA disinhibition, while glutamatergic output remains consistent. This finding indicates that GABAergic MePD output is sensitive to stimulatory inputs, while glutamatergic output is likely to play more of a balancing role. This is consistent with the established role of amygdala in reproductive function modulation, as lesions to the MePD have been shown to advance puberty [42] and prevent stress-induced suppression of LH pulses [43]. On the other hand, optogenetic stimulation of kisspeptin neurons in the MePD [12] increases LH pulse frequency, and administration of peripheral kisspeptin inhibits neuronal activation in the amygdala as well as increases LH secretion [44], which can be explained by a decrease in the activity of GABAergic efferent neurons. Here, we have modelled the effects of pharmacological interventions via partial suppression of MePD circuit functional interactions and studied how different levels of suppression affect qualitative dynamics in the model; this is in contrast to our previous work [13] where we assumed complete suppression of signalling. In reality, however, complete suppression is unlikely to be the case, as neuronal cells may respond to pharmacological interventions by upregulating receptors or modifying their signalling pathways to compensate for the inhibited receptors. When modelling partial GABA signalling suppression, the inhibitory component of MePD output increases because there is not enough GABA-GABA disinhibition, but at the same time glutamatergic activity also goes up due to decreased inhibition from the GABAergic population, balancing out the inhibitory output. This compensation mechanism could provide an alternative explanation as to why solely GABA receptor antagonism does not change LH pulsatility [13]. Feeding oscillatory MePD projections’ output into the KNDy network allowed us to consider dynamic upstream modulation of the KNDy network rather than a constant input as considered in our previous work [7, 13]. We investigated how such dynamic input changes the response of the KNDy relaxation oscillator. Specifically, we have demonstrated that when the KNDy relaxation oscillator receives a dynamic input of a significantly different frequency, this can result in a complex quasi-periodic pulse pattern, i.e. irregular-shaped pulses with no change in interpulse interval. Quasi-periodicity has also been identified in other relaxation oscillator systems subject to periodic inputs, indicating that systems characterised by quasi-periodic behaviour often possess a level of resilience against external perturbations [38, 39]. On the other hand, the transition from oscillatory to constant input, which we observe during the suppression of glutamatergic functional interactions, makes the KNDy system far more sensitive to the magnitude of the change. Previous modelling work suggests high sensitivity of the KNDy network to the magnitude of constant external stimuli that can lead to cessation of GnRH pulses [7]. However, to maintain a functional reproductive system, the GnRH pulse generator must be resilient to small perturbations that can arise from changes in the MePD (or other upstream brain regions) output due to environmental stimuli, which is more plausible under periodic input modulation. Previously published data identified that stimulation of MePD GABA projections in the ARC modulates GnRH pulse generator activity [16]. Using our extended mathematical model we interrogated the role of the glutamatergic projections in such modulation, which has not been shown previously. The model analysis predicted a possible non-monotonic response of KNDy activity. These model predictions were experimentally confirmed in vivo, suggesting a novel excitatory MePD glutamatergic projection capable of impacting the KNDy network, and hence modulating the GnRH pulse generator. Specifically, we show that stimulation at 5 Hz leads to a statistically significant decrease in LH interpulse interval, but further increase in stimulation frequency (10 Hz and 20 Hz) produces no significant change. Considering the excitatory function of glutamate, one might question why an increase in stimulation does not lead to a further decrease in LH IPI or even the complete loss of pulsatility due to potential overstimulation (or in other words depolarisation block). Using our extended mathematical model we reproduced this non-monotonic stimulus- response relationship, providing insight into potential mechanisms though which glutamatergic projections modulate GnRH pulse generator activity. In the model, as the excitatory drive increases, so too does the response of excitatory projection neurons, leading to positive correlation between excitatory input and neuronal response, similar to the effect observed during 5 Hz stimulation of glutamatergic projections. However, when excitation becomes stronger, which mimics the effects of the higher frequency optogenetic stimulation, the system reaches a point where further increase in the excitatory drive does not lead to a proportional increase in the neuronal response. Excessive stimulation may also trigger compensatory mechanisms that counteract the increased activity, leading to a limited net change [45]. As we observe in the model, the inhibitory projections could become engaged to maintain the balance, counteracting the excitatory drive and contributing to a non-monotonic response. Based on the model analysis, we argue that balanced feed- forward excitation and feed-forward inhibition ensure that .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 12ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . the overall excitability of the KNDy network is robust to the elevated excitatory input, permitting GABAergic projections to exert ’inhibitory brake’, thus constraining the LH IPI to pre-stimulation levels. Here we propose a mathematical model and explore the influence of MePD activity on GnRH pulse generator dynamics. The modelling and analysis approach we have used could be applied to gain insight into the behaviour of other brain regions involved in modulation of the GnRH pulse generator. Additionally, given the phenomenological nature of our MePD model, in the future, the extended MePD-KNDy model could be used to interrogate the effects of other stimulatory neuronal populations signalling to the MePD GABA-glutamate circuit. This could then be employed to perform in silico simulations to interpret and/or predict experimentally observed effects on the GnRH pulse generator and reproductive function. Appendix Data Accessibility. The code to reproduce the analysis and data can be found in GitHub repository. Authors’ Contributions. Kateryna Nechyporenko, Conceptualisation, Investigation, Methodology, Software, Writing – original draft; Margaritis Voliotis, Conceptualisation, Investigation, Writing – review and editing; Xiao Feng Li, Conceptualisation, Data curation, Investigation; Owen Hollings, Conceptualisation, Data curation, Investigation; Deyana Ivanova, Conceptualisation, Investigation; Jamie Walker, Conceptualisation, Investigation, Writing - review and editing; Kevin O’Byrne, Conceptualisation, Investigation, Writing - review and editing; Krasimira Tsaneva-Atanasova, Conceptualisation, Investigation, Writing - review and editing. Funding. KTA gratefully acknowledges the financial support of the EPSRC via grant EP/T017856/1. KTA and MV gratefully acknowledge the financial support of the BBSRC via grant BB/W005883/1. KOB and XFLI gratefully acknowledge the financial support of the BBSRC via grant BB/W005913/1. JJW gratefully acknowledges the financial support of the MRC via grants MR/N008936/1 and MR/T032480/1. Disclaimer. For the purpose of open access, the author has applied a ‘Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

1P . E. Belchetz, T. M. Plant, Y. Nakai, E. J. Keogh, and E. Knobil, “Hypophysial responses to continuous and intermittent delivery of hypopthalamic gonadotropin-releasing hormone”, Science 202, Belchetz, P E Plant, T M Nakai, Y Keogh, E J Knobil, E Journal Article United States 1978/11/10 Science. 1978 Nov 10;202(4368):631-3. doi: 10.1126/science.100883., 631–3 (1978). 2A. E. Herbison, “Control of puberty onset and fertility by gonadotropin-releasing hormone neurons”, Nature Reviews Endocrinology 12, 452–466 (2016). 3K. T. O’Byrne and E. Knobil, “Electrophysiological approaches to gonadotrophin releasing hormone pulse generator activity in the rhesus monkey”, Hum Reprod 8 Suppl 2, 37–40 (1993). 4P . Grachev and R. L. Goodman, “The GnRH Pulse Generator”,AIMS Medical Science 3, 359–385 (2016). 5A. E. Herbison, “The Gonadotropin-Releasing Hormone Pulse Generator”, Endocrinology 159, 3723–3736 (2018). 6M. Voliotis, X. F. Li, R. De Burgh, G. Lass, S. L. Lightman, K. T. O’Byrne, and K. Tsaneva-Atanasova, “The Origin of GnRH Pulse Generation: An Integrative Mathematical-Experimental Approach”, The Journal of Neuroscience 39, 9738–9747 (2019). 7M. Voliotis, X. F. Li, R. A. De Burgh, G. Lass, D. Ivanova, C. McIntyre, K. O’Byrne, and K. Tsaneva-Atanasova, “Modulation of pulsatile GnRH dynamics across the ovarian cycle via changes in the network excitability and basal activity of the arcuate kisspeptin network”, eLife 10, 10.7554/elife.71252 (2021). 8J. Úbeda-Contreras, I. Marín-Blasco, R. Nadal, and A. Armario, “Brain c-fos expression patterns induced by emotional stressors differing in nature and intensity”, Brain Structure and Function 223, 2213–2227 (2018). 9X. F. Li, D. A. Adekunbi, H. M. Alobaid, S. Li, M. Pilot, S. L. Lightman, and K. T. O’Byrne, “Role of the posterodorsal medial amygdala in predator odour stress-induced puberty delay in female rats”, Journal of Neuroendocrinology 31, e12719 (2019). 10D. Ivanova, X. Li, Y. Liu, C. McIntyre, C. Fernandes, G. Lass, L. Kong, and K. T. O’Byrne, “Role of Posterodorsal Medial Amygdala Urocortin-3 in Pubertal Timing in Female Mice”, Front Endocrinol (Lausanne) 13, 893029 (2022). 11D. A. Adekunbi, X. F. Li, S. Li, O. A. Adegoke, B. O. Iranloye, A. O. Morakinyo, S. L. Lightman, P . D. Taylor, L. Poston, and K. T. O’Byrne, “Role of amygdala kisspeptin in pubertal timing in female rats”,PLOS ONE 12, e0183596 (2017). 12G. Lass et al., “Optogenetic stimulation of kisspeptin neurones within the posterodorsal medial amygdala increases luteinising hormone pulse frequency in female mice”, Journal of Neuroendocrinology 32, 10.1111/jne.12823 (2020). 13G. Lass et al., “GnRH pulse generator frequency is modulated by kisspeptin and GABA-glutamate interactions in the posterodorsal medial amygdala in female mice”, Journal of Neuroendocrinology 34, 10.1111/jne.13207 (2022). 14A. M. Moore, L. M. Coolen, and M. N. Lehman, “Kisspeptin/Neurokinin B/Dynorphin (KNDy) cells as integrators of diverse internal and external cues: evidence from viral-based monosynaptic tract-tracing in mice”, Scientific Reports 9, 10.1038/s41598-019-51201-0 (2019). .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 13ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15S.-H. Yeo, V . Kyle, C. Blouet, S. Jones, and W. H. Colledge, “Mapping neuronal inputs to Kiss1 neurons in the arcuate nucleus of the mouse”, PLOS ONE 14, e0213927 (2019). 16C. McIntyre, X. F. Li, R. de Burgh, D. Ivanova, G. Lass, and K. T. O’Byrne, “GABA Signaling in the Posterodorsal Medial Amygdala Mediates Stress-induced Suppression of LH Pulsatility in Female Mice”, Endocrinology 164, 10 . 1210 / endocr/bqac197 (2022). 17K. B. J. Franklin and G. Paxinos, The Mouse Brain in Stereotaxic Coordinates (Boston, Amsterdam, 2004). 18F. J. Steyn, Y. Wan, J. Clarkson, J. D. Veldhuis, A. E. Herbison, and C. Chen, “Development of a Methodology for and Assessment of Pulsatile Luteinizing Hormone Secretion in Juvenile and Adult Male Mice”,Endocrinology 154, 4939–4945 (2013). 19R. B. McCosh, M. J. Kreisman, and K. M. Breen, “Frequent Tail-tip Blood Sampling in Mice for the Assessment of Pulsatile Luteinizing Hormone Secretion”, Journal of Visualized Experiments, 10.3791/57894 (2018). 20A. Vidal, Q. Zhang, C. Médigue, S. Fabre, and F. Clément, “DynPeak: An Algorithm for Pulse Detection and Frequency Analysis in Hormonal Time Series”, PLOS ONE 7, e39001 (2012). 21J. A. Yang, J. K. Hughes, R. A. Parra, K. M. Volk, and A. S. Kauffman, “Stress rapidly suppresses in vivo LH pulses and increases activation of RFRP-3 neurons in male mice”, Journal of Endocrinology 239, 339–350 (2018). 22G. B. Choi, H.-W. Dong, A. J. Murphy, D. M. Valenzuela, G. D. Yancopoulos, L. W. Swanson, and D. J. Anderson, “Lhx6 Delineates a Pathway Mediating Innate Reproductive Behaviors from the Amygdala to the Hypothalamus”, Neuron 46, 647–660 (2005). 23C. S. Johnson, W. Hong, and P . E. Micevych, “Posterodorsal Medial Amygdala Regulation of Female Social Behavior: GABA versus Glutamate Projections”, The Journal of Neuroscience 41, 8790–8800 (2021). 24H. R. Wilson and J. D. Cowan, “Excitatory and inhibitory interactions in localized populations of model neurons”, Biophys J 12, 1–24 (1972). 25H. R. Wilson and J. D. Cowan, “A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue”, Kybernetik 13, 55–80 (1973). 26S. Keshavarzi, R. K. P . Sullivan, D. J. Ianno, and P . Sah, “Functional Properties and Projections of Neurons in the Medial Amygdala”, The Journal of Neuroscience 34, 8699–8715 (2014). 27X. Bian, “Physiological and morphological characterization of GABAergic neurons in the medial amygdala”, Brain Res 1509, 8–19 (2013). 28K. J. Friston, Statistical Parametric Mapping : the Analysis of Functional Brain Images (Elsevier / Academic Press, Amsterdam ; Boston, 2007). 29H. R. Wilson and J. D. Cowan, “Evolution of the Wilson–Cowan equations”, Biological Cybernetics 115, 643–653 (2021). 30R. M. Borisyuk and A. B. Kirillov, “Bifurcation analysis of a neural network model”, Biological Cybernetics 66, 319–325 (1992). 31Y. Li, A. Mathis, B. F. Grewe, J. A. Osterhout, B. Ahanonu, M. J. Schnitzer, V . N. Murthy, and C. Dulac, “Neuronal Representation of Social Information in the Medial Amygdala of Awake Behaving Mice”, Cell 171, 1176–1190.e17 (2017). 32J. Qiu, C. C. Nestor, C. Zhang, S. L. Padilla, R. D. Palmiter, M. J. Kelly, and O. K. Rønnekleiv, “High-frequency stimulation- induced peptide release synchronizes arcuate kisspeptin neurons and excites GnRH neurons”,Elife 5, 10.7554/eLife. 16246 (2016). 33K. Yasuda, K. Raynor, H. Kong, C. D. Breder, J. Takeda, T. Reisine, and G. I. Bell, “Cloning and functional comparison of kappa and delta opioid receptors from mouse brain.”, Proceedings of the National Academy of Sciences 90, 6736–6740 (1993). 34G. R. Seabrook, B. J. Bowery, and R. G. Hill, “Pharmacology of tachykinin receptors on neurones in the ventral tegmental area of rat brain slices”, Eur J Pharmacol 273, 113–9 (1995). 35A. Dutton and R. E. Dyball, “Phasic firing enhances vasopressin release from the rat neurohypophysis”, J Physiol 290, 433–40 (1979). 36E. Doedel, T. Fairgrieve, B. Sandstede, A. Champneys, Y. Kuznetsov, and X. Wang,AUTO-07P: Continuation and bifurcation software for ordinary differential equations, 2007. 37J. Kim, S. J. Semaan, D. K. Clifton, R. A. Steiner, S. Dhamija, and A. S. Kauffman, “Regulation of Kiss1 expression by sex steroids in the amygdala of the rat and mouse”, Endocrinology 152, 2020–30 (2011). 38K. Kohari, T. Saito, and H. Kawakami, “On a relaxation oscillator with periodic threshold”, in 1993 IEEE International Symposium on Circuits and Systems (ISCAS) (1993), 2490–2493 vol.4. 39B. Christiansen, P . Alstrm, and M. T. Levinsen, “Routes to chaos and complete phase locking in modulated relaxation oscillators”, Phys. Rev. A 42, 1891–1900 (1990). 40D. Ivanova, X. F. Li, C. McIntyre, and K. T. O’Byrne, “Posterodorsal Medial Amygdala Urocortin-3, GABA, and Glutamate Mediate Suppression of LH Pulsatility in Female Mice”, Endocrinology 164, 10.1210/endocr/bqac196 (2022). 41W. Singer, “Neuronal synchrony: a versatile code for the definition of relations?”,Neuron 24, 49–65, 111–25 (1999). 42X. F. Li, M. H. Hu, B. P . Hanley, Y. S. Lin, L. Poston, S. L. Lightman, and K. T. O’Byrne, “The Posterodorsal Medial Amygdala Regulates the Timing of Puberty Onset in Female Rats”, Endocrinology 156, 3725–3736 (2015). .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 14ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43Y. Lin, X. Li, M. Lupi, J. S. Kinsey-Jones, B. Shao, S. L. Lightman, and K. T. O’Byrne, “The Role of the Medial and Central Amygdala in Stress-Induced Suppression of Pulsatile LH Secretion in Female Rats”, Endocrinology 152, 545–555 (2011). 44A. N. Comninos et al., “Kisspeptin signaling in the amygdala modulates reproductive hormone secretion”, Brain Structure and Function 221, 2035–2047 (2016). 45S. B. Wolff and B. P . Ölveczky, “The promise and perils of causal circuit manipulations”,Curr Opin Neurobiol 49, 84–94 (2018). .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 15ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 1. (a) Schematic diagram of the GABA-glutamate circuit in the MePD. Kisspeptin activates glutamatergic and GABAergic neuronal populations in the MePD. Both glutamatergic and GABAergic populations interact functionally with a population of GABA efferent neurons providing excitation and inhibition, respectively. The activity within the GABAergic efferent neurons and glutamatergic neurons extends to the ARC to modulate the activity of the KNDy population. (b) Network representation of the MePD circuit model depicting interactions between the populations of glutamatergic neurons( Gl), GABA interneurons (Gi) and GABAergic efferents (Ge). Both Gl and Gi receive kisspeptin excitation (Kp), that is distributed proportionally with relative glutamatergic excitation ratio α. The parameters c represent the strength of functional interactions between the populations. Parameters β1 and β2 are the interaction suppression coefficients, used to mimic the decrease in the GABAergic and glutamatergic interaction strengths, respectively, in order to simulate the effects of chemical antagonism. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 16ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 2 6 10 14 MePD excitation 0 0.1 0.2 0.3 0.4 0.5Glu population activity a HB 0 2 6 10 14 MePD excitation 0 0.1 0.2 0.3 0.4 0.5 GABAint population activity b HB 0 2 6 10 14 MePD excitation 0 0.1 0.2 0.3 0.4 0.5 GABAeff population activity c HB 0 100 200 300 Time (s) 0 0.1 0.2 0.3 0.4Averaged Activity d 2 4 6 8 10 12 14 MePD excitation 20 40 60 80 100 120 Period (s) f 0 100 200 300 Time (s) 0 0.1 0.2 0.3 0.4Averaged Activity e Glu GABAeff 0 2 4 6 8 10 12 14 16 18 20 MePD excitation 0.7 0.8 0.9 1 Relative Glu excitation ratio g BT HB SN SN HC -0.2 -0.1 0 0.1 0.2 Mean MePD output Figure 2. MePD GABA-glutamate circuit dynamics. (a-c) One-parameter bifurcation diagrams of the model for varying kisspeptin (K p) excitation of the GABA-glutamate circuit (MePD excitation). The qualitative dynamics in the populations of (a) glutamatergic neurons (Glu), (b) GABA interneurons (GABAint) and (c) GABAergic efferent neurons (GABAeff) changes as the kisspeptin excitation level provided to the circuit varies. Circular markers denote a Hopf (HB) bifurcation that gives rise to the oscillatory behaviour in the model. The red line shows the maximum amplitude of the limit cycle solution branch. (d-e) Simulation of population activity of Glu (green) and GABA eff (red) at two different levels of MePD excitation (K p = 2.3 and Kp = 6), corresponding to the yellow and blue triangular markers in panels a-c. (f) MePD excitation vs period of the system’s limit cycle oscillations. ( g) Two- parameter bifurcation diagram using the level of MePD excitation (K p) and the relative Glu excitation ratio (α) as free parameters. Superimposed heat map indicates the mean MePD output in the oscillatory region. Red, green, and purple lines represent Hopf (HB), saddle-node (SN) and homoclinic bifurcations (HC), respectively. Circular markers depict a Bogdanov-Takens (BT) point. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 17ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 3. Effect of functional interaction strength between the population of GABA interneurons (GABAint) and the population of GABA efferent neurons GABAeff on the MePD projections’ output. Heat map depicting (a) the magnitude of mean glutamatergic projections’ output, (b) the magnitude of the mean GABAergic output under varying MePD excitation (Kp) and strength of the interaction between the population of GABA interneurons and the population of GABA efferent neurons (c ie). For low strength of interaction, GABAergic output increases as MePD excitation is increased, while excessively high interaction drives the activity of the population of GABA efferent (projection) neurons to a zero level. (c) Heat map depicting the combined mean output of glutamate and GABA projections. The dashed line represents the selected functional interaction strengths used in the model (see Table 1). .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 18ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . a b Oscillatory regime Oscillatory regime Figure 4. Effect of varying the MePD functional network connectivity on the MePD projections’ output dynamics. (a) Bifurcation analysis under varying functional interaction strength between the populations of GABAergic interneurons and glutamatergic neurons (cil and cli). (b) Bifurcation analysis under varying strength of glutamatergic self-excitation (c ll) and functional interaction strength between glutamatergic neurons and GABAergic interneurons (cli). The network schematic diagrams highlight the connections under investigation. In the one-parameter bifurcation diagram, the red lines indicate the max of the limit cycle solutions branch, and circle markers denote Hopf ( HB), saddle node (SN), and homoclinic (HC) bifurcation points, respectively. These points are used for the two-parameter continuation. In the two-parameter bifurcation diagrams, red, green, and purple lines represent Hopf (HB), saddle-node (SN) and homoclinic bifurcation curves (HC), respectively. The circular markers depict the Bogdanov-Takens (BT), saddle node loop (SNL) and cusp (CP) points. The heat map represents the mean MePD projections’ output in the oscillatory region. The star marker identifies the parameter values in Table 1. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 19ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . b edc hgf Figure 5. Differential effects of MePD projections on the dynamics of the GnRH pulse generator. (a) Providing oscillatory input to the KNDy network model leads to quasi-periodic behaviour evolving on a torus; we use a colour scheme to illustrate time evolution on this torus. ( b) Increasing excitation in the MePD leads to a decrease in GABAergic tone, and as a result a decrease in interpulse interval (IPI) in the KNDy system. ( c) Two-parameter bifurcation diagram (MePD excitation (K p) and GABAergic interaction suppression coefficient (β 1)) showing the region where the MePD exhibits oscillatory dynamics and the mean MePD output superimposed on it. (d) Suppression of GABAergic interaction evenly increases both GABAergic and glutamatergic tone, leaving the overall MePD output as before the suppression, hence KNDy IPI remains the same. (e) Conjunction of the suppression of GABAergic interaction strength with increased MePD stimulation increases the proportion of GABAergic projections’ output, resulting in an increase in KNDy IPI. (f) Two-parameter bifurcation diagram (MePD excitation (Kp) and glutamatergic interaction suppression coefficient (β 2)) showing the region where the MePD exhibits oscillatory dynamics with superimposed mean MePD projections’ output. (g) Suppression of glutamatergic interaction strength leads to a loss of oscillatory dynamics in the MePD. Inhibitory GABAergic output is still higher than the excitatory glutamatergic output, leading to the return of periodicity in the KNDy system and leaving KNDy IPI as before the suppression. (h) Combining the suppression of glutamatergic interaction strength with the increased excitation of the MePD overstimulates ARC KNDy network, resulting in a transition to quiescent dynamics. A more detailed versions of the two-parameter bifurcation diagrams depicted in (c) and (f) are provided in Figure S2. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 20ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Stimulation 10 15 20 25Average KNDy IPI (min) a TR Glu projection GABA projection 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Stimulation 0 0.1 0.2 0.3Weighted MePD activity(Hz) c HB Glu GABAeff Ctrl Ctrl Virus 2 Hz 5 Hz 10 Hz 20 Hz 10 15 20LH IPI (min) b ## Prestimulation Stimulation No stim ctrl 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Stimulation 0 0.02 0.04 0.06Mean MePD output (Hz) d Pulses cease Figure 6. In silico simulation of MePD projection stimulation and non-monotonic LH interpulse interval response to stimulation of MePD glutamatergic projection in vivo. (a) Simulation of GnRH pulse generator response to stimulation of GABA and glutamate MePD projections. Mimicking the effects of GABA projection stimulation leads to the loss of oscillatory dynamics in the KNDy activity as stimulation level is increased. During glutamate projection stimulation in silico KNDy IPI initially declines and then rises back as the level of stimulation (K p) (given in arbitrary units (a.u.)) increases. A torus bifurcation (TR) denotes the transition from limit cycle to limit torus in the extended MePD-KNDy network model. The location of TR and HB align. (a) Summary of average LH interpulse interval (IPI) for non-stimulation control (n = 5), control virus (n = 4) and optical stimulation at 2 Hz (n = 5), 5 Hz (n = 6), 10 Hz (n = 5), and 20 Hz (n = 5) of MePD glutamatergic projection terminals in the ARC. We depict the mean ± SEM, and individual data points for each animal as circles in the histogram plots. Stimulation at 5 Hz results in a statistically significant suppression of LH pulse frequency ##F = 12.8517, p = 0.0050 vs control period for 5 Hz optical stimulation group. ( c) The external MePD projections’ input to the KNDy network, showing glutamatergic and GABAergic contributions as a function of the level of stimulation (in a.u.). Hopf point (HB) point denotes the location of the transition from stationary to oscillatory dynamics in the MePD circuit. (d) The variation in the weighted mean MePD output under different levels of stimulation. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 21ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 2 6 10 14 18 22 MePD excitation 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Glu population activity a HB HC SN SN 0 2 6 10 14 18 22 MePD excitation 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 GABAeff population activity b HB HC SN SN 0 5 10 15 20 MePD excitation 5 10 15 log(Period (s)) c Figure S1. Extended bifurcation diagrams of kisspeptin-mediated excitatiry input to GABA-glutamate circuit fromFigure 2 for (a) Gl and (b) Ge dynamics. Circle markers denote Hopf (HB), homoclinic (HC) and saddle node (SN) bifurcations. The red line indicates the maximum amplitude of the limit cycle solutions branch. (c) Excitatory input vs. log Period of oscillations. The period, firstly, decreases, followed by a rapid increase, signifying homoclinic bifurcation. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint 22ro yalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 2 4 6 MePD excitation 0 0.2 0.4 0.6 0.8 1 GABA suppression a CP CP BT SNL SNL SNL SN HB HC SNIC 0 1 2 3 4 5 6 7 MePD excitation 0 0.2 0.4 0.6 0.8 1 Glutamate suppression b HB Figure S2. Two parameter bifurcation diagrams, used to identify the oscillatory regions in Figure 5(c,f). (a) Two-parameter bifurcation diagram of MePD excitation (K p) and GABAergic interaction suppression coefficient (β 1). (b) Two-parameter bifurcation diagram of MePD excitation (K p) and glutamtergic interaction suppression coefficient (β2). Red, green, purple lines represent Hopf (HB), saddle-node (SN) and homoclinic bifurcation curves (HC), respectively. The circular markers depict Bogdanov-Takens (BT), saddle-node loop (SNL) and cusp (CP) points. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint