{"paper_id":"0ea98ea4-f481-461b-a45e-d54dc48277e1","body_text":"rsif.royalsocietypublishing.org\nResearch\nArticle submitted to journal\nKeywords:\nneuromodulation, amygdala,\nmathematical model\nAuthor for correspondence:\nKateryna Nechyporenko\ne-mail: kn356@exeter.ac.uk\nKrasimira Tsaneva-Atanasova\ne-mail:\nK.Tsaneva-Atanasova@exeter.ac.uk\nNeuronal Network Dynamics in the\nPosterodorsal Amygdala: Shaping\nReproductive Hormone Pulsatility\nKateryna Nechyporenko 1,2, Margaritis Voliotis1,2, Xiao Feng\nLi3, Owen Hollings3, Deyana Ivanova4, Jamie J Walker1,\nKevin T O’Byrne3, Krasimira Tsaneva-Atanasova1,2,5\n1 Department of Mathematics and Statistics, University of Exeter, Stocker Rd,\nExeter, EX4 4PY , United Kingdom\n2 Living Systems Institute, University of Exeter, Stocker Rd, Exeter, Ex4 4PY ,\nUnited Kingdom\n3 Department of Women and Children’s Health, School of Life Course and\nPopulation Sciences,King’s College London, Guy’s Campus, London, SE1\n1UL, United Kingdom\n4 Division of Endocrinology, Diabetes, and Hypertension, Brigham and\nWomen’s Hospital, Harvard Medical School, Boston, MA, USA\n5 EPSRC Hub for Quantitative Modelling in Healthcare, University of Exeter,\nStocker Rd,Exeter, Ex4 4PY , United Kingdom\nNormal reproductive function and fertility rely on the rhythmic secretion of\ngonadotropin-releasing hormone (GnRH), which is driven by the hypothalamic\nGnRH pulse generator. A key regulator of the GnRH pulse generator is the\nposterodorsal subnucleus of the medial amygdala (MePD), a brain region\nwhich is involved in processing external environmental cues, including the\neffect of stress. However, the neuronal pathways enabling the dynamic,\nstress-triggered modulation of GnRH secretion remain largely unknown.\nHere, we employ in-silico modelling in order to explore the impact of\ndynamic inputs on GnRH pulse generator activity. We introduce and analyse\na mathematical model representing MePD neuronal circuits composed of\nGABAergic and glutamatergic neuronal populations, integrating it with our\nGnRH pulse generator model. Our analysis dissects the influence of excitatory\nand inhibitory MePD projections’ outputs on the GnRH pulse generator’s\nactivity and reveals a functionally relevant MePD glutamatergic projection to\nthe GnRH pulse generator, which we probe with in vivo optogenetics. Our\nstudy sheds light on how MePD neuronal dynamics affect the GnRH pulse\ngenerator activity, and offers insights into stress-related dysregulation.\n© The Authors. Published by the Royal Society under the terms of the\nCreative Commons Attribution License http://creativecommons.org/licenses/\nby/4.0/, which permits unrestricted use, provided the original author and\nsource are credited.\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n2ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\n1. Introduction\nThe rhythmic secretion of gonadotropin-releasing hormone\n(GnRH) from the hypothalamus into the portal circulation\nis crucial in triggering the pulsatile release of gonadotropin\nhormones (luteinizing hormone (LH) and follicle-\nstimulating hormone (FSH)) from the pituitary gland [1,\n2]. This dynamic process significantly contributes to the\ninitiation of puberty and plays a pivotal role in ensuring\nfertility [2, 3]. The pulsatile release of GnRH is controlled\nby an upstream brain network known as the \"GnRH pulse\ngenerator\". This network is composed of KNDy neurons\nfound in the arcuate nucleus (ARC) of the hypothalamus,\nwhich co-express kisspeptin, neurokinin B and dynorphin\nA [4, 5]. We have previously used a mathematical model\nof the KNDy network to show that periodic pulsatile\nactivity emerges as the basal activity or external activation\nof the network is increased, and confirmed this modelling\nprediction in vivo using optogenetic stimulation of the\nKNDy network [6, 7].\nThe GnRH pulse generator is sensitive to inputs from\nneural centres relaying information about the psychological\nand physiological state of the organism. In particular, the\nposterodorsal subnuclei of the medial amygdala (MePD),\na limbic structure responsible for emotional processing\nof complex external cues, is responsive to stress [8] and\nregulates pubertal onset [9–11] as well as LH pulsatility\n[12, 13]. Whilst there is limited understanding of how\nthe MePD processes stress-related information and relays\nit to the GnRH pulse generator, it has been shown that\noptogenetic stimulation of kisspeptin neurons in the MePD\nincreases LH pulse frequency [12]. This effect is mediated\nby both GABAergic and glutamatergic signalling within\nthe MePD, since pharmacological antagonism of GABA\nreceptors within the MePD prevents the increase in LH\npulse frequency, while pharmacological antagonism of\nglutamate receptors within the MePD terminates the LH\npulses altogether [13]. It is likely that MePD regulation\nof LH pulsatility is mediated, at least in part, through\ndirect projections from the MePD to KNDy neurons within\nthe ARC. Indeed, viral-based monosynaptic tract-tracing\nin mice has shown that the amygdala provides inputs\nto ARC KNDy neurons [14, 15]. Furthermore, although\nglutamatergic projections from the MePD have not yet been\nexplored experimentally, it has been shown that stimulation\nof MePD GABAergic projection terminals in the ARC causes\na suppression of LH pulses [16].\nTaken together, the available data suggest that the MePD\ncan modulate GnRH pulse generator activity through local\nkisspeptin signalling. However, how the GABA-glutamate\nneuronal network within the MePD integrates kisspeptin\nactivity and relays this to the GnRH pulse generator is\nless clear. Here, we investigate this using a mathematical\nmean-field model of the MePD GABA-glutamate neuronal\ncircuit (see Figure 1(a)). We calibrate the model using\navailable experimental observations and analyse it to\nunderstand how changes in functional connectivity within\nthe MePD’s neuronal network affect the system’s output.\nWe then couple this model of the MePD with our earlier\nmathematical model of the GnRH pulse generator [6, 7]\nto study the effects of manipulating both GABAergic and\nglutamatergic signalling on the activity of the GnRH pulse\ngenerator.\n2. Methods\n(a) Animals\nAdult Vglut-flp mice heterozygous for the allele\nSlc17a6<tm1.1(flpo)Hze> (Strain #:030212, B6;129S-\nSlc17a6tm1.1(flpo)Hze/J; Jackson Laboratory, Bar Harbor,\nME, USA) were bred in house. Mice were genotyped by PCR\nusing the following primers: mutant reverse, 13007- ACA\nCCG GCC TTA TTC CAA G; common, 34763- GAA ACG\nGGG GAC ATC ACT C; and wildtype reverse, 34764- GGA\nATC TCA TGG TCT GTT TTG. Mice, aged7-9 weeks at time\nof initial surgery, were group housed unless chronically\nimplanted with fiber optic cannulae and kept at25°C ± 1°C,\n12:12 hr light/dark cycle (lights on 0700 h), with ad libitum\naccess to food and water. All procedures were performed in\naccordance with UK home office regulation and approved\nby The King’s College London Animal Welfare and Ethical\nReview Body.\n(b) Stereotaxic injection of viral constructs\nand implantation of fibre optic cannula\nAll surgical procedures were carried out under general\nanaesthesia using ketamine (Vetalar,100 mg/kg, i.p.; Pfizer,\nSandwich, UK) and xylazine (Rompun, 10 mg/kg, i.p.;\nBayer, Leverkusen, Germany) under aseptic conditions.\nMice were secured in a David Kopf stereotaxic frame (Model\n900, Kopf Instruments) and bilaterally ovariectomised\n(OVX). A midline incision of the scalp was used to expose\nthe skull. The periosteum was removed and two small\nbone screws were inserted into the skull. Using a robot\nstereotaxic system (Neurostar, Tubingen, Germany) two\nwindows were drilled intracranially directly above target\ncoordinates for the posterodorsal subnucleus of the medial\namygdala (MePD) (2.35 mm lateral, 1.45 mm posterior to\nbregma, at a depth of 5.49 mm below the skull surface) and\narcuate nucleus (ARC) (0.24 mm lateral, 1.51 mm posterior\nto bregma, at a depth of 5.85 mm below the skull surface)\nobtained from the mouse brain atlas of [17]. Either AAV-\nCAG-FLEXFRT-ChR2(H134R)-mCherry (n = 7, 200 nL, 3 ×\n1012 GC/mL, Serotype: 9; Addgene, MA, USA) to express\nchannelrhodopsin (ChR2) or control virus (n= 4, AAV-Ef1a-\nfDIO-mCherry, 200 mL, 3 × 1012 GC/mL; Serotype: 9;\nAddgene) was unilaterally injected over 10 minutes into the\nMePD using a 2-µL Hamilton microsyringe (Esslab, Essex,\nUK). The needle was left in place for 5 minutes, before\nwithdrawing 0.2 mm, waiting another minute, and then\nfully withdrawing over the course of 1 minute. A fiber\noptic cannula (200 µm,0.39 NA, 1.25 mm ceramic ferrule,\nDoric Lenses, Quebec, Canada) was inserted into the brain\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n3ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\ntargeting the ARC. This was fixed to the skull surface and\nbone screws using dental cement (Super-Bond Universal\nKit, Prestige Dental, UK) before suturing closed the incision.\nA one week recovery period was given post-surgery. After\nthis period, the mice were handles daily to acclimatize them\nto the tail-top blood sampling procedure [18]. Mice were\nleft for 4 weeks to achieve effective opsin expression before\nexperimentation.\n(c) Experimental design and blood sampling\nfor LH measurement\nFor measurement of LH pulsatility during optogenetic\nstimulation, the tip of the mouse’s tail was removed\nwith a sterile scalpel for tail-tip blood sampling [19]. The\nchronically implanted fiber-optic cannula was attached to\na multimode fiber-optic rotary joint patch cables (Thorlabs,\nLtd, Ely, UK) via a ceramic mating sleeve which allows\nmice to freely move while receiving blue light (473 nm\nwavelength). Laser (DPSS laser, Laserglow Technologies,\nToronto, Canada) intensity was set to10 mW at the tip of the\nfibre optic patch cable. The frequency and pattern of optical\nstimulation was controlled by software designed in house.\nAfter 1 h acclimatization, blood samples (5µl) were collected\nevery 5 min for 2 h. After 1 h controlled blood sampling\nwithout optical stimulation to determine baseline LH pulse\nfrequency, mice received patterned optical stimulation (5 s\non 5s off, 10-ms pulse width) at2, 5, 10 or 20 Hz for 1 h while\nblood sampling continued. Non-stimulation controls were\nperformed in the same manner, with no stimulation during\nthe second hour. The control virus injected animals received\n5 Hz optical stimulation. Mice received all treatments in a\nrandom order, with at least 3 but typically 5 days between\nexperiments.\nThe blood samples were snap-frozen on dry ice and\nstoring at −80°C until processed. In-house LH enzyme-\nlinked immunosorbent assay (LH ELISA) similar to\nthat described by [18] was used for processing of the\nmouse blood samples. The mouse LH standard (AFP-\n5306A; NIDDK-NHPP) was purchased from Harbor-UCLA\nalong with the primary antibody (polyclonal antibody,\nrabbit LH antiserum, AFP240580Rb; NIDDK-NHPP). The\nsecondary antibody (donkey anti-rabbit IgG polyclonal\nantibody [horseradish peroxidase]; NA934) was from VWR\nInternational. Validation of the LH ELISA was done in\naccordance with the procedure described in [18] derived\nfrom protocols defined by the International Union of Pure\nand Applied Chemistry. Serially diluted mLH standard\nreplicates were used to determine the linear detection range.\nNonlinear regression analysis was performed using serially\ndiluted mLH standards of known concentration to create\na standard curve for interpolating the LH concentration in\nwhole blood samples, as described previously [6]. The assay\nsensitivity was 0.031 ng/mL, with intra- and inter-assay\ncoefficients of variation of 4.6% and 10.2% respectively.\n(d) LH pulse detection and statistical analysis\nLH pulses were determined using the DynPeak algorithm\n[20], with settings adjusted to accommodate for the\nhigh LH pulse frequency in OVX mice as previously\noutlined by Breen and colleagues [21]. These include using\nthe programmes’ default parameters except the global\nthreshold was increased to35%, the nominal peak threshold\nwas reduced to 20 min and the 3-point peak threshold was\nremoved. Average LH inter-pulse interval (IPI) (the period\nof time between two LH pulse peaks) was calculated for the\n1 h control period and 1 h optogenetic stimulation period\nor equivalent non-stimulation period in control animals.\nStatistical significance was tested using a two-way repeated\nmeasures ANOVA and post-hoc Tukey test. Data was\nrepresented as mean ± SEM and p < 0.05 was considered\nsignificant.\n(e) Mean-Field Model of the MePD\nGiven the established presence of GABA and glutamate\nneuronal populations in the MePD [22, 23], we model\ntheir interplay employing Wilson-Cowan framework\n[24, 25]. The framework allows us to take a system-\nlevel approach to describe the dynamic evolution of\nexcitatory/inhibitory activity in neuronal populations\ndue to functional interactions within a synaptically-\ncoupled neuronal network, incorporating both cooperation\nand competition mechanisms. Rather than considering\nindividual neurons within the populations, the framework\ngives a coarse-grained account of the mean activity in\nthe network, which allows the investigation of putative\nfunctional interactions between the various populations as\nwell as the overall network output, enabling coupling of\nour MePD network model to other neuronal networks, such\nas the GnRH pulse generator as represented by the KNDy\nnetwork in the ARC[6, 7].\nThe Wilson-Cowan framework includes an inhibitory\nand excitatory dependant variables, both receiving an\nexcitatory input. However, straightforward application of\nthis framework is not sufficient to describe the MePD\nGABA-glutamate neuronal network to represent differential\ndynamics of GABA-mediated disinhibitory mechanism [26,\n27]. Therefore, we have extended the original Wilson-\nCowan model by incorporating an additional inhibitory\npopulation (GABA), that does not receive excitatory input\n(see Figure 1(b)).\nA key component of the Wilson-Cowan modelling\nframework is a sigmoid stimulus-response function\nϕ(a, F, θ), which controls the mean level of activity\ngenerated in the populations at a time t [24]:\nϕ(a, F, θ) = 1\n1 + exp(−a(F − θ)) − 1\n1 + exp(aθ) , (2.1)\nwhere F indicates the input to a given population and\nparameters a and θ define the value of maximum slope and\nhalf-maximum firing threshold, respectively. The use of the\nsigmoidal function is motivated by the fact that the majority\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n4ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\nof neurons have fluctuating membrane potential near an\nexcitability threshold, such that the probability of firing\ngrows exponentially upon depolarisation [28]. Additionally,\nthe constant 1/(1 + exp(aθ)) ensures that under the absence\nof stimulatory input to the population the firing ceases, i.e.\nϕ(a, 0, θ) = 0 [24].\nThe input F is given by the linear sum of excitatory and\ninhibitory contributions, as follows:\nFl(Gl, Gi, Ge) =(1 − β2)cllGl − (1 − β1)cilGi + αKp,\n(2.2)\nFi(Gl, Gi, Ge) =(1 − β2)cliGl + (1 − α)Kp, (2.3)\nFe(Gl, Gi, Ge) =(1 − β2)cleGl − (1 − β1)cieGi, (2.4)\nwhere the parameters c represent the strength of interaction\nfrom one population to another, as shown in Figure 1(b).\nThe parameter Kp depicts the overall kisspeptin level of\nexcitatory input to the system, which is then distributed\nto the populations of glutamatergic neurons and GABA\ninterneurons in accordance with the relative glutamatergic\nexcitation ratio parameter α ∈ [0, 1] representing the\nproportion of input directed to glutamatergic neuronal\npopulation. We minimised the number of inhibitory\ncoupling strength parameters by setting self-inhibition in\nthe GABAergic populations and functional interactions\nbetween GABAergic efferent neurons and the other two\npopulations to zero. In the absence of data that specifically\nsupports the inclusion of such inhibitory interactions, a\nmodel with fewer parameters is justified and easier to\ninterpret. As one of our aims is to investigate effects\nof GABA and glutamate receptor antagonism following\npharamcological interventions, we incorporate the terms\n(1 − β1) and (1 − β2), where β1 and β2 represent the\nproportion of suppressed functional interaction between\nGABA and glutamate neuronal populations, respectively.\nUsing the stimulus-response function and the proposed\ninteractions between the populations (Figure 1(b)) the\naveraged activity in the populations is governed by the\nfollowing functions:\nfl(Gl, Gi, Ge) = − Gl + (1 − Gl)ϕ(al, Fl(Gl, Gi, Ge), θl),\n(2.5)\nfi(Gl, Gi, Ge) = − Gi + (1 − Gi)ϕ(ai, Fi(Gl, Gi, Ge), θi),\n(2.6)\nfe(Gl, Gi, Ge) = − Ge + (1 − Ge)ϕ(ae, Fe(Gl, Gi, Ge), θe),\n(2.7)\nwhere the dependent variables Gl, Gi, Ge represent the\nmean activity in the populations of glutamatergic neurons,\nGABA interneurons and GABAergic efferent neurons at\ntime t, respectively. The model also includes refractory\ndynamics via the term (1 − G), which controls the time\nperiod during which the populations are unable to produce\na signal following an activation, and its primary effect is\ndecreasing the maximum firing rate [29]. The MePD activity\nis then governed by the following ordinary differential\nequations (ODEs):\ndGl\ndt =fl(Gl, Gi, Ge), (2.8)\ndGi\ndt =fi(Gl, Gi, Ge), (2.9)\ndGe\ndt =fe(Gl, Gi, Ge). (2.10)\nThe presented MePD system is non-dimensional w.r.t. time\n[30]. To effectively couple the system, we introduce the\ntime scaling factor d that relates arbitrary time in equations\n2.8-2.10 to time in minutes:\nt = T · δ, (2.11)\nwhere t is the original arbitrary time, T is the new time\nmeasured in minutes, and δ is the scaling factor (min −1).\nThe time-converted version of the model is as follows:\ndGl\ndT =δfl(Gl, Gi, Ge), (2.12)\ndGi\ndT =δfi(Gl, Gi, Ge), (2.13)\ndGe\ndT =δfe(Gl, Gi, Ge). (2.14)\nThe parameter values for the MePD model can be found\nin Table 1.\n(f) Calculating MePD output in the mean-field\nmodel\nThe magnitude of the mean glutamatergic and GABAergic\nMePD projections’ output is found in the same way, i.e. as\nthe integral of mean activity in the respective population\nover the integration time period T :\nMean glutamate output = 1\nT\nZ T\n0\n(Gl(t)dt, (2.15)\nMean GABAergic output = 1\nT\nZ T\n0\n(Ge(t)dt. (2.16)\nTo quantify the periodic output of the MePD GABA-\nglutamate circuit we compute the integral of the difference\nof mean activity in the populations of glutamatergic\nneurons (Gl) and GABAergic efferent neurons (Ge) over the\nintegration time period T :\nMean MePD output = 1\nT\nZ T\n0\n(Gl(t) − Ge(t))dt, (2.17)\nwhere the term 1\nT is the reciprocal of the time duration,\nallowing to normalise the output.\n(g) Coarse-grained Model of ARC KNDy\npopulation with MePD input\nBased on experimental evidence regarding MePD\nprojections to other brain regions including the ARC\n[14, 15] we couple our MePD model’s output with our\nARC KNDy network model [6, 7] aiming to explore the\neffects ofperturbations to the MePD circuit on GnRH\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n5ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\nTable 1. MePD model parameter values\nParameter Definition Value Reference\n..........................................................................................................................................................................................................................................................................................\nδ Temporal scaling factor (min−1) 3 [31]\nα relative glutamatergic excitation ratio (a.u.) 0.9 Derived\nKp Excitatory input to the MePD circuit (a.u.) [0, 14]\ncll glutamatergic self-excitation strength (Glu → Glu) (a.u.) 18 Derived\ncli interaction strength Glu → GABAint (a.u.) 16 Derived\ncil interaction strength GABAint → Glu (a.u.) 35 Derived\ncle glutamatergic excitation of GABA efferents (Glu → GABAeff) (a.u.) 40 Derived\ncie GABAergic inhibition of GABA efferent neurons (GABAint → GABAeff) (a.u.) 25 Derived\nal maximum slope of Glu (a.u.) 1.3 [24]\nai maximum slope of GABAint (a.u.) 2 [24]\nae maximum slope of GABAeff (a.u.) 2 [24]\nθl half-maximum firing threshold for Glu (a.u.) 4 [24]\nθi half-maximum firing threshold for GABAint (a.u.) 3.7 [24]\nθe half-maximum firing threshold for GABAeff (a.u.) 3.7 [24]\nβ1 GABAergic interaction suppression coefficient (a.u.) [0, 1]\nβ2 glutamatergic interaction suppression coefficient (a.u.) [0, 1]\npulse generator activity and to validate our model against\nexperimental observations in [12, 13]. The model describing\nthe dynamics in the KNDy neuronal network is given by\nthe following system of ordinary differential equations:\ndD\ndT =fD(v) − dDD, (2.18)\ndN\ndT =fN (N, v) − dN N, (2.19)\ndv\ndT =fv(N, v) − dvv, (2.20)\nwhere D and N represent the concentration of Dynorphin\nand Neurokinin B produced by the population, and v\ndescribes the averaged firing activity in the population in\nspikes/min. Parameters dD,dN , and dv control the linear\ndecay for each variable. Dynorphin and Neurokinin B\nsecretion rates are represented by functions fD and fN ,\nrespectively, while fv describes how the firing rate changes\nin response to the Neurokinin B concentration and current\nfiring rate. The neuropeptides’ secretion rates are given by\nthe following functions:\nfD(v) =kD\nv2\nv2 + K2\nv,1\n, (2.21)\nfN (N, v) =kN\nv2\nv2 + K2\nv,2\nK2\nD\nK2\nD + D2 , (2.22)\nwhere kD and kN signify the neuropeptides’ secretion rates;\nKv,1 and Kv,2 describe the frequency value for which\nthe rate of Dynorphin and Neurokinin B secretion is half-\nmaximum; and KD describes the Dynorphin concentration\nthat results in half-maximum inhibition of NKB. In the\noriginal introduction of the KNDy model [6], the function\nfv can take both positive and negative values. Here we\nmodified fv by restricting its output to be non-negative via\nvertical shift of the sigmoid function:\nfv(N, v) =v0\n1\n1 + exp (k(−I + m)) (2.23)\n(2.24)\nwhere v0 is the maximum increase of the firing rate in\nresponse to synaptic inputs I (Hz). The parameter m\nsignifies the synaptic input level at which the increase in\nthe firing rate becomes half-maximum. The parameter k\nrepresents the membrane’s time constant, which determines\nhow quickly the neuron’s membrane potential changes in\nresponse to inputs.\nThe synaptic inputs I consist of both external and\ninternal contributions:\nI = I0 + pv\nN 2\nN 2 + K2\nN\nv + jlGl − jeGe, (2.25)\nwhere I0 stands for the basal input in the population.\nThe excitatory effect of Neurokinin B on the firing rate is\naccounted via a sigmoid function with KN representing\nNeurokinin B’s half-maximal effect and pv control’s the\nstrength of the connection between the neurons in the\nKNDy. To account for the effects of the MePD output\non the KNDy, we use the terms jlGl and −jeGe,\nwhich signify the excitatory glutamatergic and inhibitory\nGABAergic contribution from the MePD, respectively.\nThe parameters jl and je are presynaptic firing rate\nconversion parameters for the corresponding populations\nin the MePD. These parameters allow converting non-\ndimensional output from the MePD to the input to the\nKNDy in Hz and assign weight to the contribution and,\nwe alter when simulating the effects of MePD projections’\nstimulation. In the case of simulating the effects of\nneurotransmitter antagonism and kisspeptin stimulation\nin the MePD we set jl = je = 1. To mimic the effects of\nstimulating glutamatergic and GABAergic projections, we\nconsider varied levels of MePD network excitation ( Kp)\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n6ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\nand alter the presynaptic firing rate conversion parameters\nto change the weight distribution of the projections. In\ncase of glutamatergic projections stimulation we increase\nthe weight of the glutamatergic contribution and decrease\nthe GABAergic weight contribution (j l = 1.5, je = 0.5),\nwhile for GABA projections we increase the weight of the\nGABAergic contribution and decrease the glutamatergic\nweight contribution (j l = 0.5, je = 1.5). For further details\non the original KNDy model, refer to [6]. The parameters\nfor the KNDy network can be found in Table 2.\n(h) Numerical simulations and bifurcation\nanalysis.\nBifurcation analysis was performed in AUTO 07-p [36],\nwhile numerical simulations were carried out in MATLAB\nusing ode45 (Runge-Kutta method) for the MePD system\nand ode15s (variable-step, variable-order (VSVO) solver) for\nthe coupled MePD-KNDy model. The codes for reproducing\nthe analysis and simulations presented in this manuscript\ncan be found in GitHub repository.\n3. Results\n(a) Modelling the MePD’s GABA-glutamate\ncircuit and its projections to the ARC\nThe diagram depicted in Figure 1(a) provides an overview\nof our model describing the functional connectivity in the\nMePD neural circuit along with the MePD’s projections\n(outputs) to the ARC. In our mode we consider an\nexcitatory population of glutamatergic neurons and an\ninhibitory population of GABAergic neurons given the\nexperimentally-established presence of glutamatergic and\nGABAergic neurons in the MePD [22,23]. These populations\nof glutamatergic and GABAergic neurons interact with\neach other, and also extend excitatory and inhibitory\nconnections, respectively, to a distinct neuronal population\nof GABAergic neurons [26, 27], which we refer to here as\nGABA efferent neurons. The activity in the populations\nof the GABA efferent neurons and glutamatergic neurons\ndefines the MePD’s output that we consider in the\nmodel to be acting on the ARC. It has previously\nbeen shown by [13] that a kisspeptin-expressing neuronal\npopulation, found in the MePD [37], provides excitatory\ninput to the populations of GABA interneurons and\nglutamatergic neurons and has been accordingly included\nin our modelling. Our mathematical model is based on\nthe established Wilson-Cowan modelling framework [24,\n25] and hence allows us to simulate the mean activity of\nthe different neuronal populations; namely glutamatergic\nneurons (Glu), GABA interneurons (GABA int) and GABA\nefferent neurons (GABAeff).\nIn this study, we couple the MePD neuronal network\nmodel to our KNDy neuronal network model [6, 7].\nIn previous work, we coupled a first-generation model\nof the MePD neuronal network with our KNDy (pulse\ngenerator) model and performed numerical simulations in\norder to reproduce the results of optogenetic stimulation\nof kisspeptin and pharmacological antagonism experiments\nin the MePD in [13]. The MePD neuronal network model\nused in [13] was based on the same framework as the model\nin this manuscript, but under the assumption of stationary\nMePD network activity and hence constant MePD output.\nIn the present study, we investigate the proposed MePD\nneuronal network in more detail, taking into account the\npossibility of dynamic (e.g. oscillatory) MePD activity. To\nthis end, we match the temporal activity in the circuit to the\ntime scales of calcium activity recorded in MePD neurons\n[31]. Such activity is now routinely used as a proxy of\nmean neuronal activity, and in our case it is mediated by\nGABA and glutamate neuronal populations in the MePD.\nFull details of the model are given in Mean-Field Model of\nthe MePD.\n(b) How does excitatory input decrease\ninhibitory tone in the MePD circuit?\nOptogenetic stimulation of MePD kisspeptin neurons has\nbeen shown to have a significant effect on LH pulses\n[12], presumably via exciting GABA-glutamate neuronal\ncircuits and their projections to the ARC. To investigate how\nthis effect could be relayed through the GABA-glutamate\nMePD neuronal network, we study the model’s behaviour\nunder various levels of excitatory kisspeptin input. Previous\nanalysis of the Wilson-Cowan model [24, 25, 30] has shown\nthat oscillatory dynamics in the model can be induced\nvia glutamatergic self-excitation and a negative-feedback\nloop between the populations of GABA interneurons and\nglutamate neurons. Accordingly, in our MePD model,\nwe assume excitatory coupling from the glutamatergic\npopulation to the population of GABA interneurons as\nwell as inhibitory coupling from the GABA interneurons\nto the glutamatergic population. To mimic experimental\noptogenetic stimulation of the GABA-glutamate network,\nwe perform a bifurcation analysis using the level of\nMePD kisspeptin (\"MePD excitation\" in Figure 2) as a free\nparameter (for further details on the numerical methods, see\nNumerical simulations and bifurcation analysis).\nThe analysis reveals that at low kisspeptin excitation,\nthe activity of all neuronal populations considered in our\nmodel is low and exhibits stationary dynamics (Figure 2( a-\nc)). This could be explained by the fact that the system does\nnot receive enough excitation to sustain oscillations, and as\na result, settles in a state of low mean activity where each\npopulation behaves as a pool of independent single neuron\noscillators. As excitation increases, the activity in all three\npopulations is amplified. Numerical continuation along the\nstable equilibrium branch reveals that the system undergoes\na change in qualitative dynamics (atKp = 1.4) due to a Hopf\nbifurcation (HB), giving rise to a branch of stable periodic\n(limit cycle) solutions. Within the parameter range where\nthe limit cycle solutions exist, at a threshold kisspeptin level\nof excitatory input (K p = 1.6) the gain in the population of\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n7ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\nTable 2. KNDy model parameter values\nParameter Definition Value Reference\n..........................................................................................................................................................................................................................................................................................\ndD Dyn degradation rate (min−1) 0.2 [6]\ndN NKB degradation rate (min−1) 1 [6]\ndv Firing rate reset rate (min−1) 10 [32]\nkD Dyn signalling strength (nMmin−1) 4 [7]\nkN NKB signalling strength (nMmin−1) 40 [7]\npv Effective strength of synaptic input (a.u.) 0.006 [7]\nv0 Maximum rate of neuronal activity increase (spikes min−2) 25000 [32]\nKD Dyn IC50 (nM) 0.3 [33]\nKN NKB IC50 (nM) 4 [34]\nKv,1 Firing rate for half-maximal Dyn secretion (spikes min−1) 600 [35]\nKv,2 Firing rate for half-maximal NKB secretion (spikes min−1) 200 [35]\nk Membrane’s time constant (min) 10 Fixed\nm Half-maximal firing rate (min−1) 0.5 Fixed\nI0 Basal activity (min−1) 0.14 Fixed\njl presynaptic firing rate conversion parameter for glutamatergic projections (Hz) {1, 1.5, 0.5}\nje presynaptic firing rate conversion parameter for GABAergic projections (Hz) {1, 0.5, 1.5}\nGABA efferent neurons switches from positive to negative\n(Figure 2(c)). As a result, further increase in excitation\nleads to decrease in the activity of the GABA efferent\nneuronal population, which approaches zero with further\nincrease in the kisspeptin excitatory input (see Figure 2(c)).\nThis can be explained by the fact that inhibitory input\nfrom GABA interneurons to GABA efferents outweighs\nthe excitatory input from the glutamatergic population\n(Figure 2(a-b)). The presence of a negative feedback loop\nin the system leads to a higher rate of increase in activity\nof the GABA interneuron population compared to the\npopulation of glutamate neurons. Therefore, by comparing\nthe activity of the populations of glutamatergic neurons and\nGABA efferents at different excitation levels (Figure 2( d-\ne)) we observe that, overall, the MePD projections’ output\nwould increase due to the reduction in the inhibitory\nGABAergic tone. In the parameter range where the model\nexhibits periodic behaviour, we also find that the oscillatory\nperiod decreases as we increase excitation (Figure 2( f)).\nThe numerical range of the oscillation period confirms\nthat the temporal activity in the model aligns well with\nthe experimentally observed average calcium oscillation\nperiod reported in [31], where the oscillations in MePD\nneuronal activity occur, approximately, in the span of\na minute. We further extend the bifurcation diagram\n(see Figure S1), identifying that the further increase in\nexcitation leads to an exponential increase in the oscillation\nperiod and subsequent destruction of the limit cycle\nsolution via a global homoclinic bifurcation. Following\nfurther increase in the kisspeptin level, the system enters\na bistable regime, induced and destroyed via saddle-\nnode bifurcations, followed by constant a high population\nactivity mode.\nIt is unknown whether MePD kisspeptin directly\nmodulates GABA interneurons and/or glutamatergic\nneurons. To investigate the role of the distribution of\nexcitatory input between the inhibitory (GABA) and\nexcitatory (glutamate) populations in the model, we\nintroduce a parameter that controls the relative kisspeptin\nexcitation ratio (α). We are thus able to continue the loci of\nthe co-dimension one Hopf and saddle-node bifurcations\n(electronic supplementary material, figure S1) in two-\nparameter space (namely, the excitatory input Kp and\nthe relative kisspeptin excitation ratio α) (Figure 2(g)).\nThe Hopf bifurcation and saddle-node bifurcation curves\ncoalesce in the two-parameter space, where a co-dimension\ntwo Bogdanov-Takens (BT) point emerges. This (BT)\npoint is also related to the appearance of a homoclinic\nbifurcation curve, representing a situation where the stable\nand unstable manifolds of a saddle equilibrium intersect,\nindicating the presence of complex dynamical behaviour in\nthe system. The bifurcation curves representing the Hopf\n(HB), saddle node (SN) and homoclinic (HC) loci allow us\nto identify regions in two-parameter space characterised\nby different qualitative dynamics in the system; and in\nparticular, regions where the system oscillates. We note\nthat for the current choice of parameters (Table 1), the\noscillations in the system can be induced only under\nthe condition that the majority of excitation is directed\nto the excitatory population of glutamatergic neurons.\nTo investigate how the distribution and different levels\nof excitation affect the output of the system within the\noscillatory region, we compute a heat map depicting\nchanges in the mean MePD projections’ output (Figure 2(g))\ndue to changes in activity of the populations of excitatory\n(glutamate) neurons and inhibitory (GABA) efferent\nneurons in our MePD model. For complete details on\nhow the mean MePD projections’ output is defined and\ncomputed, see Calculating MePD output in the mean-\nfield model. We find that an increase in the level of MePD\nkisspeptin excitation leads to a transition from inhibitory\nto excitatory MePD output. As the proportion of kisspeptin\nexcitation to the glutamatergic population is increased, the\ninhibitory tone of the MePD circuit can be maintained under\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n8ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\nstronger excitation in the model (Figure 2(g)). This suggests\nthat additional (kisspeptin) excitation of the glutamatergic\npopulation may lead to an increase in GABAergic tone,\ndepending on the functional interaction strength between\nthe glutamatergic and GABAergic neuronal populations.\nTaken together, our theoretical findings suggest that the\nreduction in GABA efferent neuron activity amid increased\nexcitation of the MePD neuronal circuit may be reliant on\nthe intricate interplay between competitive inhibitory and\nexcitatory connections to the population of GABA efferent\nneurons in the MePD circuit.\n(c) How does MePD functional network\nconnectivity affect MePD projections’\noutput?\nThe presence of GABAergic and glutamatergic neuronal\npopulations in the MePD [22,23], as well as their importance\nin the modulation of GnRH pulse generator activity,\nhas been demonstrated experimentally [13]. However, the\nrole of their functional interactions within the circuit\nremains unknown. Hence, in this section we investigate\nhow changes in the functional interaction (coupling)\nstrength affect the dynamics of the GABA-glutamate\ncircuit and the MePD output. The aim here is to\ncharacterise network interaction patterns associated with\noscillatory dynamics and the corresponding periodic MePD\nexcitatory/inhibitory projections’ output.\nThe stimulatory effect of the amygdala on the GnRH\npulse generator under the stimulation of kisspeptin\nneurons has been attributed to the activation of GABAergic\ninterneurones, which in turn inhibit GABAergic efferent\nneurons [13], forming a GABA-GABA disinhibitory\ninteraction, which is of interest in understanding functional\nmechanisms and dynamics in the MePD circuit. As\ncompeting excitatory and inhibitory signals counterbalance\neach other, we fix the interaction strength responsible\nfor glutamate input to GABA efferent neurons and\ninvestigate how MePD projections’ output changes under\nthe variation of kisspeptin (level) stimulation and the\nstrength of interaction between GABA interneurons and\nGABA efferents (c ie) Figure 3. The mean glutamatergic\nprojections’ output remains relatively unaffected by the\nstrength of interaction due to the lack of GABA efferent\nneuronal input to the glutamate neuronal population,\nbut its magnitude moderately increases as the excitation\nin the circuit increases (Figure 3(a)). Meanwhile, as the\nsystem is excited, a low strength of interaction results\nin an increase in the magnitude of the mean GABAergic\nprojections’ output (Figure 3(b)), as the GABA-GABA\ninteraction is not sufficient to induce a decrease in the\nactivity of the population of GABA efferent neurons under\nthe increased excitation in the MePD circuit. On the other\nhand, under high interaction strength, the population of\nGABA interneurons exerts an excessive inhibitory input\nto the population of GABA efferent neurons, resulting in a\nreduction in the GABA efferent neuronal population activity\nto zero (Figure 3(b)). We compute the combined MePD\nprojections’ output as the difference in magnitude between\nthe glutamatergic and the GABAergic projections’ outputs.\nWe show that a low level of functional interaction strength\nbetween GABA interneurons and GABA efferents leads\nto a predominantly inhibitory mean MePD projections’\noutput, whereas a high interaction strength results in an\nexcitatory mean MePD projection’s output that remains\nmostly unchanged under further excitation of the MePD\ncircuit (Figure 3(c)). This happens as the change in the\nactivity of the population of GABA efferent neurons is much\nhigher compared to the population of glutamate neurons\n(Figure 3(a-b)). Hence, the mean MePD projections’ output\nis heavily dependent on the strength of GABA-GABA\ndisinhibition.\nAnother critical component of the system is the\nfunctional interaction strength between the populations\nof GABA interneurons and glutamate neurons, which\nfacilitates oscillatory behaviour in the model by providing\nnegative feedback between the two populations. Having\nfixed the kisspeptin excitation level to induce oscillatory\ndynamics (Kp = 2.3), we perform one parameter bifurcation\nanalysis using the functional interaction strength between\nGABA interneurons and glutamate neurons (c il) as a\nbifurcation parameter. Our analysis reveals Hopf ( HB),\nsaddle node (SN) and homoclinic (HC) bifurcations that\nwe then continue in two parameters (using the functional\ninteraction strength between glutamate neurons and GABA\ninterneurons (c li) as a second bifurcation parameter). This\nenables us to investigate the qualitative dynamics of the\nmodel under the variation of GABA-glutamate interaction\n(Figure 4(a)). The intersections of the Hopf and saddle\nnode curves are associated with the location of Bogdanov-\nTakens (BT) points, which also gives rise to the homoclinic\ncurve, while the location where two saddle node curves\nmeet tangentially indicates the location of a co-dimension\ntwo cusp (CP) point. Our analysis confirms that the\nexistence of oscillatory dynamics in the MePD neuronal\nnetwork intricately depends on the balance between the\ninhibitory and excitatory interaction strengths. Specifically,\ninsufficient or excessive strength of functional interaction\nbetween glutamatergic neurons and the population of\nGABA interneurons causes loss of oscillatory dynamics,\nwhile high inhibitory strength of interaction does not\nprevent oscillations, but rather makes the mean MePD\nprojections’ output more inhibitory (Figure 4(a)).\nNext, we consider the combined effects of the\nglutamatergic population’s self-excitation and its excitatory\ninput to the population of GABA interneurons. We perform\none-parameter bifurcation analysis using the strength of\nfunctional interaction between glutamate neurons and\nGABA interneurons (c li) as a bifurcation parameter,\nidentifying Hopf (HB) bifurcation points (Figure 4(b)).\nIncreasing the strength of functional interaction between\nglutamatergic population and GABA interneurons causes a\ndecrease in glutamatergic population activity as the effect of\nthe negative feedback loop is amplified, resulting in higher\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n9ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\ninhibitory interaction between GABA interneurons and the\nglutamatergic population. We then continue the detected\n(HB) points in two parameters (glutamatergic self-excitation\nstrength (cll) and strength of functional interaction between\nglutamate neurons and GABA interneurons (c li)) to\ndefine the oscillatory region in two-parameter space\n(Figure 4(b)). We find that higher levels of self-excitation\nrequire higher strength of functional interaction between\nglutamate neurons and GABA interneurons in order to\ngive rise to oscillatory dynamics. Increasing the strength\nof glutamatergic functional interaction switches the mean\nMePD projections’ output from negative to positive.\nThe above analysis indicates that the intricate\nexcitatory/inhibitory MePD projections’ balance is heavily\ndependent on the MePD circuit’s functional interactions. It\nalso suggests that coordinated changes in the interaction\nstrengths between the circuit’s neuronal populations may\nbe a critical regulator of the MePD projections’ output, and\nthus its modulatory effect on the GnRH pulse generator.\n(d) MePD projections’ dynamic modulation of\nGnRH pulse generator activity\nHaving analysed the MePD GABA-glutamate neuronal\nnetwork behaviour and the effects of different model\nparameters on the mean MePD projections’ output, we\nnow investigate our coupled MePD-KNDy network model\n(see Coarse-grained Model of ARC KNDy population with\nMePD input for full model details). The aim here is to\ncharacterising the differential effects of MePD dynamic\nprojections’ output on GnRH pulsatility. In essence,\ncoupling the MePD and KNDy models results in feeding\nexternal periodic input from the MePD neuronal network\nwith the KNDy neuronal network (a.k.a. GnRH pulse\ngenerator). As the timescales of the two network models are\nsignificantly different (MePD neuronal network evolves on\na timescale of seconds while the KNDy neuronal network\noperates on a timescale of minutes), in the coupled model\nthere is more than one frequency found in the periodic\ntrajectory which now evolves on a torus, i.e. the limit cycle\nbecomes a limit torus solution (Figure 5(a)). In the case when\nthe extended system has two incommensurate frequencies,\nthe trajectory is no longer closed, leading to quasi-periodic\ndynamics. This is not surprising as it is well established\nthat in response to periodic input, relaxation oscillators\nsuch as the KNDy network model [6] can exhibit complex\ndynamics, like quasi-periodicity [38, 39].\nNext, we validate the coupled model by reproducing\nin vivo experiments where the effect of optogenetic\nstimulation of MePD kisspeptin neurons on LH pulse\nfrequency was investigated [12]. To simulate the effects of\noptogenetic stimulation, we increase the kisspeptin level of\nexcitation within the MePD, which results in activation of\nthe GABA-GABA component of the MePD neuronal circuit,\nresulting in a decrease in the activity of the population\nof GABA efferent neurons’(Figure 5(b)). Consequently, the\nreduction of the inhibitory tone in the MePD projections’\noutput under increased excitation promotes a decrease\nin the interpulse interval (IPI) in the KNDy system\n(Figure 5(b)).\nTo investigate the effects of suppression of GABAergic\ninteraction strength on the system’s dynamics, we compute\na two-parameter bifurcation diagram for a range of MePD\nexcitation (K p) and GABA functional interaction strength\nsuppression coefficient β1 (see Figure 1(b)), the latter\ndescribing the strength of GABA receptor antagonism\n(Figure 5(c)). We find that complete suppression of\nGABAergic interaction leads to the loss of oscillatory\ndynamics in the MePD circuit via its effect on the\nnegative feedback loop between the populations of GABA\ninterneurons and glutamatergic neurons. It is common\nthat the effect of pharmacological blockers is modelled by\ncomplete suppression of functional interactions [13], but in\nreality only partial suppression may occur. Here, we show\nthat partial blocking of GABAergic interaction is sufficient\nto decrease the MePD projections’ output under increased\nexcitation (Figure 5(c)). On the other hand, reduction\nof GABAergic functional interaction strength leads to a\ndecrease in the inhibitory coupling between the population\nof GABA interneurons and the glutamatergic population\nand the population of GABA efferent neurons, hence\nincreasing both glutamatergic and GABAergic tone in the\nMePD (Figure 5(d)). However, the difference between the\ninhibitory and excitatory tone remains relatively constant\nas before the suppression, resulting in an unperturbed\nKNDy network interpulse interval. Combining suppression\nof GABAergic interactions with increased excitation in the\nMePD also increases the activity in the populations of\nglutamatergic neurons and GABAergic efferent neurons,\nbut also amplifies the difference between inhibitory and\nexcitatory mean MePD projections’ output (Figure 5(e)). As\na result, the MePD projections’ input to KNDy becomes\nmore inhibitory, leading to an increase in the KNDy network\ninterpulse interval.\nSimilarly, to mimic the effects of a glutamate receptor\nantagonist, we decrease the strength of glutamatergic\ninteractions in the model and observe that oscillatory\ndynamics rapidly cease (Figure 5(f)). We also observe\nthat under partial suppression of glutamatergic interaction\nstrength and increased excitation, there is a transition\nfrom inhibitory to excitatory tone in the system, while\ncomplete abolition of glutamatergic interactions causes the\nmean MePD output to be exclusively excitatory for all\nlevels of excitation (Figure 5(f)). To preserve inhibitory\noutput from the MePD under low levels of excitation,\nwe consider a partial suppression of the functional\ninteraction strength associated with the population of\nglutamatergic neurons. For lower levels of excitation (still\nenabling oscillatory dynamics), moderate suppression of\nglutamatergic interaction strength results in a loss of\noscillatory MePD dynamics, while keeping the mean\nMePD output relatively constant (and inhibitory); hence\nthe extended system dynamics transitions from evolving\non a limit torus to a limit cycle while the IPI of the\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n10ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\nKNDy system remains unaffected (Figure 5(g)). Suppression\nof glutamatergic functional interaction strength combined\nwith an overall increase in excitation of the MePD network\ncauses a switch from inhibitory to excitatory mean MePD\noutput (due to the presence of GABA-GABA interactions).\nThis results in over-stimulation of the KNDy network\nassociated with a transition from a pulsatile to a quiescent\nmode of operation due to a depolarisation-block-type of\nphenomenon (Figure 5(h)).\n(e) How stimulation of MePD projections\nmodulates GnRH pulse generator activity\nThe suppressive effect of MePD GABA projections’ output\non the GnRH pulse generator frequency has been recently\nconfirmed experimentally by optogenetically stimulating\nMePD GABAergic terminals in the ARC [16]. Here, using\nour coupled model, we interrogate the role of both\nGABAergic and glutamatergic projections of the MePD\ncircuit in modulating KNDy dynamics. In the previous\nsection the weights of GABA and glutamate projection\noutput were set to be equal (j l = je = 1). Here in order\nto account for direct optogenetic simulation of these\nprojections as carried out in [16] we increase the weight\nof the respective projection output in the model (i.e. for\nsimulating the experiments presented in [16] we increase\nthe weight of GABA projections output in the model)\nand vary the stimulation level in the MePD circuit ( Kp).\nSpecifically, if we set je = 1.5 and jl = 0.5 we are able\nto show that increasing the stimulation level (K p) in the\ncoupled system initially produces no change in the KNDy\nIPI, followed by an exponential increase in the KNDy inter\npulse interval corresponding to cessation of the LH pulsatile\ndynamics (Figure 6( a)). These results conform with in vivo\nstimulation of MePD GABA projections, where a decrease\nin LH pulsatility (demonstrated by increase in inter pulse\ninterval) eventually leading to loss of pulsatile dynamics\nhave been observed at 10 and 20 Hz stimulation of GABA\nprojections, respectively [16].\nNow, we employ the same strategy to simulating the\neffects of glutamate projections’ stimulation (by fixing\njl = 1.5 and je = 0.5). Our model simulations show that\nunder increasing excitation of the MePD network, the\ninter pulse interval of KNDy population activity initially\ndecreases before increasing and then returning to its\ninitial IPI (Figure 6(a)). Furthermore, bifurcation analysis\nof the extended model demonstrates that the system\nundergoes a torus bifurcation (TR), associated with the\nswitch from limit cycle dynamics to dynamics evolving\non a limit torus, which we depict in (Figure 5(a)). The\nobserved non-monotonic behaviour of the coupled system\nin this case is counterintuitive, given the excitatory role\nof glutamate. Nevertheless, selective in-vivo optogenetic\nstimulation of the MePD glutamatergic projections in the\nARC with increasing levels of stimulation confirms our\nmodel predictions as shown in (Figure 6(b)). As expected,\ngiven the excitatory nature of the glutamatergic projections,\nsustained stimulation at 5 Hz results in a significant\ndecrease in LH interpulse interval from 15.42 ± 0.60 min\nto 11.81 ± 0.81 min. However, further increase in the\nfrequency of stimulation (at 10 Hz and 20 Hz) restores the\npre-stimulation IPI levels, as predicted by our modelling.\nMoreover, our modelling allows us to explore potential\nmechanisms that govern the non-monotonic response\nin the GnRH pulse generator found experimentally as\ndescribed above (Figure 6(c-d)). Continuation analysis of\nthe MePD model dynamics indicates that the decrease in\nthe KNDy network IPI observed in the model is due to\nthe amplification of glutamatergic activity in the MePD\nnetwork, while GABAergic tone remains very close to zero\n(see Figure 6(c) at stimulation level ≈ 1.2). Further increase\nin MePD network excitation, however, switches the balance\nin excitatory/inhibitory MePD projections’ output (i.e.\nMePD input to KNDy)(Figure 6(d), which in turn promotes\nan increase in the KNDy IPI. The MePD network model\nundergoes a Hopf bifurcation (HB), the location of which is\nassociated with the location of the torus bifurcation (TR) in\nthe extended MePD-KNDy network model (Figure 5( a-b)),\ndemonstrating that the transition from limit cycle to limit\ntorus occurs due to a change in the qualitative dynamics of\nthe MePD network model.\n4. Discussion\nIn our study, we have introduced and systematically\ninvestigated a model incorporating the interplay between\nGABA and glutamate neuronal populations within the\nMePD. This model was coupled to a GnRH pulse generator\nmodel [6, 7], allowing us to validate it against experimental\nfindings from [12, 13] as well as offering insights into how\nperturbations in the MePD could impact the activity of the\nGnRH pulse generator. Our model could serve as a versatile\ntool for investigating broader MePD circuit effects, such as\nfor example those stemming from the interactions between\nurocortin and the GABA/glutamate neuronal populations\n[40]. The utility of our modelling approach lies in its\nability to interrogate the neuronal mechanisms that enable\nthe MePD to modulate the dynamics of the GnRH pulse\ngenerator, hence enabling us to better understand the effects\nof environmental and psychosocial factors on reproductive\nfunction, such as pubertal timing [9–11] and modulation of\nLH secretion [12, 13].\nA model of the MePD circuit coupled to the GnRH\npulse generator has been previously studied under\nstationary MePD circuit dynamics [13]. In this mode,\nthe system functions like a collection of independent\nneuronal oscillators characterised by a constant (averaged)\npopulation level of activity. However, it is important to\nnote that while this study has offered valuable insights,\nthe actual patterns of MePD activity are likely to be more\ncomplex. Indeed, [31] shows changes in the MePD neuronal\nnetwork’s oscillatory activity in vivo associated with sex-\nspecific differences in the encoding of social stimuli and\nsexual experience. Here, we demonstrate that the extended\nmodel is able to reproduce experimental findings [12, 13],\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n11ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\nsuggesting the plausibility of an oscillatory mode of MePD\ncircuit activity. In fact, rhythmicity of neuronal populations\nis a characteristic feature of neuronal synchronisation,\nallowing the neuronal networks to manage and process\ncomplex stimuli [41]. However, to confirm or reject the\nhypothesis about the importance of oscillations in the\nMePD neuronal networks, further experiments involving\nrecordings of calcium activity in individual GABA and/or\nglutamate neurons in the MePD, and how they synchronise,\nwould be required.\nUnder oscillatory MePD circuit behaviour, our modelling\nshows that an increase in the excitatory input to the\nMePD system decreases GABAergic MePD output due\nto the activation of GABA-GABA disinhibition, while\nglutamatergic output remains consistent. This finding\nindicates that GABAergic MePD output is sensitive to\nstimulatory inputs, while glutamatergic output is likely\nto play more of a balancing role. This is consistent\nwith the established role of amygdala in reproductive\nfunction modulation, as lesions to the MePD have been\nshown to advance puberty [42] and prevent stress-induced\nsuppression of LH pulses [43]. On the other hand,\noptogenetic stimulation of kisspeptin neurons in the MePD\n[12] increases LH pulse frequency, and administration of\nperipheral kisspeptin inhibits neuronal activation in the\namygdala as well as increases LH secretion [44], which can\nbe explained by a decrease in the activity of GABAergic\nefferent neurons.\nHere, we have modelled the effects of pharmacological\ninterventions via partial suppression of MePD circuit\nfunctional interactions and studied how different levels of\nsuppression affect qualitative dynamics in the model; this\nis in contrast to our previous work [13] where we assumed\ncomplete suppression of signalling. In reality, however,\ncomplete suppression is unlikely to be the case, as neuronal\ncells may respond to pharmacological interventions by\nupregulating receptors or modifying their signalling\npathways to compensate for the inhibited receptors.\nWhen modelling partial GABA signalling suppression, the\ninhibitory component of MePD output increases because\nthere is not enough GABA-GABA disinhibition, but at\nthe same time glutamatergic activity also goes up due\nto decreased inhibition from the GABAergic population,\nbalancing out the inhibitory output. This compensation\nmechanism could provide an alternative explanation as to\nwhy solely GABA receptor antagonism does not change LH\npulsatility [13].\nFeeding oscillatory MePD projections’ output into the\nKNDy network allowed us to consider dynamic upstream\nmodulation of the KNDy network rather than a constant\ninput as considered in our previous work [7, 13]. We\ninvestigated how such dynamic input changes the response\nof the KNDy relaxation oscillator. Specifically, we have\ndemonstrated that when the KNDy relaxation oscillator\nreceives a dynamic input of a significantly different\nfrequency, this can result in a complex quasi-periodic\npulse pattern, i.e. irregular-shaped pulses with no change\nin interpulse interval. Quasi-periodicity has also been\nidentified in other relaxation oscillator systems subject to\nperiodic inputs, indicating that systems characterised by\nquasi-periodic behaviour often possess a level of resilience\nagainst external perturbations [38, 39]. On the other hand,\nthe transition from oscillatory to constant input, which we\nobserve during the suppression of glutamatergic functional\ninteractions, makes the KNDy system far more sensitive\nto the magnitude of the change. Previous modelling work\nsuggests high sensitivity of the KNDy network to the\nmagnitude of constant external stimuli that can lead to\ncessation of GnRH pulses [7]. However, to maintain a\nfunctional reproductive system, the GnRH pulse generator\nmust be resilient to small perturbations that can arise\nfrom changes in the MePD (or other upstream brain\nregions) output due to environmental stimuli, which is more\nplausible under periodic input modulation.\nPreviously published data identified that stimulation\nof MePD GABA projections in the ARC modulates\nGnRH pulse generator activity [16]. Using our extended\nmathematical model we interrogated the role of the\nglutamatergic projections in such modulation, which\nhas not been shown previously. The model analysis\npredicted a possible non-monotonic response of KNDy\nactivity. These model predictions were experimentally\nconfirmed in vivo, suggesting a novel excitatory MePD\nglutamatergic projection capable of impacting the KNDy\nnetwork, and hence modulating the GnRH pulse generator.\nSpecifically, we show that stimulation at 5 Hz leads to a\nstatistically significant decrease in LH interpulse interval,\nbut further increase in stimulation frequency (10 Hz and\n20 Hz) produces no significant change. Considering the\nexcitatory function of glutamate, one might question why\nan increase in stimulation does not lead to a further\ndecrease in LH IPI or even the complete loss of pulsatility\ndue to potential overstimulation (or in other words\ndepolarisation block). Using our extended mathematical\nmodel we reproduced this non-monotonic stimulus-\nresponse relationship, providing insight into potential\nmechanisms though which glutamatergic projections\nmodulate GnRH pulse generator activity. In the model,\nas the excitatory drive increases, so too does the response\nof excitatory projection neurons, leading to positive\ncorrelation between excitatory input and neuronal response,\nsimilar to the effect observed during 5 Hz stimulation\nof glutamatergic projections. However, when excitation\nbecomes stronger, which mimics the effects of the higher\nfrequency optogenetic stimulation, the system reaches a\npoint where further increase in the excitatory drive does not\nlead to a proportional increase in the neuronal response.\nExcessive stimulation may also trigger compensatory\nmechanisms that counteract the increased activity, leading\nto a limited net change [45]. As we observe in the\nmodel, the inhibitory projections could become engaged\nto maintain the balance, counteracting the excitatory drive\nand contributing to a non-monotonic response. Based\non the model analysis, we argue that balanced feed-\nforward excitation and feed-forward inhibition ensure that\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n12ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\nthe overall excitability of the KNDy network is robust\nto the elevated excitatory input, permitting GABAergic\nprojections to exert ’inhibitory brake’, thus constraining the\nLH IPI to pre-stimulation levels.\nHere we propose a mathematical model and explore\nthe influence of MePD activity on GnRH pulse generator\ndynamics. The modelling and analysis approach we have\nused could be applied to gain insight into the behaviour\nof other brain regions involved in modulation of the GnRH\npulse generator. Additionally, given the phenomenological\nnature of our MePD model, in the future, the extended\nMePD-KNDy model could be used to interrogate the effects\nof other stimulatory neuronal populations signalling to\nthe MePD GABA-glutamate circuit. This could then be\nemployed to perform in silico simulations to interpret\nand/or predict experimentally observed effects on the\nGnRH pulse generator and reproductive function.\nAppendix\nData Accessibility. The code to reproduce the analysis and data can be found in GitHub repository.\nAuthors’ Contributions. Kateryna Nechyporenko, Conceptualisation, Investigation, Methodology, Software, Writing – original draft;\nMargaritis Voliotis, Conceptualisation, Investigation, Writing – review and editing; Xiao Feng Li, Conceptualisation, Data curation,\nInvestigation; Owen Hollings, Conceptualisation, Data curation, Investigation; Deyana Ivanova, Conceptualisation, Investigation; Jamie\nWalker, Conceptualisation, Investigation, Writing - review and editing; Kevin O’Byrne, Conceptualisation, Investigation, Writing - review\nand editing; Krasimira Tsaneva-Atanasova, Conceptualisation, Investigation, Writing - review and editing.\nFunding. KTA gratefully acknowledges the financial support of the EPSRC via grant EP/T017856/1. KTA and MV gratefully\nacknowledge the financial support of the BBSRC via grant BB/W005883/1. KOB and XFLI gratefully acknowledge the financial support\nof the BBSRC via grant BB/W005913/1. JJW gratefully acknowledges the financial support of the MRC via grants MR/N008936/1 and\nMR/T032480/1.\nDisclaimer. For the purpose of open access, the author has applied a ‘Creative Commons Attribution (CC BY) licence to any Author\nAccepted Manuscript version arising.\nReferences\n1P . E. Belchetz, T. M. Plant, Y. Nakai, E. J. Keogh, and E. Knobil, “Hypophysial responses to continuous and intermittent\ndelivery of hypopthalamic gonadotropin-releasing hormone”, Science 202, Belchetz, P E Plant, T M Nakai, Y Keogh, E\nJ Knobil, E Journal Article United States 1978/11/10 Science. 1978 Nov 10;202(4368):631-3. doi: 10.1126/science.100883.,\n631–3 (1978).\n2A. E. Herbison, “Control of puberty onset and fertility by gonadotropin-releasing hormone neurons”, Nature Reviews\nEndocrinology 12, 452–466 (2016).\n3K. T. O’Byrne and E. Knobil, “Electrophysiological approaches to gonadotrophin releasing hormone pulse generator\nactivity in the rhesus monkey”, Hum Reprod 8 Suppl 2, 37–40 (1993).\n4P . Grachev and R. L. Goodman, “The GnRH Pulse Generator”,AIMS Medical Science 3, 359–385 (2016).\n5A. E. Herbison, “The Gonadotropin-Releasing Hormone Pulse Generator”, Endocrinology 159, 3723–3736 (2018).\n6M. Voliotis, X. F. Li, R. De Burgh, G. Lass, S. L. Lightman, K. T. O’Byrne, and K. Tsaneva-Atanasova, “The Origin of GnRH\nPulse Generation: An Integrative Mathematical-Experimental Approach”, The Journal of Neuroscience 39, 9738–9747\n(2019).\n7M. Voliotis, X. F. Li, R. A. De Burgh, G. Lass, D. Ivanova, C. McIntyre, K. O’Byrne, and K. Tsaneva-Atanasova,\n“Modulation of pulsatile GnRH dynamics across the ovarian cycle via changes in the network excitability and basal\nactivity of the arcuate kisspeptin network”, eLife 10, 10.7554/elife.71252 (2021).\n8J. Úbeda-Contreras, I. Marín-Blasco, R. Nadal, and A. Armario, “Brain c-fos expression patterns induced by emotional\nstressors differing in nature and intensity”, Brain Structure and Function 223, 2213–2227 (2018).\n9X. F. Li, D. A. Adekunbi, H. M. Alobaid, S. Li, M. Pilot, S. L. Lightman, and K. T. O’Byrne, “Role of the posterodorsal\nmedial amygdala in predator odour stress-induced puberty delay in female rats”, Journal of Neuroendocrinology 31,\ne12719 (2019).\n10D. Ivanova, X. Li, Y. Liu, C. McIntyre, C. Fernandes, G. Lass, L. Kong, and K. T. O’Byrne, “Role of Posterodorsal Medial\nAmygdala Urocortin-3 in Pubertal Timing in Female Mice”, Front Endocrinol (Lausanne) 13, 893029 (2022).\n11D. A. Adekunbi, X. F. Li, S. Li, O. A. Adegoke, B. O. Iranloye, A. O. Morakinyo, S. L. Lightman, P . D. Taylor, L. Poston,\nand K. T. O’Byrne, “Role of amygdala kisspeptin in pubertal timing in female rats”,PLOS ONE 12, e0183596 (2017).\n12G. Lass et al., “Optogenetic stimulation of kisspeptin neurones within the posterodorsal medial amygdala increases\nluteinising hormone pulse frequency in female mice”, Journal of Neuroendocrinology 32, 10.1111/jne.12823 (2020).\n13G. Lass et al., “GnRH pulse generator frequency is modulated by kisspeptin and GABA-glutamate interactions in the\nposterodorsal medial amygdala in female mice”, Journal of Neuroendocrinology 34, 10.1111/jne.13207 (2022).\n14A. M. Moore, L. M. Coolen, and M. N. Lehman, “Kisspeptin/Neurokinin B/Dynorphin (KNDy) cells as integrators of\ndiverse internal and external cues: evidence from viral-based monosynaptic tract-tracing in mice”, Scientific Reports 9,\n10.1038/s41598-019-51201-0 (2019).\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n13ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\n15S.-H. Yeo, V . Kyle, C. Blouet, S. Jones, and W. H. Colledge, “Mapping neuronal inputs to Kiss1 neurons in the arcuate\nnucleus of the mouse”, PLOS ONE 14, e0213927 (2019).\n16C. McIntyre, X. F. Li, R. de Burgh, D. Ivanova, G. Lass, and K. T. O’Byrne, “GABA Signaling in the Posterodorsal Medial\nAmygdala Mediates Stress-induced Suppression of LH Pulsatility in Female Mice”, Endocrinology 164, 10 . 1210 /\nendocr/bqac197 (2022).\n17K. B. J. Franklin and G. Paxinos, The Mouse Brain in Stereotaxic Coordinates (Boston, Amsterdam, 2004).\n18F. J. Steyn, Y. Wan, J. Clarkson, J. D. Veldhuis, A. E. Herbison, and C. Chen, “Development of a Methodology for and\nAssessment of Pulsatile Luteinizing Hormone Secretion in Juvenile and Adult Male Mice”,Endocrinology 154, 4939–4945\n(2013).\n19R. B. McCosh, M. J. Kreisman, and K. M. Breen, “Frequent Tail-tip Blood Sampling in Mice for the Assessment of Pulsatile\nLuteinizing Hormone Secretion”, Journal of Visualized Experiments, 10.3791/57894 (2018).\n20A. Vidal, Q. Zhang, C. Médigue, S. Fabre, and F. Clément, “DynPeak: An Algorithm for Pulse Detection and Frequency\nAnalysis in Hormonal Time Series”, PLOS ONE 7, e39001 (2012).\n21J. A. Yang, J. K. Hughes, R. A. Parra, K. M. Volk, and A. S. Kauffman, “Stress rapidly suppresses in vivo LH pulses and\nincreases activation of RFRP-3 neurons in male mice”, Journal of Endocrinology 239, 339–350 (2018).\n22G. B. Choi, H.-W. Dong, A. J. Murphy, D. M. Valenzuela, G. D. Yancopoulos, L. W. Swanson, and D. J. Anderson, “Lhx6\nDelineates a Pathway Mediating Innate Reproductive Behaviors from the Amygdala to the Hypothalamus”, Neuron 46,\n647–660 (2005).\n23C. S. Johnson, W. Hong, and P . E. Micevych, “Posterodorsal Medial Amygdala Regulation of Female Social Behavior:\nGABA versus Glutamate Projections”, The Journal of Neuroscience 41, 8790–8800 (2021).\n24H. R. Wilson and J. D. Cowan, “Excitatory and inhibitory interactions in localized populations of model neurons”,\nBiophys J 12, 1–24 (1972).\n25H. R. Wilson and J. D. Cowan, “A mathematical theory of the functional dynamics of cortical and thalamic nervous\ntissue”, Kybernetik 13, 55–80 (1973).\n26S. Keshavarzi, R. K. P . Sullivan, D. J. Ianno, and P . Sah, “Functional Properties and Projections of Neurons in the Medial\nAmygdala”, The Journal of Neuroscience 34, 8699–8715 (2014).\n27X. Bian, “Physiological and morphological characterization of GABAergic neurons in the medial amygdala”, Brain Res\n1509, 8–19 (2013).\n28K. J. Friston, Statistical Parametric Mapping : the Analysis of Functional Brain Images (Elsevier / Academic Press, Amsterdam\n; Boston, 2007).\n29H. R. Wilson and J. D. Cowan, “Evolution of the Wilson–Cowan equations”, Biological Cybernetics 115, 643–653 (2021).\n30R. M. Borisyuk and A. B. Kirillov, “Bifurcation analysis of a neural network model”, Biological Cybernetics 66, 319–325\n(1992).\n31Y. Li, A. Mathis, B. F. Grewe, J. A. Osterhout, B. Ahanonu, M. J. Schnitzer, V . N. Murthy, and C. Dulac, “Neuronal\nRepresentation of Social Information in the Medial Amygdala of Awake Behaving Mice”, Cell 171, 1176–1190.e17 (2017).\n32J. Qiu, C. C. Nestor, C. Zhang, S. L. Padilla, R. D. Palmiter, M. J. Kelly, and O. K. Rønnekleiv, “High-frequency stimulation-\ninduced peptide release synchronizes arcuate kisspeptin neurons and excites GnRH neurons”,Elife 5, 10.7554/eLife.\n16246 (2016).\n33K. Yasuda, K. Raynor, H. Kong, C. D. Breder, J. Takeda, T. Reisine, and G. I. Bell, “Cloning and functional comparison of\nkappa and delta opioid receptors from mouse brain.”, Proceedings of the National Academy of Sciences 90, 6736–6740\n(1993).\n34G. R. Seabrook, B. J. Bowery, and R. G. Hill, “Pharmacology of tachykinin receptors on neurones in the ventral tegmental\narea of rat brain slices”, Eur J Pharmacol 273, 113–9 (1995).\n35A. Dutton and R. E. Dyball, “Phasic firing enhances vasopressin release from the rat neurohypophysis”, J Physiol 290,\n433–40 (1979).\n36E. Doedel, T. Fairgrieve, B. Sandstede, A. Champneys, Y. Kuznetsov, and X. Wang,AUTO-07P: Continuation and bifurcation\nsoftware for ordinary differential equations, 2007.\n37J. Kim, S. J. Semaan, D. K. Clifton, R. A. Steiner, S. Dhamija, and A. S. Kauffman, “Regulation of Kiss1 expression by sex\nsteroids in the amygdala of the rat and mouse”, Endocrinology 152, 2020–30 (2011).\n38K. Kohari, T. Saito, and H. Kawakami, “On a relaxation oscillator with periodic threshold”, in 1993 IEEE International\nSymposium on Circuits and Systems (ISCAS) (1993), 2490–2493 vol.4.\n39B. Christiansen, P . Alstrm, and M. T. Levinsen, “Routes to chaos and complete phase locking in modulated relaxation\noscillators”, Phys. Rev. A 42, 1891–1900 (1990).\n40D. Ivanova, X. F. Li, C. McIntyre, and K. T. O’Byrne, “Posterodorsal Medial Amygdala Urocortin-3, GABA, and Glutamate\nMediate Suppression of LH Pulsatility in Female Mice”, Endocrinology 164, 10.1210/endocr/bqac196 (2022).\n41W. Singer, “Neuronal synchrony: a versatile code for the definition of relations?”,Neuron 24, 49–65, 111–25 (1999).\n42X. F. Li, M. H. Hu, B. P . Hanley, Y. S. Lin, L. Poston, S. L. Lightman, and K. T. O’Byrne, “The Posterodorsal Medial\nAmygdala Regulates the Timing of Puberty Onset in Female Rats”, Endocrinology 156, 3725–3736 (2015).\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n14ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\n43Y. 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\nAmygdala in Stress-Induced Suppression of Pulsatile LH Secretion in Female Rats”, Endocrinology 152, 545–555 (2011).\n44A. N. Comninos et al., “Kisspeptin signaling in the amygdala modulates reproductive hormone secretion”, Brain\nStructure and Function 221, 2035–2047 (2016).\n45S. B. Wolff and B. P . Ölveczky, “The promise and perils of causal circuit manipulations”,Curr Opin Neurobiol 49, 84–94\n(2018).\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n15ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\nFigure\n1. (a) Schematic diagram of the GABA-glutamate circuit in the MePD. Kisspeptin activates glutamatergic and GABAergic neuronal populations\nin the MePD. Both glutamatergic and GABAergic populations interact functionally with a population of GABA efferent neurons providing excitation and\ninhibition, respectively. The activity within the GABAergic efferent neurons and glutamatergic neurons extends to the ARC to modulate the activity of the\nKNDy population. (b) Network representation of the MePD circuit model depicting interactions between the populations of glutamatergic neurons( Gl),\nGABA interneurons (Gi) and GABAergic efferents (Ge). Both Gl and Gi receive kisspeptin excitation (Kp), that is distributed proportionally with relative\nglutamatergic excitation ratio α. The parameters c represent the strength of functional interactions between the populations. Parameters β1 and β2 are\nthe interaction suppression coefficients, used to mimic the decrease in the GABAergic and glutamatergic interaction strengths, respectively, in order to\nsimulate the effects of chemical antagonism.\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n16ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\n0 2 6 10 14\nMePD excitation\n0\n0.1\n0.2\n0.3\n0.4\n0.5Glu population activity\na\n HB\n0 2 6 10 14\nMePD excitation\n0\n0.1\n0.2\n0.3\n0.4\n0.5\nGABAint population activity\nb\n HB\n0 2 6 10 14\nMePD excitation\n0\n0.1\n0.2\n0.3\n0.4\n0.5\nGABAeff population activity\nc\n HB\n0 100 200 300\nTime (s)\n0\n0.1\n0.2\n0.3\n0.4Averaged Activity\nd\n2 4 6 8 10 12 14\nMePD excitation\n20\n40\n60\n80\n100\n120        Period (s)\nf\n0 100 200 300\nTime (s)\n0\n0.1\n0.2\n0.3\n0.4Averaged Activity\ne\nGlu\nGABAeff\n0 2 4 6 8 10 12 14 16 18 20\nMePD excitation\n0.7\n0.8\n0.9\n1\nRelative Glu excitation ratio\ng\nBT\nHB\nSN\nSN\nHC\n-0.2\n-0.1\n0\n0.1\n0.2\nMean MePD output\nFigure\n2. MePD GABA-glutamate circuit dynamics. (a-c) One-parameter bifurcation diagrams of the model for varying kisspeptin (K p) excitation of\nthe GABA-glutamate circuit (MePD excitation). The qualitative dynamics in the populations of (a) glutamatergic neurons (Glu), (b) GABA interneurons\n(GABAint) and (c) GABAergic efferent neurons (GABAeff) changes as the kisspeptin excitation level provided to the circuit varies. Circular markers denote\na 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\nbranch. (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),\ncorresponding 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-\nparameter bifurcation diagram using the level of MePD excitation (K p) and the relative Glu excitation ratio (α) as free parameters. Superimposed heat\nmap indicates the mean MePD output in the oscillatory region. Red, green, and purple lines represent Hopf (HB), saddle-node (SN) and homoclinic\nbifurcations (HC), respectively. Circular markers depict a Bogdanov-Takens (BT) point.\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n17ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\nFigure\n3. Effect of functional interaction strength between the population of GABA interneurons (GABAint) and the population of GABA efferent neurons\nGABAeff on the MePD projections’ output. Heat map depicting (a) the magnitude of mean glutamatergic projections’ output, (b) the magnitude of the mean\nGABAergic output under varying MePD excitation (Kp) and strength of the interaction between the population of GABA interneurons and the population\nof GABA efferent neurons (c ie). For low strength of interaction, GABAergic output increases as MePD excitation is increased, while excessively high\ninteraction drives the activity of the population of GABA efferent (projection) neurons to a zero level. (c) Heat map depicting the combined mean output\nof glutamate and GABA projections. The dashed line represents the selected functional interaction strengths used in the model (see Table 1).\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n18ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\na b\nOscillatory regime Oscillatory regime\nFigure\n4. Effect of varying the MePD functional network connectivity on the MePD projections’ output dynamics. (a) Bifurcation analysis under varying\nfunctional interaction strength between the populations of GABAergic interneurons and glutamatergic neurons (cil and cli). (b) Bifurcation analysis under\nvarying strength of glutamatergic self-excitation (c ll) and functional interaction strength between glutamatergic neurons and GABAergic interneurons\n(cli). The network schematic diagrams highlight the connections under investigation. In the one-parameter bifurcation diagram, the red lines indicate the\nmax of the limit cycle solutions branch, and circle markers denote Hopf ( HB), saddle node (SN), and homoclinic (HC) bifurcation points, respectively.\nThese points are used for the two-parameter continuation. In the two-parameter bifurcation diagrams, red, green, and purple lines represent Hopf (HB),\nsaddle-node (SN) and homoclinic bifurcation curves (HC), respectively. The circular markers depict the Bogdanov-Takens (BT), saddle node loop (SNL)\nand cusp (CP) points. The heat map represents the mean MePD projections’ output in the oscillatory region. The star marker identifies the parameter\nvalues in Table 1.\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n19ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\nb\nedc\nhgf\nFigure\n5. Differential effects of MePD projections on the dynamics of the GnRH pulse generator. (a) Providing oscillatory input to the KNDy network\nmodel leads to quasi-periodic behaviour evolving on a torus; we use a colour scheme to illustrate time evolution on this torus. ( b) Increasing excitation\nin 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\nbifurcation diagram (MePD excitation (K p) and GABAergic interaction suppression coefficient (β 1)) showing the region where the MePD exhibits\noscillatory dynamics and the mean MePD output superimposed on it. (d) Suppression of GABAergic interaction evenly increases both GABAergic and\nglutamatergic tone, leaving the overall MePD output as before the suppression, hence KNDy IPI remains the same. (e) Conjunction of the suppression\nof GABAergic interaction strength with increased MePD stimulation increases the proportion of GABAergic projections’ output, resulting in an increase in\nKNDy IPI. (f) Two-parameter bifurcation diagram (MePD excitation (Kp) and glutamatergic interaction suppression coefficient (β 2)) showing the region\nwhere the MePD exhibits oscillatory dynamics with superimposed mean MePD projections’ output. (g) Suppression of glutamatergic interaction strength\nleads to a loss of oscillatory dynamics in the MePD. Inhibitory GABAergic output is still higher than the excitatory glutamatergic output, leading to the\nreturn of periodicity in the KNDy system and leaving KNDy IPI as before the suppression. (h) Combining the suppression of glutamatergic interaction\nstrength with the increased excitation of the MePD overstimulates ARC KNDy network, resulting in a transition to quiescent dynamics. A more detailed\nversions of the two-parameter bifurcation diagrams depicted in (c) and (f) are provided in Figure S2.\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n20ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\n0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2\nStimulation\n10\n15\n20\n25Average KNDy IPI (min) \na\nTR  \nGlu projection\nGABA projection\n0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2\nStimulation\n0\n0.1\n0.2\n0.3Weighted MePD activity(Hz)\nc\n  HB \nGlu\nGABAeff\nCtrl Ctrl Virus 2 Hz 5 Hz 10 Hz 20 Hz\n10\n15\n20LH IPI (min)\nb\n##\nPrestimulation\nStimulation\nNo stim ctrl\n0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2\nStimulation\n0\n0.02\n0.04\n0.06Mean MePD output (Hz)\nd\nPulses cease\nFigure\n6. In silico simulation of MePD projection stimulation and non-monotonic LH interpulse interval response to stimulation of MePD glutamatergic\nprojection in vivo. (a) Simulation of GnRH pulse generator response to stimulation of GABA and glutamate MePD projections. Mimicking the effects of\nGABA projection stimulation leads to the loss of oscillatory dynamics in the KNDy activity as stimulation level is increased. During glutamate projection\nstimulation 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\nbifurcation (TR) denotes the transition from limit cycle to limit torus in the extended MePD-KNDy network model. The location of TR and HB align.\n(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),\n5 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\ndata points for each animal as circles in the histogram plots. Stimulation at 5 Hz results in a statistically significant suppression of LH pulse frequency\n##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,\nshowing glutamatergic and GABAergic contributions as a function of the level of stimulation (in a.u.). Hopf point (HB) point denotes the location of the\ntransition from stationary to oscillatory dynamics in the MePD circuit. (d) The variation in the weighted mean MePD output under different levels of\nstimulation.\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n21ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\n0 2 6 10 14 18 22\nMePD excitation\n0\n0.05\n0.1\n0.15\n0.2\n0.25\n0.3\n0.35\n0.4\n0.45\n0.5\nGlu population activity\na\n HB\n HC\n SN\n SN\n0 2 6 10 14 18 22\nMePD excitation\n0\n0.05\n0.1\n0.15\n0.2\n0.25\n0.3\n0.35\n0.4\n0.45\n0.5\nGABAeff population activity\nb\n HB\n HC \n SN\n SN\n0 5 10 15 20\nMePD excitation\n5\n10\n15       log(Period (s))\nc\nFigure\nS1. Extended bifurcation diagrams of kisspeptin-mediated excitatiry input to GABA-glutamate circuit fromFigure 2 for (a) Gl and (b) Ge dynamics.\nCircle markers denote Hopf (HB), homoclinic (HC) and saddle node (SN) bifurcations. The red line indicates the maximum amplitude of the limit cycle\nsolutions branch. (c) Excitatory input vs. log Period of oscillations. The period, firstly, decreases, followed by a rapid increase, signifying homoclinic\nbifurcation.\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint \n\n22ro\nyalsocietypublishing.org/journal/rsif J R Soc Interface 0000000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\n0 2 4 6\nMePD excitation\n0\n0.2\n0.4\n0.6\n0.8\n1\nGABA suppression\na\nCP\nCP \nBT SNL \nSNL\nSNL \nSN\nHB\nHC\nSNIC\n0 1 2 3 4 5 6 7\nMePD excitation\n0\n0.2\n0.4\n0.6\n0.8\n1\nGlutamate suppression\nb\nHB\nFigure\nS2. Two parameter bifurcation diagrams, used to identify the oscillatory regions in Figure 5(c,f). (a) Two-parameter bifurcation diagram of\nMePD excitation (K p) and GABAergic interaction suppression coefficient (β 1). (b) Two-parameter bifurcation diagram of MePD excitation (K p) and\nglutamtergic interaction suppression coefficient (β2). Red, green, purple lines represent Hopf (HB), saddle-node (SN) and homoclinic bifurcation curves\n(HC), respectively. The circular markers depict Bogdanov-Takens (BT), saddle-node loop (SNL) and cusp (CP) points.\n.CC-BY 4.0 International licensemade available under a \n(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 \nThe copyright holder for this preprintthis version posted March 4, 2024. ; https://doi.org/10.1101/2024.01.21.574304doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}