A computational model for gonadotropin releasing cells in the teleost fish medaka

Pituitary endocrine cells fire action potentials (APs) to regulate their cytosolic Ca2+ concentration and hormone secretion rate. Depending on animal species, cell type, and biological conditions, pituitary APs are generated either by TTX-sensitive Na+ currents (INa), high-voltage activated Ca2+ currents (ICa), or by a combination of the two. Previous computational models of pituitary cells have mainly been based on data from rats, where INa is largely inactivated at the resting potential, and spontaneous APs are predominantly mediated by ICa. Unlike in rats, spontaneous INa-mediated APs are consistently seen in pituitary cells of several other animal species, including several species of fish. In the current work we develop a computational model of gonadotropin releasing cells in the teleost fish medaka (Oryzias latipes). The model stands out from previous modeling efforts by being (1) the first model of a pituitary cell in teleosts, (2) the first pituitary cell model that fires sponateous APs that are predominantly mediated by INa, and (3) the first pituitary cell model where the kinetics of the depolarizing currents, INa and ICa, are directly fitted to voltage-clamp data. We explore the firing properties of the model, and compare it to the properties of previous models that fire ICa-based APs. We put a particular focus on how the big conductance K+ current (IBK) modulates the AP shape. Interestingly, we find that IBK can prolong AP duration in models that fire ICa-based APs, while it consistently shortens the duration of the predominantly INa-mediated APs in the medaka gonadotroph model. Although the model is constrained to experimental data from gonadotroph cells in medaka, it may likely provide insights also into other pituitary cell types that fire INa-mediated APs.

Introduction Several types of excitable cells elicit electrical pulses called action potentials (APs), which, depending on cell type, can trigger neurotransmitter release, cellular contraction, hormone release or other actions. APs are generated by a combination of ion channels in the plasma membrane, which are typically characterized by the type of ions they are permeable to, and their voltage and/or Ca 2+ dependent gating properties. The primary role of APs in endocrine pituitary cells is to regulate the cytosolic Ca 2+ concentration, which in turn controls the hormone secretion rate in these cells [1]. Hormone secretion often occurs as a response to input from the hypothalamus, peripheral endocrine glands, or from other pituitary cells. However, many endocrine cells are also spontaneously active [1][2][3][4][5][6][7][8][9][10]. The spontaneous activity is partly a means to regulate the re-filling of intracellular Ca 2+ stores, but in several cells also leads to a basal release of hormones. An understanding of the mechanisms regulating the electrodynamics of these cells is therefore fundamental for understanding their overall functioning.
While neuronal APs are predominantly mediated by TTX-sensitive Na + currents (I Na ), AP generation in endocrine cells depends strongly on high-voltage-activated Ca 2+ currents (I Ca ), which in addition to their role in affecting the voltage dynamics of the cell, also are the main source of Ca 2+ entry through the plasma membrane [3,11,12]. In some studies of endocrine cells, APs were exclusively mediated by I Ca , and the spontaneous membrane excitability was insensitive or nearly so to TTX [1,2,[13][14][15][16]. In other studies, APs were evoked by a combination of I Ca and I Na [4,7,[17][18][19]. In one of these studies, TTX was found to block single, brief action potentials, while action potentials of long duration and low amplitude persisted [18], indicating the roles and different time courses of the I Ca and I Na components. The strong involvement of I Ca could explain why pituitary APs typically last longer (typically some tens of milliseconds [8]) than neuronal APs (a few milliseconds), which are mainly mediated by I Na .
All endocrine cells express I Na [8], and TTX sensitive APs can typically be triggered by current injections from hyperpolarized holding potentials even in cells where they are not elicited spontaneously [4,17,20,21]. The reason why the spontaneous activity in many cases is TTX insensitive is likely that a major fraction of I Na is inactivated at the resting membrane potential [15,16]. The reason why this is not always the case, may be that the resting potentials vary greatly between different studies. Only for rat somatotrophs, resting potentials ranging as wide as from −30 mV [13] to −80 mV [18] have been reported.
Computational models constructed to capture the essential activity of pituitary cells have predominantly relied on rat data. The typical resting potentials for rat pituitary cells lie in the range between −50 mV and −60 mV, and at these resting levels, I Na tends to be inactivated and the spontaneous activity TTX insensitive (see reviews in [8,22]). Models based on rat data have therefore typically excluded I Na [3,9,[23][24][25][26][27][28][29][30]. As TTX-sensitive spontaneous APs have been seen in mice corticotrophs [19], I Na was included in a recent computational model of these cells [31], and in a more generic murine pituitary cell model based on the previous rat and mice models [32]. However, the role of I Na in these models was mainly in modulating the firing patterns, and it was not essential for AP firing as such [31,32]. Furthermore, I Ca and I Na were in these models described by simplified kinetics schemes that were adjusted to give the models the desired firing properties, but not explicitly adapted to voltage-clamp recordings of the respective currents in the cell species being modeled.
There are reasons to believe that the dynamical properties of the above cited pituitary cell models are not well suited to represent teleost pituitary cells. Firstly, TTX-sensitive spontaneous activity has been seen in goldfish gonadotrophs resting at −60 mV [4], and TTX sensitive APs has been evoked from a holding potential as high as −50 mV in pituitary cells in cod [7], suggesting that I Na may be more available in resting pituitary cells in fish [4]. Secondly, data from gonadotrophs and somatotrophs in goldfish [4,20] and unspecified pituitary cells in tilapia [5] show APs with very short duration (< 10 ms) compared to the APs in the previous rat models (several tens of ms), putatively indicative of a stronger involvement of I Na . A third difference between fish and rat pituitary cells may be in the role of the big conductance K + current (I BK ), which has been shown to have a paradoxical role in some rat pituitary cells, i.e., it can prolong the duration of I Ca -mediated APs, and sometimes give rise to pseudo-plateau bursts, contrary to what one would expect from a hyperpolarizing current [9,25]. I BK is almost absent in rat gonadotrophs [25], and this was proposed as an explanation to why these cells tend to be less bursty than other pituitary cell types in rats [1,9,32]. In contrast, I BK is highly expressed in medaka gonadotrophs, but without making these cells bursty [12]. The indication that there are differences between rat and fish pituitary cells are further supported by experiments presented in the current work, performed on luteinizing hormone (LH)-producing gonadotroph cells in medaka. We show that these cells elicit brief spontaneous APs that, unlike spontaneous APs in the previous murine pituitary cell models, predominantly are mediated by TTX sensitive Na + currents (I Na ). Furthermore, we show that I BK acts to make APs narrower in medaka gonadotrophs, and thus have the opposite effect from what they have in rat gonadotrophs.
As previous computational models based on murine data seem unsuited to describe the spontaneous activity of teleost pituitary cells, we here present a novel pituitary cell model constrained to data from medaka gonadotrophs. Given the putatively complex interplay between I Ca and I Na during the AP upstroke, we put extra effort into developing accurate models of these two currents, and constrained their kinetics to voltage-clamp recordings of the individual currents. In addition to I Ca and I Na , we included a leak current and three K + currents in the model. These we adopted from previous pituitary cell models, and adjusted to adapt the firing properties of the model to current-clamp recordings from medaka gonadotrophs under control conditions, after application of TTX, and after application of the I BK blocker paxilline.
By comparing the medaka gonadotroph model, which predominantly fires I Na -mediated APs, with three previous models of murine pituitary cells, which predominantly [32] or exclusively [9,27] fire I Ca -mediated APs, we explore the consequences of having different AP-generating mechanisms. We find that the medaka gonadotroph model produces spontaneous APs that are faster than those in the murine models, and thus more suited to describe the firing properties of teleost pituitary cells. Furthermore, we show that while I BK may broaden APs in the murine pituitary models [9,27,32], it consistently had a narrowing effect on APs in the medaka gonadotroph model, and propose explanations to these model differences. By this, we add to the discussion of the role played by I BK in shaping pituitary APs [8], and suggest that the effect of I BK on APs is mainly determined by the timing of I BK -activation relative to the AP peak, as also proposed previously [9].
Although the model presented here was tailored to represent gonadotroph cells in medaka, we believe that it is of a more general value for improving our understanding of I Na -based APs in the pituitary, which are elicited by several endocrine cell types and in several animal species, depending on biological conditions [4,7,17,18,[33][34][35].

Characteristic response patterns of medaka gonadotrophs
The general electrophysiological properties of gonadotroph cells in medaka were assessed through a series of voltage-clamp and current-clamp experiments. The voltage-clamp experiments used to develop kinetics models of Na + and Ca 2+ currents are presented in the Methods section. Here, we focus on the key properties of spontaneous APs as recorded in current clamp. Selected, representative experiments are shown in Fig 1. Although variations were observed, the medaka gonadotrophs typically had a resting potential around −50 mV, which is within the range found previously for goldfish [4] and cod [10] gonadotrophs. As for goldfish gonadotrophs, the majority of medaka gonadotrophs fired spontaneous APs with peak voltages slightly below 0 mV. The spontaneous APs were always regular spikes (i.e., not bursts or plateau-like) and the AP width varied between 3 and 7 ms (blue traces in Fig 1). Similar brief AP waveforms have been seen in previous studies on fish [4,5,20], while the APs reported for rat gonadotrophs are typically slower, i.e. from 10-100 ms [8]. The spontaneous AP activity was completely abolished by TTX application (Fig 1B).
Finally, we explored how paxilline (an I BK blocker) affected the spontaneous activity of medaka gonadotrophs. In the experiment shown in Fig 1D, paxilline increased the firing rate and slightly reduced the mean AP peak amplitude, but these effects were not seen consistently in experiments using paxilline application. However, in all experiments, paxilline application was followed by a small increase of the resting membrane potential (Fig 1C), and a broadening of the AP waveform (Fig 1D2 and 1D3). Similar effects have been seen in goldfish somatotrophs, where application of tetraethylammonium (a general blocker of K + currents) lead to broadening of APs [20]. The effect of I BK in goldfish and medaka gonadotrophs is thus to make APs narrower, which is the opposite of what was found in rat pituitary somatotrophs and lactotrophs, where I BK lead to broader APs and sometimes to burst-like activity [9,25].

A computational model of medaka gonadotrophs
The model for medaka gonadotrophs is described in detail in the Methods section, but a brief overview is given here. The dynamics of the membrane potential is determined by the differential equation: The three K + currents, I K , I BK and I SK , were based on previously published models ( [9,32]), but adjusted (see Methods) so that the final model had an AP shape and AP firing rate that were in better agreement with the experimental data in Fig 1. I K denotes the delayed rectifier K + channel [9], I SK denotes the small-conductance K + channel, activated by the intracellular Ca 2+ concentration [9], and I BK denotes the big-conductance K + channel. The latter was assumed to be both voltage and Ca 2+ -dependent. As it is often co-localized high-voltage-activated Ca 2+ channels, it was assumed to sense a domain-Ca 2+ concentration proportional to I Ca [32]. The depolarizing membrane currents consisted of a high-voltage activated Ca 2+ current (I Ca ) and a Na + current (I Na ), both of which were novel for this model, and adapted to new voltage-clamp data from gonadotroph cells in medaka (see Methods). I Na activated in the range between −50 mV and −10 mV, with half activation at −32 mV (black curves, Fig 2A1), quite similar to what was previously found in goldfish gonadotrophs [4]. I Na inactivated in the range between −90 mV and −40 mV, with half-inactivation at −64 mV, which was lower than (C1) Spontaneous activity before (blue) and after (red) TTX application. (C2-C3) Close-ups of two selected events before (blue) and after (red) TTX application. (D1) Spontaneous activity before (blue) and after (red) paxilline application. (D2-D3) Close-ups of two selected events before (blue) and after (red) paxilline application. The firing rates in the various recordings were 0.64 Hz (A1), 0.57 Hz (B1), 1.22 Hz (C1, before TTX), 0.17 Hz (D1, before paxilline) and about 0.35 Hz (D1, after paxilline). AP width (defined as the time between the upstroke and downstroke crossings of the voltage midways between −50 mV and the peak potential) varied between 3 and 7 ms, with mean width of 3.7 ms (A1), 4.9 ms (B1), 3.7 ms (C1, before TTX). In (D1), average AP widths were 4.2 ms before paxilline. After paxilline, the AP shapes varied, with AP widths ranging from 9 ms to 90 ms, with a mean of 25 ms. AP peak voltages varied between -11.8 mV and +5.5 mV, with mean peak values of −0.4 mV (A1), −3.4 mV (B1), −5.1 mV (C1, before TTX), −3.1 mV (D1, before paxilline), and −6.4 mV (D1, after paxilline). AP width was calculated at half max amplitude between −50 mV and AP peak. The experiments were performed on gonadotroph LH-producing cells in medaka. All depicted traces were corrected with a liquid junction potential of −9 mV. The time indicated below each panel refers to the duration of the entire trace shown. in goldfish, where the half-inactivation was found to be around −50 mV [4]. With the activation kinetics adapted to medaka data, only 6% of I Na was available at the typical resting potential of −50 mV. The fact that medaka still showed TTX-sensitive spontaneous activity thus suggests that I Na is highly expressed in these cells. In comparison, in the generic murine pituitary model [32], I Na activation required depolarization to voltages far above the resting potential (red curves, Fig 2A1)), meaning that this model could not elicit spontaneous I Na -based APs.
Both I Na and I Ca had fast activation in medaka gonadotrophs, I Na being slightly faster with a time constant of about 0.5-0.8 ms in the critical voltage range (Fig 2A2), whereas I Ca had a time constant >1 ms in the critical voltage range (Fig 2B2). I Ca activated in the range between −40 mV and +10 mV, with a half activation at 16 mV (red curve in Fig 2B1). This activation curve was much steeper than I Ca in the rat models (colored curves). The high activation in medaka gonadotrophs threshold suggests that I Ca is unsuitable for initiating spontaneous APs, making spontaneous activity critically dependent on I Na , unlike in all rat models, where I Ca was highly available around the resting potential.

BK currents cause briefer APs in medaka gonadotrophs
With the right tuning, the medaka gonadotroph model reproduced the essential firing patterns seen in the experiments (Fig 1). In control conditions, it fired sharp APs (AP width was 6 ms) with relatively low peak voltages (around -10 mV), and had a spontaneous firing rate slightly below 1 Hz (Fig 3A). AP firing was completely abolished when the Na + -conductance g Na was set to zero, mimicking the effect of TTX application in the experiments.  A series of previous models have shown that I BK may act to broaden APs and promote bursting in rat pituitary cells [9,25,27,32]. In contrast, a high I BK expression in medaka gonadotrophs [12] does not make these cells bursty. On the contrary, the experiments in Fig 1D showed that medaka gonadotrophs fired broader APs when BK channels were blocked. This was also seen in the model simulations, when the BK-channel conductance (g BK ) was reduced relative to its value during control (Fig 3B-3F). Generally, a reduction in g BK lead to a broadening of the AP event. The resemblance with paxilline data was strongest in simulations with partial reduction, such as in Fig 3C and 3D, where g BK had been reduced to 16% and 15% of its control value, respectively. Then, AP events were about 60 ms wide, and included plateau potentials that presumably reflected an interplay between I Ca and repolarizing currents activating/inactivating after the initial AP peak. When g BK was set to zero, the oscillations were not seen, and the AP was prolonged by a flatter and more enduring > 100 ms plateau. It is reasonable to assume that also in the experiments, the blockage of BK by paxilline was not complete.

Membrane mechanisms controlling the AP width
To explore in further detail how the various membrane mechanisms affected the AP firing, we performed a feature-based sensitivity analysis of the medaka gonadotroph model (Fig 4). We then assigned the maximum conductances of all included currents uniform distributions within intervals ±50% of their default values (Table 1), and quantified the effect that this parameter variability had on selected aspects of the model output (see Methods). An exception was made for g BK , which was assigned a uniform distribution between 0 and the maximum value given in Table 1), i.e., from fully available to fully blocked, motivated by the fact that this was the possible range spanned in the paxilline experiments ( Fig 1D). We note the total-order Sobol sensitivity indices considered in the current analysis reflects complex interactions between several nonlinear mechanisms, and that mechanistic interpretations Computational model for medaka gonadotrophs therefore are difficult. Below, we have still attempted to extract the main picture that emerged from the analysis.
Three features of the model responses were considered: (i) IsBursting, (ii) IsRegular, and (iii) IsNotSpiking. Following the definition used by Tabak et al. [9], plateau potentials of duration longer than 60 ms (such as those in Fig 3C2-3F2) were defined as bursts. For simplicity, we used the definition loosely, and referred to enduring plateau potentials as bursts even in cases where they did not contain any oscillations (such as in Fig 3C2-3F2)). APs of shorter duration than this (such as those in Fig 3A2 and 3B2) were defined as regular spikes. All the features (i-iii) were binary, meaning e.g., that IsBursting was equal to 1 in a given simulation if it contained one or more bursts, and equal to 0 if not. The mean value of a feature (taken over all simulations) then represented the fraction of simulations that possessed this feature. For example, IsBursting had a mean value of 0.11, IsRegular had a mean value of 0.35, and IsNot-Spiking had a mean value of 0.55, which means, respectively, that 11% of the model parameterizations fired bursts, 35% fired regular APs, and 55% did not fire any kind of APs. AP activity was thus seen in less than half of the parameterizations. This reflects that the default configuration had a resting potential only slightly above the AP generation threshold, so that any parameter re-sampling that would make the cell slightly less excitable, would abolish its ability for AP generation. We note that the mean values of the three features sum up to 1.01 and not to unity. This was because a few of the parameterizations fired both bursts and regular APs within the same simulation.
The total-order Sobol indices, shown in histograms in Fig 4, quantify how much of the variability (between different simulations) in the response features that are explained by the variation of the different model parameters, i.e., the maximum conductances. When interpreting these results, we should keep in mind that the feature sensitivities are not independent, i.e., if Isbursting equals 0 for a given implementation, it means that either IsRegular or IsNotSpiking must equal 1. When the sensitivity to g Na was small for IsBursting, but quite large for IsRegular and IsNotspiking, it then means that g Na was important for switching the model between not firing and regular firing, while it contributed less to prolonging the APs into possible bursts. In contrast,IsNotSpiking was almost insensitive to g BK , while IsBursting and IsRegular had a high sensitivity to g BK . A little simplified, we can thus say that g Na determined whether the model fired an AP (Fig 4C), while g BK determined whether the AP, if fired, became a burst or a regular spike (Fig 4A and 4B).
All three features had a high sensitivity to g K , which indicates that g K played multiple roles for the firing properties of the model. Firstly, I K had a nonzero activity level around rest (cf. Fig 2C1), and was important (along with I Na and the leakage current, I l ) for determining whether the resting potential was above the AP threshold, hence the high sensitivity of IsNot-Spiking to g K . Having a broad activation range, I K was also important for repolarizing the Computational model for medaka gonadotrophs membrane potential after the AP peak, and thus for the duration of the AP. Therefore, also IsBursting had a high sensitivity to g K . The mechanisms behind burst generation are reflected in the sensitivity of IsBursting to g BK , g K and g Ca . Here, I Ca is responsible for mediating the plateaus that prolong APs into possible bursts, while g BK and g K may prevent bursts by facilitating a faster down-stroke. The interaction between g BK and g K in mediating the AP downstroke is complex, as we comment on further in the next subsection. The sensitivity to the last K + -channel, g SK , was very low in all the features considered here. g SK had very little impact on the AP shape or the ability of the model to elicit APs, but was included in the model since it was important for regulating the firing rate.

BK currents have opposite effects on AP shape in different cells
As we have seen, I BK consistently had a narrowing effect on APs in the medaka gonadotroph model, while it has previously been reported to broaden APs in several models based on data from murine pituitary cells [9,27,32]. To gain insight into this dual role, we here explore the relationship between g BK expression and AP shape in four different models, including (i) the medaka gonadotroph model (Fig 5A), (ii) a previously published models of a rat lactotroph (Fig 5B), (iii) a generic pituitary cell model based on data from rats and mice (Fig 5C), and (iv) a model of a rat somatotroph (Fig 5D). In the medaka gonadotroph model, APs were predominantly mediated by I Na , while in the murine models, APs were predominantly mediated by I Ca . Despite several differences, all models contained I BK and I K , which were the most important ion channels for mediating the AP downstroke.
In all models, an increase in g BK consistently lead to a reduction of the AP-peak voltage ( Fig  5A1-5D1), as one might expect from a hyperpolarizing current. Additional simulations on the medaka gonadotroph model revealed that the magnitude of the reduction depended strongly on the I BK -activation time constant (τ BK ). In the default configuration, τ BK was set to 3 ms (orange curve in Fig 5A1). With a slower τ BK , I BK activation remained low during the AP upstroke, and its effect on the AP-peak voltage was low (red curve in Fig 5A1). Contrarily, when τ BK was faster, I BK activation largely occurred during the AP upstroke, and I BK then had a larger effect on the AP-peak value (blue curve in Fig 5A1). In general, I BK could affect both the upstroke (reducing the peak voltage) and downstroke (repolarizing the membrane) of the AP, and the relative contribution to the two phases would depend on the relative timing of I BK activation and the AP peak.
The effect on g BK on AP width (Fig 5A2-5D2) and afterhyperpolarization (AHP) depth, defined simply as the minimum voltage reached between two APs (Fig 5A3-5D3), was more complex and less intuitive. To gain an insight into the mechanisms at play, we start by exploring how g BK affected the AHP depth (Fig 5A3-5D3). The AHP was predominantly due to the joint effect of the two hyperpolarizing currents, g BK and g K . It may therefore seem counterintuitive that, for low g BK , an increase in g BK actually decreased the AHP (less hyperpolarization means higher AHP voltages). The explanation lies in the simultaneous effect that g BK had on reducing the AP-peak voltage (Fig 5A1-5D1). Lower AP-peak values generally meant less I K activation [9,25], as this current activates at high voltage levels. Hence, g BK had a dual affect on the AHP. It could facilitate AHP through its direct, hyperpolarizing effect, but at the same time counteract AHP indirectly, by limiting g K activation. When g BK was small, the AHP was predominantly mediated by g K , and the indirect effect dominated, so that an increase in g BK reduced the AHP up to a certain point, where the direct effect to took over, so that a further increase in g BK enhanced the AHP.
The dual role of g BK is also reflected in the effect that g BK had on the AP width (Fig 5A2-5D2). An increase in g BK could either lead to briefer APs, through the direct hyperpolarizing effect of I BK , or broader APs, through the indirect effect of g BK reducing I K activation. This dual role of g BK is seen most clearly in the rat lactotroph model (Fig 5B2). For low values of g BK , the direct effect dominated, and the AP width decayed monotonically with increasing g BK up to a certain threshold value, where a further increase in g BK gave a sharp transition to very broad APs (pseudo-plateau bursts). The paradoxical role of I BK as a burst-promoter in the lactotroph model was explored in detail in the original study [9], and in a later re-implementation of the model [37]. The same dual role of g BK on the AP shape was seen in the generic pituitary . The inset in A2 shows the same curves as the main panel, but with a wider range on the y-axis. (B) The rat lactotroph model was taken from [9]. (C) The generic murine pituitary cell model was taken from [32]. Two versions were considered, one with (orange curve) and one without (blue curve) sodium conductances added in the model. (D) The rat somatotroph model was taken from [27]. (A1-E1) The AP-peak voltage decayed monotonically with g BK in all models. (A2-D2) An increase in g BK could lead to both a broadening and narrowing of APs, depending on conditions. AP width was defined as the time between the upstroke and downstroke crossings of the voltage midways between −50 mV and the peak potential. (A3-D3) An increase in g BK could cause both stronger or weaker afterhyperpolarization (AHP), depending on conditions. AHP was defined as the minimum voltage reached between two spikes. (A-D) In all panels, the x-axis showed g BK relative to a reference value g BKref , which was taken to be the default BK-conductance in the respective models.
In the rat somatotroph model, the indirect effect dominated for all values of g BK , and the AP width increased monotonically with increasing g BK (Fig 5D2). Oppositely, in the default parameterization of the medaka gonadotroph model, the direct effect dominated for all values of g BK , and the AP width decreased monotonically with increasing g BK (orange curve, Fig  5A2). Only by decreasing τ BK to unrealistically low values, g BK could have a broadening effect on APs in the medaka gonadotroph model (blue curve, inset in Fig 5A2).
We note that the complex interplay of mechanisms is only partly captured by the simplified, heuristic explanations presented above. In particular, for high g BK values, neither the AP width or AHP increased monotonically with g BK in all models (Fig 5B2 and 5C3). This non-monotonic behavior putatively reflects a complex and highly sensitive interplay between several mechanisms active in the aftermath of the AP peak, and we did not attempt to explore it in further detail.
As we have seen, I BK facilitated bursting in all the considered models based on murine data, but not in the default parameterization of the medaka gonadotroph model. As the I BK kinetics in the medaka gonadotroph model was essentially the same as in the generic murine pituitary cell model [32], we hypothesized that the different role played by I BK in the medaka gonadotroph model versus the murine pituitary cell models was due to differences in AP shape, rather than I BK kinetics. By comparing the AP upstrokes of the different models, we see that the fastest upstroke was found in the medaka gonadotroph model where APs were predominantly I Na mediated (blue curve in Fig 6). The second fastest upstroke was seen in the generic murine model in the case where its APs were mediated by a combination of I Ca and I Na (purple curve in Fig 6). Hence, the addition of I Na to this model made the AP upstroke steeper, and, as we saw in Fig 5C2, this made the model less susceptible to bursting, i.e., the transition to bursting occurred for a much higher value of g BK .
In the remaining models, APs were mediated solely by I Ca , and had slower upstroke. Hence, in the murine models, I BK had more time to activate during the AP upstrokes, which explains why I BK indirectly could promote bursting in these models, by reducing AP-peak value and thereby I K activation. In order to have the same effect in the medaka gonadotroph model, τ BK needed to be speeded up dramatically, as illustrated in the blue curve in Fig 5A2. This scenario is hypothetical, as no experiments suggest that I BK does promote bursting in medaka gonadotrophs. However, it is interesting to note that this parameterization of the model fired bursts (i.e., APs with width > 60 ms) both for very low and very high values for g BK , while it fired regular AP for intermediate g BK value. Also a previous model of a rat corticotroph showed such bursting behaviour in two disjoint regions in g BK -space (see Fig 3 in [31]).
In summary, I BK had an inhibitory effect on I K by reducing the AP amplitude, and a collaborative effect with I K in mediating the AP downstroke. With slow AP upstrokes, as in the murine pituitary cell models, the inhibitory effect of I BK on I K could dominate, and I BK could result in a net reduction of hyperpolarization and as such promote broader APs and sometimes bursts. With faster AP upstrokes, such as in the medaka gonadotroph models, the collaborative effect of I BK always dominated, and I BK consistently facilitated narrower APs. In a broader scope, this suggests that I BK can act as a mechanism that reduces the duration of already fast APs, and prolongs the duration of already slow APs.

Discussion
TTX-sensitive Na + currents (I Na ) are present in all pituitary cells, but are in many cases inactive during spontaneous activity [8]. Previous models of the electrical activity of pituitary cells have focused on conditions where I Na is of lesser importance, and where AP generation is predominantly mediated by high-voltage activated Ca 2+ currents [3,9,23,[25][26][27][28]. To our knowledge, we have in the current work presented the first models that describe pituitary cells under conditions where AP generation is I Na -mediated. The model was adapted to experimental data from LH-producing gonadotrophs in medaka, whose spontaneous activity is highly I Na -dependent. Voltage-clamp data was used to develop models for the activation kinetics for I Na and I Ca currents, and the firing properties of the model were further adapted to currentclamp data from spontaneously active cells (under control conditions, and after application of TTX and paxilline).
To examine the consequences of having different AP generation mechanisms, we performed a comparison between the the medaka gonadotroph model, which fired I Na -mediated APs, and three models of murine pituitary cells which fired APs that were exclusively mediated Computational model for medaka gonadotrophs by I Ca [9,27], or by a combination of I Ca and I Na [32]. The most interesting result that came out of this comparison was that I BK had a dual role on AP shape, and could under some conditions broaden APs and promote bursting, and under other conditions make them narrower. We suggested that the broadening effect could only occur in scenarios where I BK had sufficient time to activate during the AP upstroke, and thus required either a very fast I BK activation time constant, or a relatively slow AP upstroke. In the murine models, the AP upstrokes were slower than in the medaka gonadotroph model, and we suggest that this explains why increased I BK can promote bursting in many murine pituitary cells [8,9,25,32], but not in medaka gonadotrophs [12]. Also other K + channels have been shown to have such a burst-promoting role in murine pituitary cells [38,39]. It should be noted that the effect on reducing AP width is a commonly reported role for I BK in many excitable cells [40][41][42][43][44], while the AP-broadening and burst promoting effect that I BK is less conventional.
The role of I BK as a burst promoter has not been found consistently in rat lactotrophes. In the study by Miranda et al. 2003, AP width in rat GH 3 , a widely used model for pituitary lactotrophs, was instead found to increase when I BK was blocked with paxilline [43], i.e., similar to what we found for medaka gonadotrophs (Fig 1D). The different effects of I BK on AP width observed in different laboratories [25,40,43] was addressed by Tabak et al. 2011 [9], who proposed possible explanations that could reconcile the conflicting results. One possible explanation could be there is a variability in terms of how BK channels are localized in various cells, and that BK channels that are co-localized with Ca 2+ channels will respond rapidly to voltage fluctuations and promote bursting, while BK channels that are not co-localized with Ca 2+ channels will react more slowly to voltage fluctuations and have the opposite effect [9]. A second possible explanation, also suggested by Tabak et al. 2011, was that I BK might have different kinetic properties in different cells due to variations in their phosphorylation state [9]. A third explanation could be that different cells have different BK splice variants [45], or different regulatory sub-units.
The model comparison in Fig 5 provides additional possible explanations to the conflicting conclusions regarding the role of I BK in lactotrophs. Firstly, the fact that I BK has affects the AP shape differently in different cells does not by necessity reflect differences in I BK kinetics or localization. Simulations on the rat lactotroph model showed that the same I BK could both have a broadening and narrowing effect on the APs within one and the same model ( Fig 5B). That is, APs could be made broader either by reducing g BK to very low values, or by increasing it to very high values. This dual effect of I BK was even more pronounced in a version of the medaka gonadotroph model (blue curve in Fig 5B), and a previous model of rat corticotrophs [31], which elicited bursts both under full g BK blockage and for large g BK expression, while they fired regular APs for intermediate g BK expression. Hence, in general, g BK blockage could affect the AP width in either way, depending on the initial level of g BK expression, and conflicting conclusions regarding the role of I BK could reflect that variations in g BK expression under control conditions. Secondly, the way on which I BK will affect the APs in a given cell can not be predicted from g BK kinetics/expression alone, but also depends on the AP generating mechanisms in the cell. Our simulations suggested that APs with a steep upstroke were prone to be made briefer by I BK , while APs with a slower upstroke were prone to be prolonged by I BK , as also suggested in the dynamic clamp experiments by Tabak et al. [9]. Putatively, I Na mediated APs will generally have a steeper upstroke than I Ca mediated APs, as was the case in the models studied here. If this is the case, the role of I BK in a given cell may thus largely be determined by which membrane mechanisms that mediate its AP upstroke, and especially the degree to which I Na is involved, which is highly resting-potential dependent, and likely to depend strongly on experimental conditions. Differences in I Na involvement could in principle explain the conflicting experiments on rat lactotrophs [25,43]. In the experiment by , where I BK was found to broaden APs, APs were predominantly mediated by I Ca [9,25]. In the experiment by Miranda et al. 2003, where I BK was found to narrow APs (i.e., blocking I BK lead to broader APs), it was reported that this only occurred under conditions in which short APs were present. It is likely that the events described in that work as short APs were largely I Na -mediated, so that the differences between the studies by Van Goor et al. Miranda et al. 2003 suggest that different AP generation mechanisms dominated in the two experiments.
Although the medaka gonadotroph model captured the essential firing properties of medaka gonadotrophs, the agreement between model and data was not perfect. For example, the AP width during control conditions (Fig 3A2) was in the upper range of that seen in the experiments, while AP peak voltage was in the lower range of what was seen in the experiments (Fig 1A and 1B). We were not able to obtain briefer APs with larger peak values without compromising the agreement between the experimental data and other model features, such as afterhyperpolarization, firing rate and response to I BK -blockage. The conductances selected for the default model (Table 1) were thus a compromise made to obtain an acceptable match to several features simultaneously. The fact that we were not able to obtain a more accurate match between model and data likely reflect that some of the ion channels present in the model are imperfect representations of the ion channels present in the real cell. For example, the simplified kinetics schemes used for I K , I BK and I SK were adopted from models of rat pituitary cells [9,32], and were not constrained to data from medaka gonadotrophs. In addition, the biological cell is likely to contain a variety of additional ion channels that were not included in the model.
To our knowledge, the medaka gonadotroph model is the first computational model of an endocrine cell that fires APs that are predominantly I Na -based. Although it was adapted to experimental recordings from LH-producing gonadotrophs in medaka, we believe that the model has a more general value. Different types of pituitary cells in several different species share many of the same membrane mechanisms [8]. In particular, I Na -based APs are elicited by several pituitary endocrine cell types and in several animal species, depending on biological conditions [4,7,17,18,[33][34][35]. It is thus likely that the response patterns of related cell types may be captured by up-or down-regulation of selected mechanisms included in medaka gonadotroph model. Insight in which parameters that should be adjusted in order to obtain desired changes in the model's firing properties could then be obtained through a sensitivity analysis, such as that presented in Fig 4, or in previous, more comprehensive studies that consider a larger number of model features [37,39].

Experimental procedures
The electrophysiological experiments were conducted using the patch-clamp technique on brain-pituitary slices from adult female medaka (as described in [46]). To record spontaneous action potentials and Ca 2+ currents we used amphotericin B perforated patch configuration, while for Na + currents we used whole cell configuration. Extracellular (EC) solution used for recording spontaneous action potentials (current clamp) contained 134 mM NaCl, 2.9 mM KCl, 2 mM CaCl 2 , 1.2 mM MgCl 2 1.2, 10 mM HEPES, 4.5 mM glucose. The solution was adjusted to a pH of 7.75 with NaOH and osmolarity adjusted to 290 mOsm with mannitol before sterile filtration. Before use, the EC solution was added 0.1% bovine serum albumin (BSA). For Na + current recordings (voltage clamp) we used a Ca 2+ free and Na + fixed (140 mM) EC solution, pH adjusted with trizma base. In addition, 10 μM nifedipine, 2 mM 4-Aminopyridine (4-AP) and 4 mM Tetraethylammonium (TEA) was added to the EC solution just before the experiments. To record Ca 2+ currents, we substituted NaCl with 120 mM choline-Cl and added 20 mM Ca 2+ , 2 mM 4-AP and 4 mM TEA. The patch pipettes were made from thick-walled borosilicate glass using a horizontal puller (P 100 from Sutter Instruments). The resistance of the patch pipettes was 4-5 MO for perforated patch recordings and 6-7 MO for whole-cell recordings. For recordings of spontaneous action potentials, the following intracellular (IC) electrode solution was added to the patch pipette: 120 mM KOH, 20 mM KCl, 10 mM HEPES, 20 mM Sucrose, and 0.2 mM EGTA. The pH was adjusted to 7.2 using C 6 H 13 NO 4 S (mes) acid, and the osmolality to 280 mOsm using sucrose. The solution was added 0.24 mg/ml amphotericin B to perforate the cell membrane (see [46] for details). In voltage clamp experiments the K + was removed from the intracellular solution to isolate Na + and Ca 2+ currents. This was achieved by substituting KOH and KCl with 130 mM Cs-mes titrated to pH 7.2 with CsOH. The electrode was coupled to a Multiclamp 700B amplifier (Molecular Devices) and recorded signal was digitized (Digidata 1550 with humsilencer, Molecular Devices) at 10 KHz and filtered at one-third of the sampling rate. In selected experiments, voltage-gated Na + channels were blocked using 5 μM TTX, and BK channels were blocked using 5 μM paxilline. Both drugs were dissolved in EC solution and applied using 20 kPa puff ejection through a 2 MO pipette, 30-40 μm from the target cell.
Under the experimental (voltage-clamp) conditions used for recording Na + currents, and under the experimental current-clamp conditions, a liquid junction potential of about −9 mV was calculated and corrected for in the data shown in Fig 1, and in the kinetics model for I Na (Fig 2A). A liquid junction potential of about −15 mV was calculated for the experimental (voltage-clamp) conditions used for recording Ca 2+ currents, and was corrected for in the kinetics model for I m Ca (Fig 2B).

Ethics statement
Handling, husbandry and use of fish were in accordance with the guidelines and requirements for the care and welfare of research animals of the Norwegian Animal Health Authority and of the Norwegian University of Life Sciences.

Model of medaka gonadotroph
As stated in the Results-section, the medaka gonadotroph model was described by the equation: The membrane capacitance was set to the standard value C m = 1μF/cm 2 , and the leak conductance was described by with a reversal potential E leak = −45 mV. Due to a nonzero-activation of I K around the resting level, this gave an effective resting potential around −50 mV, similar to that in the experimental data in Fig 1. The kinetics of all ion channels were summarized in Fig 2, but are described in further detail here. I Na was modeled using the standard Hodgkin and Huxley-form [47]: with a reversal potential E Na = 50 mV, and gating kinetics defined by: The steady-state activation and time constants (q 1 , h 1 , τ q and τ h ) were fitted to voltage-clamp data from medaka gonadotrophs, as described below, in the subsection titled "Model for the voltage-gated Na + channels". I Ca was modelled using the Goldman-Hodgkin-Katz formalism, which accounts for dynamics effect on Ca 2+ reversal potentials [48]: with Here, R = 8.314J/(mol � K is the gas constant, F = 96485.3C/mol is the Faraday constant, T is the temperature, which was set to 293.15 K in all simulations. [Ca] and [Ca] e were the cytosolic and extracellular Ca 2+ concentrations, respectively. The former was explicitly modelled (see below), while the latter was assumed to be constant at 2 mM. As Eq 6 shows, we used two activation gates m. This is typical for models of L-type Ca 2+ channels (see e.g. [49][50][51][52]), which are the most abundantly expressed HVA channels in the cells studied here [11]. The steady-state activation and time constant (m 1 and τ h ) were fitted to voltage-clamp data from medaka gonadotrophs, as described below, in the subsection titled "Model for high-voltage activated Ca 2+ channels". We note that g Ca in the Goldman-Hodgkin-Katz formalism (Eq 6) is not a conductance, but a permeability with units cm/s. It is proportional to the conductance, and for simplicity, we have referred to it as a conductance in the text. The delayed rectifyer K + channel was modelled as with reversal potential E K = −75 mV, and a time dependent activation gate described by The steady-state activation was described by: with a slope parameter s n = 10 mV, and half-activation v n = −5 mV. The model for I K was identical to that in a previous rat lactotroph model [9], with the exception that the constant (voltage-independent) time constant τ K was made faster (5 ms) in the medaka gonadotroph model to account for the more rapid APs elicited by these cells.
The model for the BK-channel kinetics was a was taken from a previous model (where it was called BK-near) of murine corticotrophs [29,30], which has also been used in a generic rat pituitary cell model [32]: The activation kinetics was given by: The constant (voltage-independent) activation time constant τ BK was set to 3 ms. The steadystate activation was given by [29]: with As BK channels are often co-localized with high-voltage activated Ca 2+ channels, BK activation was assumed to depend on a domain concentration c dom in Ca 2+ nanodomains, which in turn was assumed to be proportional to the instantaneous Ca 2+ influx through I Ca . We therefore set: where A = 1.21 mmol�cm −1 �C −1 is a parameter that converts a current density into a concentration, and c ref = 2μM is a reference concentration. The parameter A was not taken from previous studies, but tuned so that the model obtained suitable BK activation in simulations on the medaka gonadotroph model. Finally, the SK channel was the same as in the previous model of a rat lactotroph [9], and was modelled as: with an instantaneous, Ca 2+ dependent, steady-state activation: where [Ca] denotes the cytosolic Ca 2+ concentration, and k s is a half-activation concentration of 0.4 μM. I Ca and I SK were dependent on the global cytosolic Ca 2+ concentration. This was modelled as a simple extrusion mechanism, receiving a source through I Ca , and with a concentration dependent decay term assumed to capture the effects of various ion pumps and buffering mechanisms: Here, f c = 0.01 is the assumed fraction of free Ca 2+ in the cytoplasm, α = 0.015mM � cm 2 /μC converts an incoming current to a molar concentration, and k c = 0.12ms −1 is the extrusion rate [9]. The conductances used in the default parameterization of the model are given in Table 1.

Model for the voltage-gated Na + channels
The steady-state values and time courses of the gating kinetics were determined using standard procedures (see e.g. [35,47,53,54]), and was based on the experiments summarized in Fig 7. To determine activation, the cell was held at −60 mV for an endured period, and then stepped to different holding potentials between −80 to 100 mV with 5 mV increments (Fig 7A2), each for which the response current (I Na ) was recorded ( Fig 7A1). The inactivation properties of Na + were investigated using stepwise pre-pulses (for 500 ms) between −90 and 55 mV with 5 mV increments before recording the current at −10 mV (Fig 7B2). The resulting Na + current then depended on the original holding potential (Fig 7B1). Finally, the recovery time for the Na + current was explored by exposing the cell to a pair of square pulses (stepping from a holding potential of −60 mV to −10 mV for 10 ms) separated by a time interval Δt (Fig 7C2). The smallest Δt was 10 ms, and after this, Δt was increased with 100 ms in each trial. The cell responded to both pulses by eliciting Na + -current spikes (Fig 7C2). When Δt was small, the peak voltage of the secondary spike was significantly reduced compared to the first spike, and a full recovery required a Δt in the order of 1/2-1 s.

Steady-state activation and inactivation.
To determine steady-state activation, the peak current (I max ) was determined for each holding potential in Fig 7A, and the maximum peak was observed at about −10 mV. For inactivation, the peak current (I max ) was recorded for each holding potential in Fig 7B. In both cases, the maximal conductance (g max ) for each holding potential was computed by the equation: . In (C), the cell was exposed to a pair of square 10 ms pulses arriving with various inter-pulse intervals (Δt). The first pulse always arrived after 40 ms (and coincided in all experiments), while each secondary spike represents a specific experiment (i.e., a specific Δt). The traces were normalized so that the first spike had a peak value of −1 (corresponding to approximately −0.25 pA). https://doi.org/10.1371/journal.pcbi.1006662.g007

Computational model for medaka gonadotrophs
Under the experimental conditions, the intra-and extracellular Na + concentrations were 4 mM and 140 mM, respectively, and the temperature was 26 degrees Celsius, which gives a reversal potential (E Na = RT/F � ln([Na] ex )/ln([Na] in ) of 92 mV. The estimates of g max for activation and inactivation are indicated by the markers 'x' in Fig 8A and 8B, respectively, and markers 'o' indicate a second experiment. The dependency of g max on V hold was fitted by a Boltzmann curve: where � g max corresponds to g max estimated for the largest peak in the entire data set (i.e. at about −0.22 nA in Fig 7A and −0.18 nA in Fig 7B). The factor k determines the slope of the Boltzman curve, the exponent a corresponds to the number of activation or inactivation gates, and V � determines the voltage range where the curve rises. When a = 1 (as for inactivation), and V � listed in Table 2.

Time constants for activation and inactivation.
With three opening gates (q) and one closing gate (h), the time constants for activation and inactivation were derived by fitting the function [47] to the response curves in Fig 7A1. For low step potentials (< −40 mV), the response was too small and noisy to reveal any clear trend, and we were unable to obtain meaningful fits using the functional form of Eq 19. For this reason, only the experiments with a step potential of −40 mV and higher were used when fitting the time constants. The fitting procedure resulted in a pair of time constants (τ q and τ h ) for each step potential in the protocol, as indicated by the data points ('x' and 'o') in Fig 8C and 8D. The data points obtained by fitting Eq 21 to the traces in Fig 7A1 were sufficient to obtain a clear picture of the voltage dependence of the activation time constant (τ q ), which had a peak value at −24 mV, i.e. within the voltage-range for which there was suitable data (Fig 8C). The inactivation time constant (τ h ) was, however, monotonously decreasing over the voltage range for which there was good data. We therefore needed additional data points for the voltage dependence of τ h in the range V < −40 mV. Based on the insight from the recovery-experiments (Fig 7C), we expected inactivation to be very slow at the resting potential and below. To account for this, we introduced the three additional data points marked by '�' in Fig 8D, which assure a recovery time in the correct order of magnitude. The data points for the time constants were fitted with curves on the functional form proposed by Traub et al. [55]: Good fits to the data points were obtained with the parameter values in Table 2.

Model for high-voltage activated Ca 2+ channels
When estimating the steady-state values and time constant we followed procedures inspired from previous studies of L-type Ca 2+ channel activation, we did not use Eq 6, but used the simpler kinetics scheme I Ca = g HV A m 2 (V − E Ca ) (see e.g. [50])) assuming a constant reversal potential.
Steady-state activation. The steady-state value and time constant for m were determined from the experiments summarized in Fig 9A and 9B. To study steady-state activation, the cell was held at −60 mV for an endured period, and then stepped to different holding potentials, each for which the response current (I Ca ) was recorded (Fig 9A). Due to the small cellular size, perforated patches was used for recording the Ca 2+ currents, and the recorded currents were small and noisy. As Fig 9A shows, the I Ca responses did not follow a characteristic exponential curve towards steady state, as seen in many other experiments. Likely, this was due to I Ca comprising a complex of different HVA channels (e.g., P, Q, R, L-type) which have different activation kinetics [49][50][51][56][57][58]. In addition, some in some of the weaker responses I Ca even switched from an inward to an outward current, something that could indicate effects of ER release on the calcium reversal potential. Due to these complications, only the early part of the response was used, i.e., from stimulus onset and to the negative peak value in interval indicated by dashed vertical lines). Voltage-dependent deactivation of Ca 2+ currents (Fig 9B) was examined by measuring the tail current that followed after a 5 ms step to 10 mV when returning to voltages between −10 mV and −60 mV. The deactivation protocol was used to provide additional data points for the activation time constants (see below).
The peak current (I max ) was recorded for each holding potential in Fig 9A, and the maximum peak was observed at about 20 mV. By observations, (I Ca ) became an outward current when step potentials were increased beyond 70 mV, and based on this we assumed a reversal potential of E Ca = 70 mV. Similar to what we did for the Na + channel, the maximal conductance (g max ) for each holding potential was computed by the equation: The estimates of g max for activation are indicated by the crosses in Fig 9C. The dependency of g max on V hold was then fitted by a Boltzmann curve (Eq 20), with a = 2 activation gates, and a good fit was obtained with values for a, k and V � as in Table 3. Time constants for I Ca activation. Assuming two opening gates (m), the time constant for I Ca -activation was derived by fitting the function [47]: to the response curves in Fig 9A and 9B. The activation protocol was used to determine τ m at high step potentials (from −5 mV and upwards), where the response was not too small and noisy to reveal any clear trend (blue data points in Fig 9D). The deactivation protocol was used to determine τ m for lower step potentials (red data points in Fig 9D). Like for the Na + channel, the voltage dependence of the time constants were fitted using the functional form in Eq 22. Good fits to the data points were obtained with the parameter values in Table 3. Table 3. Parameters for Ca 2+ activation. The parameters p1-p6 are used together with Eq 22 to yield the time constants for steady state activation (in units of ms). The remaining parameters are used together with Eq 20 to obtain the steady state activation function. The curves obtained in this way describe the voltage dependence under experimental conditions, and was afterwards corrected with a liquid junction potential of −15 mV (see Fig 2B).

Implementation details
Software. The medaka gonadotroph model and the rat lactotroph model [9] were implemented using the Python interface for the NEURON simulator [59]. For the rat lactotroph model, we used a previous implementation [37]. All simulations on these two models were run using adaptive time stepping provided by the NEURON simulator. The generic murine pituitary cell model was taken from the supplementary data in the original pulblication [32], and the rat somatotroph model [27] was downloaded from https://lbm.niddk.nih.gov/sherman/ gallery/neural/somatotroph//figure4.ode. These models were coded in XPP [60], and we exported data from XPP simulations and analyzed them in Matlab to obtain the metrics plotted in Fig 5. Experimental current-clamp data (Fig 1), experimental voltage-clamp data (Figs 1, 7, 8, and 9)) and fitted ion-channel kinetics (Fig 2) were plotted using Matlab (http://se.mathworks. com/). When comparing AP shapes, simulations performed in XPP and NEURON were exported to Matlab arrays, and plotted in Matlab (Fig 6). All other plots were made in Python (http://www.python.org).
The sensitivity analysis (Fig 4) was performed by aid of the Python-based toolbox Uncertainpy [36]. The features considered (IsBursting, IsRegular, and IsNotSpiking were custom made for the analysis in the current work. Uncertainpy was run using polynomial chaos with the point collocation method (the default of Uncertainpy) and a polynomial order of five. The sensitivity analysis was based on calculating Sobol indices. Only the total-order Sobol indices were presented. A total-order Sobol index quantifies the sensitivity of a feature to a given parameter, accounting for all higher order co-interactions between the parameter and all other parameters (see [61] or the brief overview in Appendix B of [62]). For the sensitivity analysis (Fig 4), the model was run for 60,000 ms, and the first 10,000 ms were discarded to eliminate initial transients, while the remaining 50,000 ms were used in the uncertainty analysis.
The medaka gonadotroph model, and the code for generating Figs 3 and 4 are available for download (doi: 10.5281/zenodo.3359635).
Construction and tuning of the medaka gonadotroph model. The model construction and tuning of the medaka gonadotroph model were done in several steps. First, we took the membrane capacitance (C m ), the models from I leak , I K and I SK , and the reversal potentials for E Na and E K from the previous model of a rat lactotroph [9], and implemented them in the Python interface for the NEURON simulator [59]. As a starting point, we used the same parameterization as the original model [9]. However, in the original model, C m and the conductances (g X ) were given as a total capacitance and total conductances, and the implementation in NEURON required that we converted them from total currents/capacitance to currents/capacitance per membrane area. To do this, we divided them by scaling factor F scale ¼ 10:05 � 10 À 6 cm 2 ; which was selected so that the final capacitance per membrane area equal to 1μF/cm 2 . The same scaling was used when adopting I BK from [30]. Furthermore, the internal calcium handling (Eq 18) was also taken from the rat lactotroph model. However, to preserve the internal Ca 2+ dynamics under the scaling, the factors α (in Eq 18) and A (in Eq 15), which converted Ca 2+ currents into Ca 2+ concentrations were multiplied by F scale . A similar scaling was used in [37]. Second, I Na and I Ca were added to the model, with arbitrarily chosen initial conductances. Third, we tuned the model by adjusting a selection of parameters in order to obtain a firing pattern in closer resemblance with the experimental current-clamp data (Fig 1). The adjustments included: (i) τ BK was set to 5 ms, a compromise between the time constants of 3 ms and