Advertisement

  • Loading metrics

Open Access

Peer-reviewed

Research Article

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

  • Geir Halnes ,

    Roles Conceptualization, Formal analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing – original draft, Writing – review & editing

    geir.halnes@nmbu.no

    Affiliations Faculty for Science and Technology, Norwegian University of Life Sciences, Ås, Norway, Centre for Integrative Neuroplasticity, University of Oslo, Oslo, Norway

    ⨯
  • Simen Tennøe,

    Roles Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – review & editing

    Affiliations Centre for Integrative Neuroplasticity, University of Oslo, Oslo, Norway, Department of Informatics, University of Oslo, Oslo, Norway

    ⨯
  • Trude M. Haug,

    Roles Conceptualization, Methodology, Supervision, Validation, Writing – review & editing

    Affiliation Institute of Oral Biology, University of Oslo, Oslo, Norway

    ⨯
  • Gaute T. Einevoll,

    Roles Conceptualization, Funding acquisition, Project administration, Supervision, Writing – review & editing

    Affiliations Faculty for Science and Technology, Norwegian University of Life Sciences, Ås, Norway, Department of Physics, University of Oslo, Oslo, Norway

    ⨯
  • Finn-Arne Weltzien,

    Roles Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing

    Affiliation Department of Basic Sciences and Aquatic Medicine, Norwegian University of Life Sciences, Campus Adamstuen, Oslo, Norway

    ⨯
  • Kjetil Hodne

    Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Validation, Writing – original draft, Writing – review & editing

    Affiliation Department of Basic Sciences and Aquatic Medicine, Norwegian University of Life Sciences, Campus Adamstuen, Oslo, Norway

    ⨯

Abstract

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.

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. Since one type of ion channels is permeable to Ca2+ ions, there is typically an influx of Ca2+ during an AP. Pituitary cells therefore use AP firing to regulate their cytosolic Ca2+ concentration, which in turn controls their hormone secretion rate. The amount of Ca2+ that enters during an AP depends strongly on how long it lasts, and it is therefore important to understand the mechanisms that control this. Pituitary APs are generally mediated by a combination of Ca2+ channels and Na+ channels, and the relative contributions of from the two vary between cell types, animal species and biological conditions. Previous computer models have predominantly been adapted to data from pituitary cells which tend to fire Ca2+-based APs. Here we develop a new model, adapted to data from pituitary cells in the fish medaka, which APs that are predominantly Na+-based, and compare its dynamical properties to the previous models that fire Ca2+-based APs.

Citation: Halnes G, Tennøe S, Haug TM, Einevoll GT, Weltzien F-A, Hodne K (2019) A computational model for gonadotropin releasing cells in the teleost fish medaka. PLoS Comput Biol 15(8): e1006662. https://doi.org/10.1371/journal.pcbi.1006662

Editor: Duncan J. MacGregor, University of Edinburgh, UNITED KINGDOM

Received: November 21, 2018; Accepted: August 1, 2019; Published: August 22, 2019

Copyright: © 2019 Halnes et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the manuscript.

Funding: This work was funded by the Research Council of Norway (https://www.forskningsradet.no) via (1) the BIOTEK2021 Digital Life project`DigiBrain’, grant no 248828 (received by GTE and F-AW), and (2) the Aquaculture program, grant no 244461 (received by KH). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

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 Ca2+ dependent gating properties. The primary role of APs in endocrine pituitary cells is to regulate the cytosolic Ca2+ 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–10]. The spontaneous activity is partly a means to regulate the re-filling of intracellular Ca2+ 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 (INa), AP generation in endocrine cells depends strongly on high-voltage-activated Ca2+ currents (ICa), which in addition to their role in affecting the voltage dynamics of the cell, also are the main source of Ca2+ entry through the plasma membrane [3, 11, 12]. In some studies of endocrine cells, APs were exclusively mediated by ICa, and the spontaneous membrane excitability was insensitive or nearly so to TTX [1, 2, 13–16]. In other studies, APs were evoked by a combination of ICa and INa [4, 7, 17–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 ICa and INa components. The strong involvement of ICa 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 INa.

All endocrine cells express INa [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 INa 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, INa tends to be inactivated and the spontaneous activity TTX insensitive (see reviews in [8, 22]). Models based on rat data have therefore typically excluded INa [3, 9, 23–30]. As TTX-sensitive spontaneous APs have been seen in mice corticotrophs [19], INa 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 INa in these models was mainly in modulating the firing patterns, and it was not essential for AP firing as such [31, 32]. Furthermore, ICa and INa 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 INa 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 INa. A third difference between fish and rat pituitary cells may be in the role of the big conductance K+ current (IBK), which has been shown to have a paradoxical role in some rat pituitary cells, i.e., it can 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 [9, 25]. IBK 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, IBK 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 (INa). Furthermore, we show that IBK 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 ICa and INa 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 ICa and INa, 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 IBK blocker paxilline.

By comparing the medaka gonadotroph model, which predominantly fires INa-mediated APs, with three previous models of murine pituitary cells, which predominantly [32] or exclusively [9, 27] fire ICa-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 IBK 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 IBK in shaping pituitary APs [8], and suggest that the effect of IBK on APs is mainly determined by the timing of IBK-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 INa-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–35].

Results

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 Ca2+ 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.

thumbnail
Download:
Fig 1. Experimental voltage recordings.

(A1-B1) Spontaneous AP firing in two selected cells. (A2-A3) Close-ups of selected APs from the cell in A1. (B2-B3) Close-ups of selected APs from the cell in A2. (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.

https://doi.org/10.1371/journal.pcbi.1006662.g001

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 IBK 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 IBK 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 IBK 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: (1) The three K+ currents, IK, IBK and ISK, 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. IK denotes the delayed rectifier K+ channel [9], ISK denotes the small-conductance K+ channel, activated by the intracellular Ca2+ concentration [9], and IBK denotes the big-conductance K+ channel. The latter was assumed to be both voltage and Ca2+-dependent. As it is often co-localized high-voltage-activated Ca2+ channels, it was assumed to sense a domain-Ca2+ concentration proportional to ICa [32].

The depolarizing membrane currents consisted of a high-voltage activated Ca2+ current (ICa) and a Na+ current (INa), both of which were novel for this model, and adapted to new voltage-clamp data from gonadotroph cells in medaka (see Methods). INa 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]. INa inactivated in the range between −90 mV and −40 mV, with half-inactivation at −64 mV, which was lower than 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 INa was available at the typical resting potential of −50 mV. The fact that medaka still showed TTX-sensitive spontaneous activity thus suggests that INa is highly expressed in these cells. In comparison, in the generic murine pituitary model [32], INa activation required depolarization to voltages far above the resting potential (red curves, Fig 2A1)), meaning that this model could not elicit spontaneous INa-based APs.

thumbnail
Download:
Fig 2. Ion channel kinetics.

(A1) INa in models of a medaka gonadotroph and a generic murine pituitary cell (G). INa had three activation gates (q3) and one inactivation gate (h). (B1) ICa activation in models of medaka gonadotrophs, generic murine pituitary cells (G), rat lactotrophs (RL), and rat somatotrophs (RS). Two activation gates were used in medaka (m2), and one in the other models. (C1) IK had one activation gate (n). (D1) IBK, had one activation gate (f), depending on voltage and domain [Ca2+], the latter assumed to be proportional to ICa. Results shown for ICa = 0, 10, 30 and 100 pA. (E1) ISK was Ca2+ activated with one activation gate s. (A2-B2) Voltage-dependent activation time constants were computed for INa (A2) and ICa (black curves), while all murine model used fixed time constants (red curves). Voltage-independent activation-time constants were used for IK (C2), IBK (D2) and ISK (E2). (A-E).

https://doi.org/10.1371/journal.pcbi.1006662.g002

Both INa and ICa had fast activation in medaka gonadotrophs, INa being slightly faster with a time constant of about 0.5–0.8 ms in the critical voltage range (Fig 2A2), whereas ICa had a time constant >1 ms in the critical voltage range (Fig 2B2). ICa 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 ICa in the rat models (colored curves). The high activation in medaka gonadotrophs threshold suggests that ICa is unsuitable for initiating spontaneous APs, making spontaneous activity critically dependent on INa, unlike in all rat models, where ICa 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 gNa was set to zero, mimicking the effect of TTX application in the experiments.

thumbnail
Download:
Fig 3. Effects of IBK on AP shape.

The spontaneous activity of the medaka gonadotroph model for different levels of BK-expression. (A1) Spontaneous activity under control conditions, where gBK was maximally expressed (blue curve, gBK = 0.31mS/cm2). AP firing was completely abolished by setting gNa = 0, mimicking the effect of TTX (red curve). (B1-F1) Simulations with gBK reduced to fractions (indicated above panels) of the maximum value. Reductions in gBK consistently lead to broader APs. (A2-F2) Close-ups of the first APs in seen in A1-F1.

https://doi.org/10.1371/journal.pcbi.1006662.g003

A series of previous models have shown that IBK may act to broaden APs and promote bursting in rat pituitary cells [9, 25, 27, 32]. In contrast, a high IBK 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 (gBK) was reduced relative to its value during control (Fig 3B–3F). Generally, a reduction in gBK 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 gBK 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 ICa and repolarizing currents activating/inactivating after the initial AP peak. When gBK 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 gBK, 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 therefore are difficult. Below, we have still attempted to extract the main picture that emerged from the analysis.

thumbnail
Download:
Fig 4. Feature-based sensitivity analysis.

Sensitivity to variations in the maximum ion-channel conductances. The analysis summarizes a large number of simulations where the maximum conductances of all ion channels were varied within intervals ±50% of their original values. An exception was made for gBK, which was varied between 0 and 0.31 mS/cm2. Features were binary (0 or 1), and (A) IsBursting = 1 for simulations that elicited one or more bursts, (B) IsRegular = 1 for simulations that elicited one or more regular APs, and (C) IsNotSpiking = 1 for simulations that did not elicit any AP events. (A-C) Histograms depict the total-order Sobol sensitivity indices which quantify how much of the variability in a response features that is explained by the variation of a given model parameter, including all its co-variances with other parameters. The analysis was performed by aid of the recently developed toolbox Uncertainpy [36] (see Methods for details).

https://doi.org/10.1371/journal.pcbi.1006662.g004

thumbnail
Download:
Table 1. Conductances in the default parameterization of the medaka gonadotroph model.

gBK was varied between simulations, and had values between 0 and the (maximum) value given the table. * gCa had the units of a permeability.

https://doi.org/10.1371/journal.pcbi.1006662.t001

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 IsNotSpiking 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 gNa was small for IsBursting, but quite large for IsRegular and IsNotspiking, it then means that gNa 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 gBK, while IsBursting and IsRegular had a high sensitivity to gBK. A little simplified, we can thus say that gNa determined whether the model fired an AP (Fig 4C), while gBK 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 gK, which indicates that gK played multiple roles for the firing properties of the model. Firstly, IK had a nonzero activity level around rest (cf. Fig 2C1), and was important (along with INa and the leakage current, Il) for determining whether the resting potential was above the AP threshold, hence the high sensitivity of IsNotSpiking to gK. Having a broad activation range, IK was also important for repolarizing the membrane potential after the AP peak, and thus for the duration of the AP. Therefore, also IsBursting had a high sensitivity to gK.

The mechanisms behind burst generation are reflected in the sensitivity of IsBursting to gBK, gK and gCa. Here, ICa is responsible for mediating the plateaus that prolong APs into possible bursts, while gBK and gK may prevent bursts by facilitating a faster down-stroke. The interaction between gBK and gK in mediating the AP downstroke is complex, as we comment on further in the next subsection. The sensitivity to the last K+-channel, gSK, was very low in all the features considered here. gSK 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, IBK 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 gBK 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 INa, while in the murine models, APs were predominantly mediated by ICa. Despite several differences, all models contained IBK and IK, which were the most important ion channels for mediating the AP downstroke.

thumbnail
Download:
Fig 5. Effects of IBK on AP shape in different pituitary cell models.

(A) Four versions of the medaka gonadotroph model were studied, differing in terms of the BK activation time constant (τBK). The default parameterization had τBK = 3 ms (orange curve). 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 gBK in all models. (A2-D2) An increase in gBK 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 gBK 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 gBK relative to a reference value gBKref, which was taken to be the default BK-conductance in the respective models.

https://doi.org/10.1371/journal.pcbi.1006662.g005

In all models, an increase in gBK 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 IBK-activation time constant (τBK). In the default configuration, τBK was set to 3 ms (orange curve in Fig 5A1). With a slower τBK, IBK 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, IBK activation largely occurred during the AP upstroke, and IBK then had a larger effect on the AP-peak value (blue curve in Fig 5A1). In general, IBK 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 IBK activation and the AP peak.

The effect on gBK 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 gBK affected the AHP depth (Fig 5A3–5D3). The AHP was predominantly due to the joint effect of the two hyperpolarizing currents, gBK and gK. It may therefore seem counterintuitive that, for low gBK, an increase in gBK actually decreased the AHP (less hyperpolarization means higher AHP voltages). The explanation lies in the simultaneous effect that gBK had on reducing the AP-peak voltage (Fig 5A1–5D1). Lower AP-peak values generally meant less IK activation [9, 25], as this current activates at high voltage levels. Hence, gBK 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 gK activation. When gBK was small, the AHP was predominantly mediated by gK, and the indirect effect dominated, so that an increase in gBK reduced the AHP up to a certain point, where the direct effect to took over, so that a further increase in gBK enhanced the AHP.

The dual role of gBK is also reflected in the effect that gBK had on the AP width (Fig 5A2–5D2). An increase in gBK could either lead to briefer APs, through the direct hyperpolarizing effect of IBK, or broader APs, through the indirect effect of gBK reducing IK activation. This dual role of gBK is seen most clearly in the rat lactotroph model (Fig 5B2). For low values of gBK, the direct effect dominated, and the AP width decayed monotonically with increasing gBK up to a certain threshold value, where a further increase in gBK gave a sharp transition to very broad APs (pseudo-plateau bursts). The paradoxical role of IBK 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 gBK on the AP shape was seen in the generic pituitary model, although the effect of gBK on narrowing APs for low gBK was there very small (Fig 5C2). In the rat somatotroph model, the indirect effect dominated for all values of gBK, and the AP width increased monotonically with increasing gBK (Fig 5D2). Oppositely, in the default parameterization of the medaka gonadotroph model, the direct effect dominated for all values of gBK, and the AP width decreased monotonically with increasing gBK (orange curve, Fig 5A2). Only by decreasing τBK to unrealistically low values, gBK 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 gBK values, neither the AP width or AHP increased monotonically with gBK 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, IBK facilitated bursting in all the considered models based on murine data, but not in the default parameterization of the medaka gonadotroph model. As the IBK 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 IBK in the medaka gonadotroph model versus the murine pituitary cell models was due to differences in AP shape, rather than IBK 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 INa 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 ICa and INa (purple curve in Fig 6). Hence, the addition of INa 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 gBK.

thumbnail
Download:
Fig 6. Action potential upstrokes in different pituitary cell models.

The AP upstroke was steepest in the medaka gonadotroph model (where it was mediated by INa), and slower in the murine pituitary models where they were predominantly mediated by ICa. The models considered were (blue) the default parameterization of the medaka gonadotroph model, (red) the rat lactotroph model, a generic murine pituitary cell model with (purple) and without (orange) sodium conductances, and (green) the rat somatotroph model. In all models, the BK conductance was set to zero.

https://doi.org/10.1371/journal.pcbi.1006662.g006

In the remaining models, APs were mediated solely by ICa, and had slower upstroke. Hence, in the murine models, IBK had more time to activate during the AP upstrokes, which explains why IBK indirectly could promote bursting in these models, by reducing AP-peak value and thereby IK 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 IBK 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 gBK, while it fired regular AP for intermediate gBK value. Also a previous model of a rat corticotroph showed such bursting behaviour in two disjoint regions in gBK-space (see Fig 3 in [31]).

In summary, IBK had an inhibitory effect on IK by reducing the AP amplitude, and a collaborative effect with IK in mediating the AP downstroke. With slow AP upstrokes, as in the murine pituitary cell models, the inhibitory effect of IBK on IK could dominate, and IBK 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 IBK always dominated, and IBK consistently facilitated narrower APs. In a broader scope, this suggests that IBK 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 (INa) 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 INa is of lesser importance, and where AP generation is predominantly mediated by high-voltage activated Ca2+ currents [3, 9, 23, 25–28]. To our knowledge, we have in the current work presented the first models that describe pituitary cells under conditions where AP generation is INa-mediated. The model was adapted to experimental data from LH-producing gonadotrophs in medaka, whose spontaneous activity is highly INa-dependent. Voltage-clamp data was used to develop models for the activation kinetics for INa and ICa currents, and the firing properties of the model were further adapted to current-clamp 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 INa-mediated APs, and three models of murine pituitary cells which fired APs that were exclusively mediated by ICa [9, 27], or by a combination of ICa and INa [32]. The most interesting result that came out of this comparison was that IBK 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 IBK had sufficient time to activate during the AP upstroke, and thus required either a very fast IBK 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 IBK 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 IBK in many excitable cells [40–44], while the AP-broadening and burst promoting effect that IBK is less conventional.

The role of IBK as a burst promoter has not been found consistently in rat lactotrophes. In the study by Miranda et al. 2003, AP width in rat GH3, a widely used model for pituitary lactotrophs, was instead found to increase when IBK was blocked with paxilline [43], i.e., similar to what we found for medaka gonadotrophs (Fig 1D). The different effects of IBK 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 Ca2+ channels will respond rapidly to voltage fluctuations and promote bursting, while BK channels that are not co-localized with Ca2+ 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 IBK 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 IBK in lactotrophs. Firstly, the fact that IBK has affects the AP shape differently in different cells does not by necessity reflect differences in IBK kinetics or localization. Simulations on the rat lactotroph model showed that the same IBK 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 gBK to very low values, or by increasing it to very high values. This dual effect of IBK 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 gBK blockage and for large gBK expression, while they fired regular APs for intermediate gBK expression. Hence, in general, gBK blockage could affect the AP width in either way, depending on the initial level of gBK expression, and conflicting conclusions regarding the role of IBK could reflect that variations in gBK expression under control conditions. Secondly, the way on which IBK will affect the APs in a given cell can not be predicted from gBK 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 IBK, while APs with a slower upstroke were prone to be prolonged by IBK, as also suggested in the dynamic clamp experiments by Tabak et al. [9]. Putatively, INa mediated APs will generally have a steeper upstroke than ICa mediated APs, as was the case in the models studied here. If this is the case, the role of IBK in a given cell may thus largely be determined by which membrane mechanisms that mediate its AP upstroke, and especially the degree to which INa is involved, which is highly resting-potential dependent, and likely to depend strongly on experimental conditions. Differences in INa involvement could in principle explain the conflicting experiments on rat lactotrophs [25, 43]. In the experiment by Van Goor et al. 2001, where IBK was found to broaden APs, APs were predominantly mediated by ICa [9, 25]. In the experiment by Miranda et al. 2003, where IBK was found to narrow APs (i.e., blocking IBK 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 INa-mediated, so that the differences between the studies by Van Goor et al. 2011 and 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 IBK-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 IK, IBK and ISK 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 INa-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, INa-based APs are elicited by several pituitary endocrine cell types and in several animal species, depending on biological conditions [4, 7, 17, 18, 33–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].

Methods

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 Ca2+ 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 CaCl2, 1.2 mM MgCl2 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 Ca2+ 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 Ca2+ currents, we substituted NaCl with 120 mM choline-Cl and added 20 mM Ca2+, 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 MΩ for perforated patch recordings and 6-7 MΩ 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 C6H13NO4S (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 Ca2+ 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 MΩ 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 INa (Fig 2A). A liquid junction potential of about −15 mV was calculated for the experimental (voltage-clamp) conditions used for recording Ca2+ currents, and was corrected for in the kinetics model for (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: (2) The membrane capacitance was set to the standard value Cm = 1μF/cm2, and the leak conductance was described by (3) with a reversal potential Eleak = −45 mV. Due to a nonzero-activation of IK 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. INa was modeled using the standard Hodgkin and Huxley-form [47]: (4) with a reversal potential ENa = 50 mV, and gating kinetics defined by: (5) The steady-state activation and time constants (q∞, h∞, τ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”.

ICa was modelled using the Goldman-Hodgkin-Katz formalism, which accounts for dynamics effect on Ca2+ reversal potentials [48]: (6) with (7) 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 Ca2+ 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 Ca2+ channels (see e.g. [49–52]), which are the most abundantly expressed HVA channels in the cells studied here [11]. The steady-state activation and time constant (m∞ and τh) were fitted to voltage-clamp data from medaka gonadotrophs, as described below, in the subsection titled “Model for high-voltage activated Ca2+ channels”. We note that gCa 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 (8) with reversal potential EK = −75 mV, and a time dependent activation gate described by (9) The steady-state activation was described by: (10) with a slope parameter sn = 10 mV, and half-activation vn = −5 mV. The model for IK 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]: (11) The activation kinetics was given by: (12) The constant (voltage-independent) activation time constant τBK was set to 3 ms. The steady-state activation was given by [29]: (13) with (14) As BK channels are often co-localized with high-voltage activated Ca2+ channels, BK activation was assumed to depend on a domain concentration cdom in Ca2+ nanodomains, which in turn was assumed to be proportional to the instantaneous Ca2+ influx through ICa. We therefore set: (15) where A = 1.21 mmol⋅cm−1⋅C−1 is a parameter that converts a current density into a concentration, and cref = 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: (16) with an instantaneous, Ca2+ dependent, steady-state activation: (17) where [Ca] denotes the cytosolic Ca2+ concentration, and ks is a half-activation concentration of 0.4 μM.

ICa and ISK were dependent on the global cytosolic Ca2+ concentration. This was modelled as a simple extrusion mechanism, receiving a source through ICa, and with a concentration dependent decay term assumed to capture the effects of various ion pumps and buffering mechanisms: (18) Here, fc = 0.01 is the assumed fraction of free Ca2+ in the cytoplasm, α = 0.015mM ⋅ cm2/μC converts an incoming current to a molar concentration, and kc = 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 (INa) 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.

thumbnail
Download:
Fig 7. Experimentally recorded Na+ currents.

(A) Na+ currents evoked by the activation protocol. (B) Na+ currents evoked by the inactivation protocol. (C) Na+ currents used to determine the time constant for recovery from inactivation. (A-C) Voltage protocols are shown below the recorded currents, and all panels show a series of experiments (traces). 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

Steady-state activation and inactivation.

To determine steady-state activation, the peak current (Imax) was determined for each holding potential in Fig 7A, and the maximum peak was observed at about −10 mV. For inactivation, the peak current (Imax) was recorded for each holding potential in Fig 7B. In both cases, the maximal conductance (gmax) for each holding potential was computed by the equation: (19) 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 (ENa = RT/F ⋅ ln([Na]ex)/ln([Na]in) of 92 mV. The estimates of gmax 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 gmax on Vhold was fitted by a Boltzmann curve: (20) where corresponds to gmax 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), V* equals V1/2, i.e. the voltage where fbz has reached its half-maximum value. With a higher number of gates, V* = V1/2 + k ⋅ ln(21/a − 1). Eq 20 gave a good fit for the steady state activation (Fig 8A) and the steady state inactivation hinf (Fig 8B) with the parameter values for a, k and V* listed in Table 2.

thumbnail
Download:
Fig 8. Fitted kinetics for the Na+ current.

Voltage dependence of steady-state activation (A), steady-state inactivation (B), activation time constant (C) and inactivation time constant (D). The data points and curves in (A-B) were normalized so that activation/inactivation curves had a maximum value of 1 (assuming fully open channels). Dashed lines represent the same model when corrected for a liquid junction potential of 9 mV.

https://doi.org/10.1371/journal.pcbi.1006662.g008

thumbnail
Download:
Table 2. Parameters for Na+ activation.

The parameters p1-p6 are used together with Eq 22 to yield the time constants for steady state activation and inactivation (in units of ms). The remaining parameters are used together with Eq 20 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).

https://doi.org/10.1371/journal.pcbi.1006662.t002

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] (21) 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]: (22)

Good fits to the data points were obtained with the parameter values in Table 2.

Model for high-voltage activated Ca2+ channels

When estimating the steady-state values and time constant we followed procedures inspired from previous studies of L-type Ca2+ channel activation, we did not use Eq 6, but used the simpler kinetics scheme ICa = gHV Am2(V − ECa) (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 (ICa) was recorded (Fig 9A). Due to the small cellular size, perforated patches was used for recording the Ca2+ currents, and the recorded currents were small and noisy. As Fig 9A shows, the ICa responses did not follow a characteristic exponential curve towards steady state, as seen in many other experiments. Likely, this was due to ICa comprising a complex of different HVA channels (e.g., P, Q, R, L-type) which have different activation kinetics [49–51, 56–58]. In addition, some in some of the weaker responses ICa 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 Ca2+ 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).

thumbnail
Download:
Fig 9. Fitted kinetics for the Ca2+ current.

(A) ICa evoked by the activation protocol. (B) ICa evoked by the deactivation protocol. (A-B) Voltage protocols are shown below the recorded currents, and all panels show a series of experiments (traces). The current-traces were low-pass filtered with a cutoff-frequency of 300 Hz. (C) Voltage dependence of steady-state activation, normalized so that the activation curve had a maximum value of 1 (assuming fully open channels). (D) Activation time constant. Red data points were estimated from the deactivation protocol (B), while blue data points were estimated using the activation protocol (A). Dashed lines in (C-D) denote the kinetics scheme when corrected for a liquid junction potential of −15 mV for the experimental conditions used for recording Ca2+ currents.

https://doi.org/10.1371/journal.pcbi.1006662.g009

The peak current (Imax) was recorded for each holding potential in Fig 9A, and the maximum peak was observed at about 20 mV. By observations, (ICa) became an outward current when step potentials were increased beyond 70 mV, and based on this we assumed a reversal potential of ECa = 70 mV. Similar to what we did for the Na+ channel, the maximal conductance (gmax) for each holding potential was computed by the equation: (23) The estimates of gmax for activation are indicated by the crosses in Fig 9C. The dependency of gmax on Vhold 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.

thumbnail
Download:
Table 3. Parameters for Ca2+ 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).

https://doi.org/10.1371/journal.pcbi.1006662.t003

Time constants for ICa activation.

Assuming two opening gates (m), the time constant for ICa-activation was derived by fitting the function [47]: (24) 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.

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 (Cm), the models from Ileak, IK and ISK, and the reversal potentials for ENa and EK 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, Cm and the conductances (gX) 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 (25) which was selected so that the final capacitance per membrane area equal to 1μF/cm2. The same scaling was used when adopting IBK from [30]. Furthermore, the internal calcium handling (Eq 18) was also taken from the rat lactotroph model. However, to preserve the internal Ca2+ dynamics under the scaling, the factors α (in Eq 18) and A (in Eq 15), which converted Ca2+ currents into Ca2+ concentrations were multiplied by Fscale. A similar scaling was used in [37].

Second, INa and ICa 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 20 ms used for two different BK channels in the original study [30], and the same value as that used in [32], where a similar model for IBK was used. (ii) τK was set to 5 ms, which was faster than in the rat lactotroph model, where τK was based on experimental data where IK activation had a fast (3.7 ms) and a slow (30 ms) component [63]. The rat lactotroph model used the time constant of the slow component, while we chose a value closer to that of the fast component, as this was due to the faster INa-mediated APs seen in medaka gonadotroph. (iii) Through manual trial and error, the conductances of all ion channels were tuned freely in order to adapt the medaka gonadotroph model to the experimental data, both in terms of the control conditions and under application of paxilline. The final set of conductances were those given in Table 1.

References

  1. 1. Stojilkovic SS, Zemkova H, Van Goor F. Biophysical basis of pituitary cell type-specific Ca2+ signaling-secretion coupling; 2005.
  2. 2. Kidokoro Y. Spontaneous calcium action potentials in a clonal pituitary cell line and their relationship to prolactin secretion. Nature. 1975;258:741. pmid:813152
  3. 3. Li YX, Rinzel J, Vergara L, Stojilković SS. Spontaneous electrical and calcium oscillations in unstimulated pituitary gonadotrophs. Biophysical Journal. 1995;69(3):785–795. pmid:8519979
  4. 4. Van Goor F, Goldberg JI, Chang JP. Electrical membrane properties and ionic currents in cultured goldfish gonadotrophs. Can J Physiol Pharmacol. 1996;74(6):729–743. pmid:8909786
  5. 5. Levavi-Sivan B, Bloch CL, Gutnick MJ, Fleidervish IA. Electrotonic coupling in the anterior pituitary of a teleost fish. Endocrinology. 2005;146(3):1048–1052. pmid:15604206
  6. 6. Xu SH, Cooke IM. Voltage-gated currents of tilapia prolactin cells. General and Comparative Endocrinology. 2007;150(2):219–232. pmid:17045992
  7. 7. Haug TM, Hodne K, Weltzien FA, Sand O. Electrophysiological Properties of Pituitary Cells in Primary Culture from Atlantic Cod (Gadus morhua). Neuroendocrinology. 2007;86(1):38–47. pmid:17565196
  8. 8. Stojilkovic SS, Tabak J, Bertram R. Ion channels and signaling in the pituitary gland; 2010.
  9. 9. Tabak J, Tomaiuolo M, Gonzalez-Iglesias aE, Milescu LS, Bertram R. Fast-Activating Voltage- and Calcium-Dependent Potassium (BK) Conductance Promotes Bursting in Pituitary Cells: A Dynamic Clamp Study. Journal of Neuroscience. 2011;31(46):16855–16863. pmid:22090511
  10. 10. Hodne K, Strandabø RAU, von Krogh K, Nourizadeh-Lillabadi R, Sand O, Weltzien FA, et al. Electrophysiological Differences Between fshb- and lhb-Expressing Gonadotropes in Primary Culture. Endocrinology. 2013;154(9):3319–3330. pmid:23836032
  11. 11. Stojilkovic SS. Molecular mechanisms of pituitary endocrine cell calcium handling. Cell Calcium. 2012;51(3):212—221.
  12. 12. Strandabø RAU, Hodne K, Ager-Wick E, Sand O, Weltzien FA, Haug TM. Signal transduction involved in GnRH2-stimulation of identified LH-producing gonadotropes from lhb-GFP transgenic medaka (Oryzias latipes). Molecular and Cellular Endocrinology. 2013;372(1-2):128–139.
  13. 13. Chen C, Israel Jm, Vincent Jd. Electrophysiological Responses of Rat Pituitary Cells in Somatotroph-Enriched Primary Cultures. J Physiol. 1989;408:493–510.
  14. 14. Kukuljan M, Stojilković SS, Rojas E, Catt KJ. Apamin-sensitive potassium channels mediate agonist-induced oscillations of membrane potential in pituitary gonadotrophs. FEBS Letters. 1992;301(1):19–22. pmid:1333410
  15. 15. Tse A, Hille B. Role of voltage-gated Na+ and Ca2+ channels in gonadotropin-releasing hormone-induced membrane potential es in identified rat gonadotropes. Endocrinology. 1993;132(4):1475–1481.
  16. 16. Van Goor F, Zivadinovic D, Stojilkovic SS. Differential expression of ionic channels in rat anterior pituitary cells. Mol Endocrinol. 2001;15(7):1222–1236. pmid:11435620
  17. 17. Biales B, Dichter MA, Tischler A. Sodium and calcium action potential in pituitary cells. Nature. 1977;267:172. pmid:16073436
  18. 18. Chen C, Zhang J, Vincent JD, Israel JM. Sodium and calcium currents in action potentials of rat somatotrophs: Their possible functions in growth hormone secretion. Life Sciences. 1990;46(14):983–989. pmid:2157930
  19. 19. Zemkova H, Kucka M, Tomić M, Stojilkovic SS, Aguilera G. Spontaneous and CRH-Induced Excitability and Calcium Signaling in Mice Corticotrophs Involves Sodium, Calcium, and Cation-Conducting Channels. Endocrinology. 2016;157(4):1576–1589. pmid:26901094
  20. 20. Yu Y, Ali DW, Chang JP. Characterization of ionic currents and electrophysiological properties of goldfish somatotropes in primary culture. General and Comparative Endocrinology. 2010;169(3):231–243. pmid:20850441
  21. 21. Chang JP, Habibi HR, Yu Y, Moussavi M, Grey CL, Pemberton JG. Calcium and other signalling pathways in neuroendocrine regulation of somatotroph functions. Cell Calcium. 2012;51(3-4):240–252. pmid:22137240
  22. 22. Stojilkovic SS, Bjelobaba I, Zemkova H. Ion channels of pituitary gonadotrophs and their roles in signaling and secretion. Frontiers in Endocrinology. 2017;8(JUN):1–9.
  23. 23. Li YX, Stojilković SS, Keizer J, Rinzel J. Sensing and refilling calcium stores in an excitable cell. Biophysical Journal. 1997;72(3):1080–1091. pmid:9138557
  24. 24. LeBeau AP, Robson AB, McKinnon AE, Donald RA, Sneyd J. Generation of action potentials in a mathematical model of corticotrophs. Biophysical Journal. 1997;73(3):1263—1275. pmid:9284294
  25. 25. Van Goor F, Li YX, Stojilkovic SS. Paradoxical role of large-conductance calcium-activated K+ (BK) channels in controlling action potential-driven Ca2+ entry in anterior pituitary cells. J Neurosci. 2001;21(16):5902–5915. pmid:11487613
  26. 26. Tabak J, Toporikova N, Freeman ME, Bertram R. Low dose of dopamine may stimulate prolactin secretion by increasing fast potassium currents. J Comput Neurosci. 2007;22(2):211–222. pmid:17058022
  27. 27. Tsaneva-Atanasova K, Sherman A, van Goor F, Stojilkovic SS. Mechanism of Spontaneous and Receptor-Controlled Electrical Activity in Pituitary Somatotrophs: Experiments and Theory. Journal of Neurophysiology. 2007;98(1):131–144. pmid:17493919
  28. 28. Toporikova N, Tabak J, Freeman ME, Bertram R. A-Type K+ Current Can Act as a Trigger for Bursting in the Absence of a Slow Variable. Neural Computation. 2008;20(2):436–451. pmid:18047413
  29. 29. Duncan PJ, Engül S, Tabak J, Ruth P, Bertram R, Shipston MJ. Large conductance Ca 2+ -activated K + (BK) channels promote secretagogue-induced transition from spiking to bursting in murine anterior pituitary corticotrophs. J Physiol. 2015;5935(5935):1197–1211.
  30. 30. Duncan PJ, Shipston MJ, Tabak J, Ruth P, Bertram R. Glucocorticoids Inhibit CRH/AVP-Evoked Bursting Activity of Male Murine Anterior Pituitary Corticotrophs. Endocrinology. 2016;157(8):3108–3121. pmid:27254001
  31. 31. Fletcher PA, Zemkova H, Stojilkovic SS, Sherman A. Modeling the diversity of spontaneous and agonist-induced electrical activity in anterior pituitary corticotrophs. Journal of Neurophysiology. 2017;117(6):2298–2311. pmid:28228586
  32. 32. Fletcher PA, Sherman A, Stojilkovic SS. Common and diverse elements of ion channels and receptors underlying electrical activity in endocrine pituitary cells. Molecular and Cellular Endocrinology. 2018;463:23–36. pmid:28652171
  33. 33. Adler M, Wong BS, Sabol SL, Busis N, Jackson MB, Weight FF. Action potentials and membrane ion channels in clonal anterior pituitary cells. Proceedings of the National Academy of Sciences. 1983;80(7):2086–2090.
  34. 34. Cobbett P, Ingram CD, Mason WT. Sodium and potassium currents involved in action potential propagation in normal bovine lactotrophs. The Journal of Physiology. 1987;392:273–299. pmid:2451724
  35. 35. Mason WT, Sikdar SK. Characterization of Voltage-Gated Sodium Channels in Ovine Gonadotrophs: Relationship to Hormone Secretion. Journal of Physiology. 1988;399:493–517. pmid:2457092
  36. 36. Tennøe S, Halnes G, Einevoll GT. Uncertainpy: A Python Toolbox for Uncertainty Quantification and Sensitivity Analysis in Computational Neuroscience. Frontiers in Neuroinformatics. 2018;12:49. pmid:30154710
  37. 37. Tennøe S, Hodne K, Haug TM, Weltzien FA, Einevoll GT, Halnes G. [Re] Fast-Activating Voltage- and Calcium-Dependent Potassium (BK) Conductance Promotes Bursting in Pituitary Cells: A Dynamic Clamp Study. ReScience. 2019;5(1):1–15.
  38. 38. Toporikova N, Tabak J, Freeman ME, Bertram R. A-Type K+ Current Can Act as a Trigger for Bursting in the Absence of a Slow Variable. Neural Computation. 2007;20(2):436–451.
  39. 39. Fletcher P, Bertram R, Tabak J. From global to local: exploring the relationship between parameters and behaviors in models of electrical excitability. Journal of Computational Neuroscience. 2016;40(3):331–345. pmid:27033230
  40. 40. Lang DG, Ritchie AK. Tetraethylammonium blockade of apamin-sensitive and insensitive Ca2(+)-activated K+ channels in a pituitary cell line. The Journal of Physiology. 1990;425(1):117–132.
  41. 41. Sah P, Faber ESL. Channels underlying neuronal calcium-activated potassium currents. Prog Neurobiol. 2002;66(5):345–353. pmid:12015199
  42. 42. Faber ESL, Sah P. Calcium-Activated Potassium Channels: Multiple Contributions to Neuronal Function. The Neuroscientist. 2003;9(3):181–194. pmid:15065814
  43. 43. Miranda P, de la Peña P, Gómez-Varela D, Barros F. Role of BK potassium channels shaping action potentials and the associated [Ca(2+)](i) oscillations in GH(3) rat anterior pituitary cells. Neuroendocrinology. 2003;77 3:162–76. pmid:12673050
  44. 44. Gu N, Vervaeke K, Storm JF. BK potassium channels facilitate high-frequency firing and cause early spike frequency adaptation in rat CA1 hippocampal pyramidal cells. 2007;3:859–882.
  45. 45. Latorre E, Harries LW. Splicing regulatory factors, ageing and age-related disease. Ageing Research Reviews. 2017;36:165—170. https://doi.org/10.1016/j.arr.2017.04.004.pmid:28456680
  46. 46. Fontaine R, Hodne K, Weltzien FA. Healthy Brain-pituitary Slices for Electrophysiological Investigations of Pituitary Cells in Teleost Fish. JoVE. 2018;(138):e57790.
  47. 47. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol. 1952;117:500–544. pmid:12991237
  48. 48. Hodgkin AL, Katz B. The effect of sodium ions on the electrical activity of the giant axon of the squid. The Journal of Physiology. 1949;108(1):37–77.
  49. 49. Kay AR, Wong RK. Calcium current activation kinetics in isolated pyramidal neurones of the Ca1 region of the mature guinea-pig hippocampus. J Physiol. 1987;392:603–616. pmid:2451732
  50. 50. Borst JGG, Sakmann B. Calcium current during a single action potential in a large presynaptic terminal of the rat brainstem. Journal of Physiology. 1998;506(1):143–157. pmid:9481678
  51. 51. Jaffe DB, Ross WN, Lisman JE, Lasser-Ross N, Miyakawa H, Johnston D. A model for dendritic Ca2+ accumulation in hippocampal pyramidal neurons based on fluorescence imaging measurements. J Neurophysiol. 1994;71(3):1065–1077. pmid:8201402
  52. 52. Halnes G, Augustinaite S, Heggelund P, Einevoll GT, Migliore M. A multi-compartment model for interneurons in the dorsal lateral geniculate nucleus. PLoS Comput Biol. 2011;7(9):e1002160. pmid:21980270
  53. 53. Martina M, Jonas P. Functional differences in {Na+} channel gating between fast-spiking interneurons and principal neurones of rat hippocampus. J Physiol. 1997;505:593–603. pmid:9457638
  54. 54. Reckziegel G, Beck H, Schramm J, Elger CE, Urban BW. Electrophysiological characterization of Na+ currents in acutely isolated human hippocampal dentate granule cells. J Physiol. 1998;509.1:139–150. pmid:9547388
  55. 55. Traub RD, Miles R. Neuronal Networks of the Hippocampus. New York, NY, USA: Cambridge University Press; 1991.
    • 56. Migliore M, Cook EP, Jaffe DB, Turner DA, Johnston D. Computer simulations of morphologically reconstructed CA3 hippocampal neurons. J Neurophysiol. 1995;73(3):1157–1168. pmid:7608762
    • 57. Wolf JA. NMDA/AMPA Ratio Impacts State Transitions and Entrainment to Oscillations in a Computational Model of the Nucleus Accumbens Medium Spiny Projection Neuron. Journal of Neuroscience. 2005;25(40):9080–9095. pmid:16207867
    • 58. Lindroos R, Dorst MC, Du K, Filipović M, Keller D, Ketzef M, et al. Basal Ganglia Neuromodulation Over Multiple Temporal and Structural Scales-Simulations of Direct Pathway MSNs Investigate the Fast Onset of Dopaminergic Effects and Predict the Role of Kv4.2. Frontiers in Neural Circuits. 2018;12(February):1–23.
    • 59. Hines ML, Davison AP, Muller E. NEURON and Python. Frontiers in neuroinformatics. 2009;3(January):1. pmid:19198661
    • 60. Ermentrout B. Simulating, analyzing, and animating dynamical systems: a guide to XPPAUT for researchers and students. vol. 14. Siam; 2002.
      • 61. Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, et al. Global Sensitivity Analysis. The Primer. Wiley; 2007. Available from: http://doi.wiley.com/10.1002/9780470725184.ch6{%}5Cnhttp://doi.wiley.com/10.1002/9780470725184.
        • 62. Halnes G, Ulfhielm E, Eklöf Ljunggren E, Kotaleski JH, Rospars JP. Modelling and sensitivity analysis of the reactions involving receptor, G-protein and effector in vertebrate olfactory receptor neurons. Journal of Computational Neuroscience. 2009;27(3):471–491. pmid:19533315
        • 63. Herrington J, Lingle CJ. Multiple components of voltage-dependent potassium current in normal rat anterior-pituitary-cells. Journal Of Neurophysiology. 1994;72(2):719–729.