BK channels have opposite effects on sodium versus calcium-mediated action potentials in endocrine pituitary cells

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 exclusively mediated by ICa. As a part of the previous modeling studies, a paradoxical role was identified for the big conductance K+ current (IBK), which was found to prolong the duration of ICa-mediated APs, and sometimes give rise to pseudo-plateau bursts, contrary to what one would expect from a hyperpolarizing current. 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 the, to our knowledge, first computational model of a pituitary cell that fires INa-mediated APs. Although we constrain the model to experimental data from gonadotrope cells in the teleost fish medaka (Oryzias latipes), it may likely provide insights also into other pituitary cell types that fire INa-mediated APs. In the current work, we use the model to explore how the effect of IBK depends on the AP generating mechanisms of pituitary cells. We do this by comparing simulations on the medaka gonadotrope model (two versions thereof) with simulations on a previously developed model of a rat pituitary cell. Interestingly, we find that IBK has the opposite effect on APs in the two models, i.e. it reduces the duration of already fast INa-mediated APs in the medaka model, and prolongs the duration of already slow ICa-mediated APs in the rat model. Author summary Excitable cells elicit electrical pulses called action potentials (APs), which are generated and shaped by a combination of ion channels in the cell membrane. While neurons use APs for interneuronal communication and heart cells use them to generate heart-beats, pituitary cells use APs to regulate their cytosolic Ca2+ concentration, which in turn controls their hormone secretion rate. The amount of Ca2+ that enters the pituitary cell during an AP depends strongly on how long it lasts, and it is therefore important to understand the mechanisms that control this. Depending on animal species and biological conditions, pituitary APs may be initiated either by Ca2+ channels or Na+ channels. Here, we explore the differences between the two scenarios by comparing simulations on two different computer models: (i) a previously developed model which fires Na+-based APs, adapted to data from pituitary cells in rats, and (ii) a novel model that fires Ca2+-based APs, adapted to data from pituitary cells in the fish medaka. Interestingly, we find that the role of big conductance K+ (BK) channels, which are known to affect the duration of the AP, are opposite in the two models, i.e., they act to prolong Ca2+-based APs while they act to shorten Na+-based APs.

The electrodynamics of excitable cells is generated by a combination of ion channels in 2 the plasma membrane, which are typically characterized by their voltage and/or Ca 2+ 3 dependence. While neurons primarily use action potentials as a means of interneuronal 4 communication, and cardiac cells use them to generate heartbeats, the primary role of 5 APs in endocrine pituitary cells is to regulate the cytosolic Ca 2+ concentration, which 6 in turn controls the hormone secretion rate in these cells [1]. Hormone secretion often 7 occurs as a response to hormonal stimuli from the hypothalamus, peripheral endocrine 8 glands, and other types of pituitary cells. However, many endocrine cells are also 9 spontaneously active [1][2][3][4][5][6][7][8][9][10]. The spontaneous activity is partly a means to regulate the 10 re-filling of intracellular Ca 2+ stores, but in several cells also leads to a basal release of 11 hormones. An understanding of the mechanisms regulating the electrodynamics of these 12 cells is therefore fundamental for understanding their overall functioning. 13 While neuronal APs are are predominantly mediated by TTX-sensitive Na + currents 14 (I N a ), AP generation in endocrine cells depends strongly on high-voltage-activated Ca 2+ 15 currents (I Ca ), which in addition to their role in affecting the voltage dynamics of the 16 cell, also are the main source of Ca 2+ entry through the plasma membrane [3,11,12]. In 17 some studies of endocrine cells, APs were exclusively mediated by I Ca , and the 18 spontaneous membrane excitability was insensitive or nearly so to TTX [1,2,[13][14][15][16]. In 19 other studies, APs were evoked by a combination of I Ca and I N a [4,7,17,18]. The 20 strong involvement of I Ca could explain why pituitary APs typically last longer 21 (typically some tens of milliseconds [8]) than neuronal APs (a few milliseconds), which 22 are mainly mediated by I N a . 23 All endocrine cells express I N a [8], and TTX sensitive APs can typically be triggered 24 by current injections from hyperpolarized holding potentials even in cells where they are 25 not elicited spontaneously [4,17,19,20]. The reason why the spontaneous activity is 26 TTX insensitive is likely that a major fraction of I N a is inactivated at the resting 27 membrane potential [15,16]. The reason why this is not always the case, may be that 28 I N a [3,9,[22][23][24][25][26][27]. When I N a was included in a recent modeling study, its role was mainly 40 in modulating the firing patterns but was not essential for AP firing as such [28]. The 41 main focus of these models was thus on exploring the interplay between the depolarizing 42 I Ca and various K + currents responsible for shaping the repolarization following an AP. 43 The essentials of this interplay were nicely captured in a relatively simple computational 44 model by Tabak et al. [9]. In this model, the AP upstroke was exclusively mediated by 45 I Ca , while the repolarizing phase was mediated by three K + currents (I K , I BK and 46 I SK ). The model elicited spontaneous APs with a duration of about 40 ms, which is 47 quite representative for what is seen experimentally in rat pituitary cells. In addition, 48 the model was able to explain the important and paradoxical role that the big 49 conductance Ca 2+ -activated K + current (I BK ) has in rat endocrine cells, where it was 50 found to make APs broader. The explanation to this broadening effect, which is the 51 opposite of what one would expect from a hyperpolarizing current [23], was that I BK 52 reduces the AP amplitude, and thereby prevents the (high-threshold) activation of the 53 delayed rectifying K + current (I K ), which otherwise would become a more powerful 54 hyperpolarizing current [9,23]. By inhibiting I K , I BK indirectly becomes facilitating, 55 and leads to broader APs and sometimes to voltage plateaus and so-called pseudo 56 plateau bursts [8,23], which are believed to be more efficient than regular APs in 57 evoking hormone secretion (see e.g. [29]). By manipulating the BK expression 58 experimentally and in computational models, it was further shown that the difference 59 between bursty endocrine cell types such as somatotropes and lactotropes, and regularly 60 spiking cell types such as gonadotropes and corticotrophs, could be explained by the 61 different cell types having different levels of I BK expression [9,28]. 62 The modeling studies cited above were predominantly based on data from rats, and 63 there are reasons to believe that teleost pituitary cells are different in terms of their 64 dynamical properties. Firstly, TTX-sensitive spontaneous activity has been seen in 65 goldfish resting at −60 mV [4], and TTX sensitive APs has been evoked from a holding 66 potential as high as −50 mV in pituitary cells in cod [7], suggesting that I N a may be 67 more available in resting pituitary cells in fish [4]. Secondly, data from goldfish [4] and 68 tilapia [5] indicate that the AP duration is shorter in fish (< 10 ms) compared to rats 69 (several tens of ms), and also this could indicate a stronger involvement of I N a . A third 70 difference between fish and rat pituitary cells is in the role of I BK . I BK is almost 71 absent in rat gonadotropes [23], and this was proposed as an explanation to why these 72 cells tend to be less bursty than other pituitary cell types [1,9,28]. In contrast, I BK is 73 highly expressed in medaka gonadotropes, but without making these cells bursty [12]. 74 The indication that there are differences between rat and fish pituitary cells are 75 supported by experiments presented in the current work, performed on gonadotrope 76 cells in medaka. We show that these cells elicit brief spontaneous APs that (unlike 77 spontaneous APs in rats) to a large degree are mediated by TTX sensitive Na + currents 78 (I N a ). Furthermore, we show that I BK acts to make APs narrower in medaka 79 gonadotropes, and thus have the opposite effect of what they do in rat pituitary cells. 80 Since previous computational models based on rat data seem unsuited to describe 81 the spontaneous activity of fish pituitary cells, we here present novel pituitary cell 82 models constrained to data from luteinizing hormone-producing medaka gonadotropes. 83 The main, and more general, aim of this study is to use computational modeling to 84 explore the differences between pituitary cells that fire APs exclusively mediated by I Ca

85
(like the previous rat models) and pituitary cells that fire APs that are predominantly 86 mediated by I N a (like the here developed medaka models), with a special focus on the 87 role that I BK has in the two cases. To do so, we compare simulations run on three 88 different computational models: RAT. The first model, which we refer to as RAT, is a reproduced version of the 90 model by Tabak et al. [9].  Although variations were observed, the medaka gonadotropes typically had a resting 121 potential around −50 mV, which is within the range found previously for goldfish [4]and 122 cod [10] gonadotropes. As for goldfish gonadotropes, the majority of medaka 123 gonadotropes fired spontaneous APs, and with peak voltages slightly below 0 mV. The 124 spontaneous APs were always regular spikes (i.e., not bursts) and typically had an 125 average duration between 3 and 7 ms (blue traces in Fig 1). Similar brief AP waveforms 126 have been seen in previous studies on fish [4,5,19], while the APs reported for rat 127 gonadotrope cells are typically slower, i.e. from 10-100 ms [8]. The spontaneous AP 128 activity was completely abolished by TTX application (Fig 1B), as has previously also 129 been seen for goldfish somatotropes [19]. 130 Finally, we explored how paxilline (an I BK blocker) affected the spontaneous 131 activity of gonadotrope medaka cells. In the experiment shown in Fig 1D, paxilline   132 increased the firing rate and slightly reduced the mean AP peak amplitude, but these 133 effects were not seen consistently in experiments using paxilline application. However, 134 in all experiments, paxilline application was followed by a small increase of the resting 135 membrane potential (Fig 1C), and a broadening of the AP waveform (Fig 1D2-D3). 136 Similar effects have been seen in goldfish somatotropes, where application of 137 tetraethylammonium (a general blocker of Ca 2+ gated K + currents) lead to broadening 138 of APs [19]. The effect of I BK in goldfish and medaka gonadotropes is thus to make (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 durations varied between 3 and 7 ms, with mean durations of 3.7 ms (A1), 4.9 ms (B1), 3.7 ms (C1, before TTX). In (D1), average AP durations were 4.2 ms before paxilline, and 25 ms after paxilline. Mean AP peak values were −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 gonadotrope luteinizing hormone-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.
APs narrower, which is the opposite of what was found in rat pituitary cells, where I BK 140 lead to broader APs and sometimes to burst-like activity [9,23].

141
Computational models of pituitary cells

142
In this work, we have compared simulations performed on three different pituitary cell 143 models. The first model, RAT, is a replication of the previously published model by 144 Tabak et al. [9]. This model was originally adapted to electrophysiological data from rat 145 pituitary cells. As summarized by Eq 1, RAT contains a leakage current I leak , three K + 146 currents I K , I BK and I SK , and a Ca 2+ current I r Ca . In addition, a gaussian noise 147 stimulus was added in selected simulations (see Methods).

148
C m dV dt = −(I r Ca + I K + I BK + I SK + I leak + I noise ).
The two other models are novel for this work and are the to date first pituitary cell 149 models based on teleost data, and the first to elicit APs that are predominantly 150 mediated by Na + currents. In these models, the Ca 2+ current I r Ca from RAT was 151 replaced by a pair of depolarizing currents, i.e., a novel Ca 2+ current (I m Ca ) and a Na + 152 current (I N a ), which were adapted to voltage clamp data from gonadotrope cells in 153 medaka (see Methods). We present two versions of the medaka model, both of which are 154 summarized by Eq 2: In the first version, MEDAKA 1, the three K + currents I K , I BK and I SK were identical 156 to those in RAT, so that only the depolarizing, AP generating mechanisms were 157 different. In the second teleost model, MEDAKA 2, we adjusted I K , I BK and I SK so 158 that the model had an AP shape and AP firing rate that were in better agreement with 159 the experimental data in Fig. 1. By comparing RAT to MEDAKA 1, we could then 160 explore the difference between a model with Ca 2+ APs and one with Na + /Ca 2+ APs 161 with all other mechanisms being the same. By comparing RAT to MEDAKA 2, we could 162 explore the difference two models that were more representative for experimental data 163 from rat versus medaka. The differences between RAT, MEDAKA 1 and MEDAKA 2 164 are summarized in Table 1, which lists all the parameter values that were not identical 165 in all three models, and Fig. 2 which shows the kinetics of all included ion channels.

166
For a full description of the models, we refer to the Methods section, but a brief 167 overview is given here. I N a activated in the range between −50 mV and −10 mV, with 168 half activation at −32 mV (Fig. 2A1), quite similar to what was previously found in 169 goldfish gonadotropes [4]. I N a inactivated in the range between −90 mV and −40 mV, 170 with half-inactivation at −64 mV, which was lower than in goldfish, where the 171 half-inactivation was found to be around −50 mV [4]. With the activation kinetics 172 adapted to medaka data, only 6 % of the I N a was available at the typical resting 173 potential of −50 mV. The fact that medaka still showed TTX-sensitive spontaneous 174 activity thus suggests that I N a is very highly expressed in these cells. which was based on data from rat lactotropes [33]. The high activation threshold 181 suggests that I m Ca is unsuitable for initiating spontaneous APs in medaka gonadotropes, 182 making spontaneous activity critically dependent on I N a .

183
The remaining currents (I K , I BK and I SK ) were all adopted from the RAT model 184 by Tabak et al. [9], using simplified kinetics descriptions with voltage independent time 185 constants. As explained above, these were used in their original form in RAT and 186 MEDAKA 1, but were adapted in MEDAKA 2.

187
BK currents have opposite effects on Na + -versus 188 Ca 2+ -mediated action potentials 189 It has previously been shown that I BK acts to broaden APs and promote bursting in 190 rat pituitary cells [9,23,28]. In contrast, a high I BK expression in medaka gonadotropes [12] does not make these cells bursty. On the contrary, the experiments in 192 Fig. 1D showed that medaka gonadotropes became more bursty when BK channels were 193 blocked. We hypothesized that opposite effects of BK channels in rats versus medaka 194 are not due to the BK channels being different, but rather due to (the same) BK 195 channels being activated in different ways in the two systems. 196 Here, we have used computational modeling to test this hypothesis, and have 197 compared simulations on the models RAT, MEDAKA 1 and MEDAKA 2 for different 198 levels of I BK expression (Fig. 3). Following the definition used by Tabak et al. [9], 199 plateau-proceeded broad APs of duration longer than 60 ms were defined as bursts (see 200 Methods). When BK was fully expressed in the models (i.e. had the values for g BK 201 given in Table 1), MEDAKA 1 and MEDAKA 2 fired relatively narrow APs that were 202 not classified as bursts, while RAT consistently fired bursts (Fig. 3A). A reduction of AP-events in RAT were enduring enough to be classified as bursts (Fig. 3B). When 208 g BK was set to zero, RAT became consistently non-bursty, while the two medaka 209 models became consistent bursters (Fig. 3C).

210
As summarized in Fig. 3D, the relationship between g BK and the propensity for 211 eliciting bursts in RAT was the opposite of that in the two medaka models. As RAT and 212 MEDAKA 1 differed only in terms of their AP generating mechanisms, these findings 213 suggest that BK has opposing effects on Ca 2+ generated versus Na + generated APs. The simulations in Fig. 3 had noise added to them. In order to explore why the same 216 I BK could have such different effects in RAT versus MEDAKA 1 and MEDAKA 2, we 217 ran the corresponding simulations without noise. In this way, we could ensure that all 218 aspects of the simulated traces reflected membrane mechanisms (and not random 219 fluctuations). In most of the noise-less simulations, all APs within a given spike train 220 were identical (Fig. 4). The models were quite different in terms of their AP firing 221 frequencies, but we here focus predominantly on how BK affected the shape of single 222 APs.

223
The explanation of how BK can act to broaden APs in RAT was given in previous 224 studies [9,23]. In brief, I BK acts to reduce the peak amplitude of an AP, as the 225 simulations shown here also demonstrated (compare grey curves in e.g., Fig. 4A2 and 226 E2)). This, in turn, affects I K . Since I K activates at high voltage levels, the 227 I BK -reduced AP amplitude leads to less I K activation, and via that to a net reduction 228 in the hyperpolarizing currents active during the AP downstroke. Accordingly, the 229 downstroke occurs more slowly, and the AP event gets broader. In the case with full 230 g BK expression, the interplay between the depolarizing I r Ca and hyperpolarizing I BK in 231 RAT also lead to transient oscillations during the downstroke (Fig. 4B2)).

232
To explain why I BK acted so differently in the two medaka models, we start by 233 noting that the I N a -mediated upstrokes of their APs were much faster than the 234 I r Ca -mediated AP in RAT (blue and red curves have faster upstrokes than the grey 235 curves in Fig. 4A2-E2)). I BK , having a time constant of 5 ms in all models, therefore 236 had less time to respond during the rapid AP upstrokes in MEDAKA 1 and MEDAKA 237 2, and its effect on the AP amplitude was, therefore, smaller than in RAT. Although 238 I BK did reduce the AP amplitude in all models, the amplitude reduction when going 239 from max g BK to zero g BK was by about 16 mV in RAT against only 8 mV and 7 mV 240 in MEDAKA 1 and MEDAKA 2, respectively. Another consequence of the faster AP 241 upstroke in the medaka models was that a major part of the action of I BK occurred 242 after the AP peak, i.e., during the AP downstroke. In MEDAKA 1 the effect of I BK 243 was thus shifted away from reducing the AP peak towards contributing to the AP 244 downstroke, and, in this way, I BK acted to make AP events narrower.

245
For the case with full g BK expression (Fig. 4A), MEDAKA 1 fired APs with an 246 unrealistically long duration and at a higher frequency than seen in the experimental 247 data (Fig 1A-B). Even for full BK expression, the AP peaks in MEDAKA 1 were The relationship between g BK and the propensity for eliciting bursts in the three models.
(A-D) The BK conductance is given relative to the maximum value g BK used in the respective models and was four times bigger in MEDAKA 2 than in RAT and MEDAKA 1. Simulations were otherwise the same as in Fig. 1 of the study by Tabak et al. [9], and a Gaussian noise stimulus was added.
followed by a plateau potential mediated by I m Ca (Fig 4A2). Such I Ca -mediated plateau 249 potentials have been seen in pituitary cells of Atlantic cod [7], but not in goldfish [4] or 250 tilapia [5], and was not seen during spontaneous firing in medaka ( Fig. 1A-B). In  (A2-E2) Close-ups of selected AP events. The events were aligned, i.e. t = 0 (in A2-E2) was defined independently for each trace (in A1-E1) to a point directly before the onset of the selected event. The BK conductance is given relative to the maximum value g BK used in the respective models and was four times bigger in MEDAKA 2 than in RAT and MEDAKA 1.
also responded to I BK blockage in a way that resembled that seen in experiments, i.e., 257 I BK blockage lead to a broadening of the APs (Fig 1D). The resemblance with data was 258 strongest for the simulations with partial blockage, when the AP plateaus underwent 259 oscillations that presumably reflected an interplay between I m Ca and repolarizing 260 currents activating/inactivating during the downstroke (red lines in Fig 9C2-D2)). It is 261 reasonable to assume that also in the experiments, the blockage of BK by paxilline was 262 not complete. Despite differences between MEDAKA 1 and MEDAKA 2, the effect of 263 g BK on the AP shape was similar in the two medaka models, and the opposite of that 264 in RAT.

265
In summary, I BK had an inhibitory effect on I K by reducing the AP amplitude, and 266 a collaborative effect with I K in mediating the AP downstroke. With fast AP upstrokes, 267 as in RAT, the inhibitory effect of I BK on I K was stronger than the collaborative, while 268 with slow AP upstrokes, as in MEDAKA 1 and MEDAKA 2, the collaborative effect 269 was stronger than the inhibitory. In this way, I BK acted as a mechanism that reduced 270 the duration of already brief (Na + -mediated) APs, and prolonged the duration of 271 already slow (Ca 2+ -mediated) APs.

272
Membrane mechanisms responsible for burstiness 273 In order to explore in further detail how the various membrane mechanisms affected the 274 AP firing in RAT, MEDAKA 1 and MEDAKA 2, we performed a feature-based 275 sensitivity analysis of the three models. We then assigned the maximum conductances of 276 all included currents uniform distributions within intervals ±50% of their default values 277 (Table 1), and quantified the effect that this parameter variability had on selected 278 aspects of the model output (see Methods). An exception was made for g BK , which was 279 assigned a uniform distribution between 0 and the maximum values given in Table 1), 280 i.e., from fully available to fully blocked, motivated by the fact that g BK did not have a 281 default value in RAT. We note the total order Sobol sensitivity indices considered in the 282 current analysis reflects complex interactions between several nonlinear mechanisms, 283 and that mechanistic interpretations therefore are difficult. Below, we have still 284 attempted to extract the main picture that emerged from the analysis.

285
Three features of the model responses were considered: (i) Isbursting, (ii) Isregular, 286 and (iii) Isfiring. All these features were binary, meaning e.g., that Isbursting was equal 287 to 1 in a given simulation if it contained one or more bursts, and equal to 0 if not. The 288 mean value of a feature (taken over all simulations) then represented the fraction of 289 simulations that had this feature. For example, in RAT Isbursting had a mean value of 290 0.57, Isregular had a mean value of 0.37, and Isfiring had a mean value of 0.91. This 291 means, respectively, that 57% of the parameterizations of RAT fired bursts, 37 % fired 292 regular APs, and 91% fired APs that were either bursts or regular spikes. We note that 293 the mean values for the Isbursting and isregular features in RAT and MEDAKA 2 did 294 not sum up exactly to the mean value for the AP firing feature. This was because some 295 parameterizations of these two models fired both bursts and regular APs within the 296 same simulation.

297
AP firing was seen quite robustly in RAT, which fired APs in 91 % of the sampled 298 parameter combinations (mean value of 0.91 in Fig 5A3), while APs were fired in only 299 77 % and 48% of the parameterizations of MEDAKA 1 and MEDAKA 2, respectively 300 ( Fig 5B3-C3). In MEDAKA 2 there was thus AP activity in less than half of the 301 parameterizations. This reflects that the default configuration had a resting potential 302 only slightly above the AP generation threshold, so that any parameter sample that 303 would make the cell slightly less excitable, would abolish its ability for AP generation. 304 Within the explored parameter range, RAT and MEDAKA 1 tended to be bursting, 305 meaning that bursts were seen in a larger fraction of the simulations than regular spikes 306 (mean values in Fig 5A1 and B1 are larger than in Fig 5A2 and B2). MEDAKA 2, on 307 Total-order Sobol indices   the contrary, was adapted to experimental data where AP firing under control 308 conditions ( Fig. 1A-B) was regular and with very narrow APs. MEDAKA 2 thus had a 309 high propensity for firing regular APs and elicited bursts only in (14%) of the parameter 310 combination (Fig 5C1-C2).

311
The total order Sobol indices (shown in histograms in Fig 5) quantify how much of 312 the variability (between different simulations) in the response features that were 313 explained by the variation of the different model parameters (i.e., maximum 314 conductances). For example, Isfiring in MEDAKA 1 and MEDAKA 2 were most 315 sensitive to g K and g N a (Fig 5B3 and C3)), meaning that these conductances were most 316 important for whether the model was capable of generating AP events. This result was 317 not surprising, since I K , I N a (along with the leakage current, I l ) were the only currents 318 with a nonzero activity level around rest, and thus were the ones that determined 319 whether the resting potential was above the AP firing threshold (cf. Fig 2).

320
Having other AP generation mechanisms, Isfiring in RAT was instead sensitive 321 PLOS 12/30 predominantly to g r Ca and g SK (Fig 5A3)). Due to the wide activation range of I r Ca , 322 generating the APs in RAT, this model was in principle always above the AP firing 323 threshold. The instances when AP firing was not seen in RAT then reflected lack of 324 hyperpolarizing mechanisms, such as g SK , that could repolarize the membrane potential 325 after the AP peak.

326
In general, Isbursting and Isregular were not sensitive to the same parameters as 327 Isfiring. This implies that the mechanisms determining whether the cell fired an AP or 328 not were not the same as those determining the duration of the AP once it had been 329 fired. For example, Isfiring in MEDAKA 2 was nearly insensitive to g BK , while 330 Isbursting and Isregular spiking were highly sensitive to this parameter. A little 331 simplified, we may say that g K and g N a thus determined whether MEDAKA 2 fired an 332 AP (Fig 5C3), while g BK (and to some degree g m Ca ) determined whether the AP became 333 a burst or a regular spike (Fig 5C1-C2). In MEDAKA 1, the situation was more 334 complex, and its Isbursting and Isregular were sensitive to almost all model parameters 335 (Fig 5B1-B2).

336
In RAT, Isbursting was most sensitive to g K and second most to g BK (Fig 5A1).

337
The lower sensitivity to g BK compared to g K might seem surprising, given the 338 previously proclaimed role of g BK as a burst promoter in this model [9]. However, we 339 can interpret this finding in light of the above-established understanding of burst 340 generation in RAT. We remember that I BK promoted bursting indirectly by inhibiting 341 I K , so that final determinant for whether an event became a burst was the magnitude 342 (or rather the reduced magnitude) of I K . The results from the sensitivity analysis then 343 simply indicate the direct effect on I K obtained by a variation of g K was larger than 344 the indirect effect on I K obtained by variation of g BK . 345 We note that the three different feature sensitivities in Fig 5 were not independent. 346 For example, if variations in a given parameter tended to switch the cell between 347 bursting and not firing at all, Isbursting and Isfiring would share a high sensitivity to 348 this parameter. Likewise, if variations in a given parameter instead tended to switch the 349 cell between bursting and regularly spiking, Isbursting and Isregular would share a high 350 sensitivity to this parameter. In this regard, g BK was the parameter with the cleanest 351 role as a switch between bursting and regular spiking, as in all models, the Isbursting 352 and Isregular had a high sensitivity to g BK , while Isfiring had a quite low sensitivity to 353 g BK .

354
A more general insight from the sensitivity analysis was that the propensity for AP 355 firing, regular spiking and bursting in all three models depended on a complex interplay 356 between several mechanisms. All models could be shifted from regularly spiking to 357 bursty by changes, not only in g BK (as demonstrated in Fig 3 and Fig 4), but also in 358 other conductances (such as g K or g Ca ). In all models, however, bursts were facilitated 359 by g Ca and counteracted by g K , while g BK was the more fascinating mechanism, having 360 the opposite effect on the burstiness in RAT versus the two medaka models.

374
An important goal of this work was to perform a comparison between the response 375 properties of pituitary cells that elicited I Ca -mediated APs and those that elicited 376 I N a -mediated APs. For this, we used a model based on data from rat, which elicited 377 I Ca -mediated APs, and compared it to two models (MEDAKA 1 and MEDAKA 2) 378 based on data from medaka, both of which elicited I N a . mediated APs. In this context, 379 MEDAKA 1 had the (methodological) advantage that all hyperpolarizing membrane 380 mechanisms were kept identical to those in RAT, so that RAT and MEDAKA 1 differed 381 only in terms of AP generation mechanisms. MEDAKA 2 was more strongly adapted to 382 experimental data, and had the advantage that its firing properties were more 383 representative for real medaka gonadotropes. The most interesting result that came out 384 of this comparison, was that BK currents (I BK ) had a diametrically opposite role in 385 terms of how they affect the AP shape in RAT versus the medaka models. When APs 386 were generated by I Ca (as in RAT), they were relatively slow, and I BK then acted to 387 broaden the AP events and promote bursting behavior [9,23]. On the contrary, when 388 APs predominantly were generated by I N a (as in MEDAKA 1 and MEDAKA 2), they 389 had a rapid upstroke, and I BK then acted to make the AP narrower, and prevent 390 bursting behavior (Fig 3). This suggests that I BK acts as a mechanism that 391 distinguishes between rapid and slow APs, and amplifies the difference between the two 392 by narrowing down already narrow APs while broadening already broad APs.

393
It should be noted that the effect seen in medaka gonadotropes, i.e., that I BK acted 394 to reduce AP duration, is a commonly reported role for I BK in many excitable 395 cells [35][36][37][38][39], while the burst promoting effect that I BK had in rat pituitary cells [9,23] 396 is less conventional. Furthermore, the role of I BK as a burst promoter has not been 397 found consistently in rat pituitary cells. In the study by Miranda et al. 2003, AP 398 duration in rat pituitary cells was instead found to increase when I BK was blocked with 399 paxilline [38], i.e. similar to what we found for medaka gonadotropes (Fig 1D). The 400 different effects of I BK on spike duration observed in different laboratories [23,35,38] 401 was addressed by Tabak et al. 2011 [9], who proposed possible explanations that could 402 reconcile the conflicting results. One possible explanation could be there is a variability 403 in terms of how BK channels are localized in various cells, and that BK channels that 404 are co-localized with Ca 2+ channels will respond rapidly to voltage fluctuations and 405 promote bursting, while BK channels that are not co-localized with Ca 2+ channels will 406 react more slowly to voltage fluctuations and have the opposite effect [9]. A second predominantly mediated by I Ca [9,23]. In the experiment by Miranda et al. 2003, 421 where I BK was found to shorten the AP duration (i.e., that blocking I BK lead to longer 422 APs), it was reported that this only occurred under conditions in which short APs were 423 present. It is likely that the events described in that work as short APs, were 424 I N a -mediated APs, so that the differences between the two studies by Van   Although MEDAKA 2 captured the essential firing properties of medaka 428 gonadotropes, the agreement between model and data was not perfect. For example, the 429 AP duration in MEDAKA 2 during control conditions (Fig 9C) was in the upper range 430 of that seen in the experiments (Fig 1A-B), and we were not able to obtain briefer APs 431 in the model without compromising the agreement between the experimental data and 432 other model features, such as the AP peak amplitude, afterhyperpolarization, and 433 response to I BK -blockage. The conductances selected for MEDAKA 2 ( Table 1) were 434 thus a compromise made to obtain an acceptable match to several features 435 simultaneously. The fact that we were not able to obtain a more accurate match 436 between model and data likely reflect that some of the ion channels present in the 437 model are imperfect representations of the ion channels present in the real cell. For 438 example, the simplified kinetics schemes used for I K , I BK and I SK were adopted from a 439 model of rat pituitary cells [9], and were not constrained to data from medaka 440 gonadotropes. In addition, the biological cell is likely to contain a variety of additional 441 ion channels [8] that were not included in the model. The electrophysiological experiments were conducted using the patch-clamp technique 457 on brain-pituitary slices from adult female medaka (as described in [41]). To record 458 spontaneous action potentials and Ca 2+ currents we used amphotericin B perforated 459 patch configuration, while for Na + currents we used whole cell configuration.
The passive membrane properties were the same in all models, with specific membrane 497 capacitance C m = 3.2 µF/cm 2 , and a leak conductance 0.032 mS/cm 2 .
RAT had a passive reversal potential of E leak = −50 mV. In RAT, I r Ca was quite active 499 at the resting potential, and counteracted the hyperpolarizing effect of I K which was 500 also had a non-zero activity level around the resting potential. Replacing  Below, we give a detailed description of the active ion channels. The kinetics of all 507 ion channels were summarized in Fig. 2, the model parameters that differed between the 508 models were listed in Table 1, while new parameters will be defined when introduced.

509
I N a (MEDAKA 1 and MEDAKA 2) was modeled using the standard Hodgkin and 510 Huxley-form [42]: with a reversal potential E N a = 50 mV, and gating kinetics defined by: PLOS

16/30
The steady-state activation and time constants (q ∞ , h ∞ , τ q and τ h ) were fitted to 513 voltage-clamp data from medaka gonadotropes, as described below, in the subsection 514 titled "Model for the voltage-gated Na + channels".

515
I m Ca (MEDAKA 1 and MEDAKA 2) was modelled using the Goldman-Hodgkin-Katz 516 formalism, which accounts for dynamics effect on Ca 2+ reversal potentials [43]: Here, R = 8.314J/(mol · K is the gas constant, F = 96485.3C/mol is the Faraday 519 constant, T is the temperature, which was set to 293.15 K in all simulations.
[Ca] and 520 [Ca] e were the cytosolic and extracellular Ca 2+ concentrations, respectively. The former 521 was explicitly modelled (see below), while the latter was assumed to be constant at 2 522 mM. As Eq 8 shows, we used two activation variables m. This is typical for models of 523 L-type Ca 2+ channels (see e.g. [44][45][46][47]), which are the most abundantly expressed HVA 524 channels in the cells studied here [11]. The steady-state activation and time constant 525 (m ∞ and τ h ) were fitted to voltage-clamp data from medaka gonadotropes, as described 526 below, in the subsection titled "Model for high-voltage activated Ca 2+ channels". We with instantaneous steady-state activation: where s m = is a slope parameter and v m the potential at half-activation.

534
The directly rectifying K + channel (all models) was modelled as with reversal potential E K = −75 mV, and a time dependent activation variable 536 described by The constant (voltage-independent) time constant τ K was 30 ms in RAT and MEDAKA 538 1, and 5 ms in MEDAKA 2, and the steady-state activation was described by: with a slope parameter s n = 10 mV, and half-activation v n = −5 mV.

540
The BK-channel (all models) was modelled as The activation kinetics was: with a constant (voltage-independent) activation time constant τ BK = 5 ms. The 543 steady-state activation was given by: with a slope parameter s m = 12 mV, and half-activation v m with values −20 mV in 545 RAT and MEDAKA 1 and −15 mV in MEDAKA 2. We note that BK was modeled as 546 a voltage (and not Ca 2+ ) dependent current. The rationale behind this simplification 547 was that BK activation depends on the Ca 2+ concentration in highly localized Ca 2+ 548 nanodomains, where it is co-localized with, and depends on influx through, high-voltage 549 activated Ca 2+ channels. This influx is in turn is largely determined by the membrane 550 potential and reaches an equilibrium within microseconds [9], making the BK channel 551 effectively voltage dependent.

552
Finally, the SK channel (all models) was modeled as: with an instantaneous, Ca 2+ dependent, steady-state activation: where [Ca] denotes the cytosolic Ca 2+ concentration, and k s was a half-activation 555 concentration of 0.4 µM.

556
I m Ca and I SK were dependent on the global cytosolic Ca 2+ concentration. This was 557 modelled as a simple extrusion mechanism, receiving a source through I m Ca (for 558 MEDAKA 1 and MEDAKA 2) and I r Ca for RAT, and with a concentration dependent 559 decay term assumed to capture the effects of various extrusion and buffering 560 mechanisms: Here, f c = 0.01 is the assumed fraction of free Ca 2+ in the cytoplasm, 562 α = 0.004725 mM · cm 2 /µC converts an incoming current to a molar concentration, and 563 k c = 0.12 ms −1 is the extrusion rate [9]. Due to the requirements from the NEURON 564 simulator, the parameters from the original model by Tabak et al. [9] were converted 565 from total currents/capacitance to currents/capacitance per membrane area, using a cell 566 body with a membrane area of π · 10 −6 cm 2 (see 567 https://github.com/ReScience/ReScience-submission/pull/53).

568
In the simulations in Fig 3, a noise input was added, described by where the noise amplitude (A noise ) was 4 pA, and η was a random process drawn from 570 a normal distribution [9].

571
Model for the voltage-gated Na + channels 572 The steady-state values and time courses of the gating kinetics were determined using 573 standard procedures (see e.g. [32,42,48,49]), and was based on the experiments 574 summarized in Fig 6. To determine activation, the cell was held at −60 mV for an 575 endured period, and then stepped to different holding potentials between −80 to 100 576 mV with 5 mV increments (Fig. 6A2), each for which the response current (I N a ) was 577 recorded (Fig. 6A1). The inactivation properties of Na + were investigated using 578 stepwise pre-pulses (for 500 ms) between −90 and 55 mV with 5 mV increments before 579 recording the current at −10 mV (Fig. 6B2). The resulting Na + current then depended 580 on the original holding potential (Fig. 6B1). Finally, the recovery time for the Na + 581 current was explored by exposing the cell to a pair of square pulses (stepping from a 582 holding potential of −60 mV to −10 mV for 10 ms) separated by a time interval ∆t 583 (Fig. 6C2). The smallest ∆t was 10 ms, and after this, ∆t was increased with 100 ms in 584 each trial. The cell responded to both pulses by eliciting Na + -current spikes (Fig. 6C2). 585 When ∆t was small, the peak voltage of the secondary spike was significantly reduced 586 compared to the first spike, and a full recovery required a ∆t in the order of 1/2-1 s. . 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).

588
To determine steady-state activation, the peak current (I max ) was determined for each 589 holding potential in Fig. 6A, and the maximum peak was observed at about −10 mV. 590 For inactivation, the peak current (I max ) was recorded for each holding potential in Fig. 591 6B. In both cases, the maximal conductance (g max ) for each holding potential was 592 computed by the equation: Under the experimental conditions, the intra-and extracellular Na + concentrations dependency of g max on V hold was fitted by a Boltzmann curve: whereḡ max corresponds to g max estimated for the largest peak in the entire data set (i.e. 600 at about −0.22 nA in Fig. 6A and −0.18 nA in Fig. 6B). The factor k determines the 601 slope of the Boltzman curve, the exponent a corresponds to the number of activation or 602 inactivation gates, and V * determines the voltage range where the curve rises. When 603 a = 1 (as for inactivation), V * equals V 1/2 , i.e. the voltage where f bz has reached its half 604 maximum value. With a higher number of gates, V * = V 1/2 + k · ln(2 1/a − 1). Eq 23 605 gave a good fit for the steady state activation m 3 inf (Fig 7A) and the steady state 606 inactivation h inf (Fig 7B) with the parameter values for a, k and V * listed in Table 2.

607
Time constants for activation and inactivation 608 With three opening gates (q) and one closing gate (h), the time constants for activation 609 and inactivation were derived by fitting the function [42] 610 to the response curves in Fig. 6A1. For low step potentials (< −40 mV), the response 611 was too small and noisy to reveal any clear trend, and we were unable to obtain which there was suitable data (Fig. 7C). The inactivation time constant (τ h ) was, 620 however, monotonously decreasing over the voltage range for which there was good data. 621 We therefore needed additional data points for the voltage dependence of τ h in the 622 range V < −40 mV. Based on the insight from the recovery-experiments (Fig. 6C), we 623 expected inactivation to be very slow at the resting potential and below. To account for 624 this, we introduced the three additional data points marked by '*' in Fig. 7D, which 625 assure a recovery time in the correct order of magnitude.

626
The data points for the time constants were fitted with curves on the functional 627 form proposed by Traub et al. [50]: Good fits to the data points were obtained with the parameter values in Table 2 Table 2. Parameters for Na + activation. The parameters p1-p6 are used together with Eq. 25 to yield the time constants for steady state activation and inactivation (in units of ms). The remaining parameters are used together with Eq. 23 to obtain the steady-state activation and inactivation functions. The curves obtained in this way describe the voltage dependence under experimental conditions, and was afterwards corrected by subtracting the liquid junction potential of −9 mV (see Fig 2A).
Model for high-voltage activated Ca 2+ channels 630 When estimating the steady-state values and time constant we followed procedures 631 inspired from previous studies of L-type Ca 2+ channel activation, we did not use Eq 8, 632 but used the simpler kinetics scheme I Ca = g HV A m 2 (V − E Ca ) (see e.g. [45])) assuming 633 a constant reversal potential.

634
Steady-state activation 635 The steady-state value and time constant for m were determined from the experiments 636 summarized in Fig. 8A-B. To study steady-state activation, the cell was held at −60 637 mV for an endured period, and then stepped to different holding potentials, each for 638 which the response current (I Ca ) was recorded (Fig. 8A). Due to the small cellular size, 639 perforated patches was used for recording the Ca 2+ currents, and the recorded currents 640 were small and noisy. As Fig. 8A shows, the I Ca responses did not follow a 641 characteristic exponential curve towards steady state, as seen in many other 642 experiments. Likely, this was due to I Ca comprising a complex of different HVA 643 channels (e.g., P, Q, R, L-type) which have different activation kinetics [44][45][46][51][52][53]. 644 In addition, some in some of the weaker responses I Ca even switched from an inward to 645 an outward current, something that could indicate effects of ER release on the calcium 646 reversal potential. Due to these complications, only the early part of the response was 647 used, i.e., from stimulus onset and to the negative peak value in interval indicated by 648 dashed vertical lines). Voltage-dependent deactivation of Ca 2+ currents (Fig. 8B) was 649 examined by measuring the tail current that followed after a 5 ms step to 10 mV when 650 returning to voltages between −10 mV and −60 mV. The deactivation protocol was 651 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. 8A, and the 653 maximum peak was observed at about 20 mV. By observations, (I Ca ) became an 654 Table 3. Parameters for Ca 2+ activation. The parameters p1-p6 are used together with Eq. 25 to yield the time constants for steady state activation (in units of ms). The remaining parameters are used together with Eq. 23 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).
Adjustments of I K , I BK and I SK made in MEDAKA 2

673
In MEDAKA 2, I K , I BK and I SK were adjusted in order to obtain a model whose  Firstly, I BK activation was shifted by 5 mV relative to RAT and MEDAKA 1. This 678 caused I BK activation to occur closer to the peak, it reduced the effect that blockage of 679 Under the control conditions, the AP firing rate was 0.7 Hz (A). Partial BK-blockage by paxilline made the cell bursty (B), and the burst shape depended on how large fraction of BK channels that were blocked (gray, yellow, pink, and green curves in C). Regular APs peaked at −4.4 mV (blue line in C) and had a duration (taken at half max amplitude between −50 mV and AP peak) of 6.4 ms. These values were within the range of peak and duration values observed in the data. As in the TTX-data (Fig 1C), blockage of g N a abolished AP generation completely (red line in A).
I BK had on the AP amplitude, and was necessary in order to obtain oscillations during 680 the plateaus following APs (green, yellow and pink curves in Fig 9C). 681 Secondly, the kinetics of I K in RAT was based on data from rat lactotropes, where 682 I K activation had a fast (3.7 ms) and a slow (30 ms) component [54]. The time constant 683 τ K in RAT was 30 ms, and thus reflected the slow component. For I K to have an effect 684 on the repolarization of faster I N a -mediated APs, τ K was reduced to 5 ms in MEDAKA 685 2, a value closer to the fast component seen in the data [54]. 686 Thirdly, with the maximum value of g BK in RAT or MEDAKA 1, APs were 687 exceeded by a plateau. I Ca -mediated plateau potentials have not been observed in 688 goldfish [4] or tilapia [5], and was not seen during spontaenous firing in medaka ( Fig.   689 1A-B). To remove the plateaus, g K and g BK were increased by factors 1.4 and 4, 690 respectively, compared to RAT and MEDAKA 1, which gave MEDAKA 2 narrow 691 regular APs under control conditions (blue curve in Fig 9C). 692 Fourthly, the faster repolarization obtained by increasing g BK and g K increased the 693 firing rate of MEDAKA 2 to unrealistically high values. Similar effects of I BK on 694 increasing firing rates have been seen in other systems [39]. A lower and more realistic 695 firing rate was obtained by increasing g SK by a factor 3.

696
With the adjustments described above, MEDAKA 2 responded to partial I BK 697 blockage in a way that resembled the experiments with paxilline application (compare 698 Fig 9C and Fig 1D).

707
All simulations (used in Figures 3, 4, 5 and 9) were run for 60,000 ms (although 708 briefer intervals were shown in the figures). The first 10,000 ms were discarded to 709 eliminate initial transients, while the remaining 50,000 ms were used in the uncertainty 710 analysis (Fig 5) and to calculate the burstiness factor (Fig 3). Simulations with noise 711 (Fig 3) were run using a fixed time step dt = 0.25 ms, while simulations without noise 712 were run using adaptive time stepping provided by the NEURON simulator.

713
The sensitivity analysis in Fig 5 was performed by aid of the Python-based toolbox 714 Uncertainpy [34]. The features considered (Bursting, Regular spiking, and AP firing as 715 defined below) were custom made for the analysis in the current work. Uncertainpy was 716 run using polynomial chaos with the point collocation method (the default of 717 Uncertainpy) and a polynomial order of five. The sensitivity analysis was based on 718 calculating Sobol indices. Only the total order Sobol indices were presented in this work. 719 A total order Sobol index quantifies the sensitivity of a feature to a given parameter, 720 accounting for all higher order co-interactions between the parameter and all other 721 parameters (see [56] or the brief overview in Appendix B of [57]).

722
The models (RAT, MEDAKA 1 and MEDAKA 2), and the code for generating Below, the various metrics used throughout this article are defined.

726
• AP width: In control conditions, AP width was defined as the time between the 727 upstroke and downstroke crossings of the voltage midways between −50 mV and 728 the peak potential.

729
• Event duration: A metric proposed by [9], tailored to capture duration of 730 longer events such as bursts. The voltage trace was first normalized so that the 731 minimum voltage was set to 0 and the maximum voltage to 1. The start of an 732 event was defined when the voltage crossed an onset threshold 0.55, and the end 733 of the event was defined when the voltage then crossed a termination threshold 734 0.45. The duration of the event was the time from onset to termination.

735
• Burst: An AP event with a duration of 60 ms or more was classified as a burst. 736 • Burstiness Factor: The fraction of AP events within a given simulation that 737 were bursts. This metric was most relevant in simulations with noise input (Fig 3). 738 In simulations without noise added, all APs within a given simulation tended to 739 be close to identical, so that the burstiness factor was either 0 or 1 (although 740 there were exceptions).

741
• Isbursting : In the sensitivity analysis (Fig 5 ), the feature Isbursting was a 742 binary variable that was 1 for simulations that contained bursts, and 0 for 743 simulations that did not.

744
• Isregular : In the sensitivity analysis, the feature Isregular was a binary variable 745 that was 1 for simulations that contained regular APs, and 0 for simulations that 746 did not.

25/30
• Isfiring : In the sensitivity analysis, the feature Isfiring was a binary variable 748 that was 1 for simulations that contained any APs (regular or bursts) and 0 for 749 simulations that did not.