Original Research ARTICLE
Modeling the Influence of Ion Channels on Neuron Dynamics in Drosophila
- 1School of Life Sciences, Arizona State University, Tempe, AZ, USA
- 2School of Mathematical and Statistical Sciences, Arizona State University, Tempe, AZ, USA
Voltage gated ion channels play a major role in determining a neuron's firing behavior, resulting in the specific processing of synaptic input patterns. Drosophila and other invertebrates provide valuable model systems for investigating ion channel kinetics and their impact on firing properties. Despite the increasing importance of Drosophila as a model system, few computational models of its ion channel kinetics have been developed. In this study, experimentally observed biophysical properties of voltage gated ion channels from the fruitfly Drosophila melanogaster are used to develop a minimal, conductance based neuron model. We investigate the impact of the densities of these channels on the excitability of the model neuron. Changing the channel densities reproduces different in situ observed firing patterns and induces a switch from integrator to resonator properties. Further, we analyze the preference to input frequency and how it depends on the channel densities and the resulting bifurcation type the system undergoes. An extension to a three dimensional model demonstrates that the inactivation kinetics of the sodium channels play an important role, allowing for firing patterns with a delayed first spike and subsequent high frequency firing as often observed in invertebrates, without altering the kinetics of the delayed rectifier current.
It is well-known that different neuron types exhibit distinct characteristic features under standardized or similar conditions such as constant current injection, due in part to the influence of differing ion channel kinetics and distributions (Shepherd, 2004). Experimental and theoretical studies show that differences in spiking patterns can be related to different combinations of ion channel densities (Goldman et al., 2001; Zeberg et al., 2010) with different channel kinetics. These differences affect the timing of action potentials (APs) and influence subthreshold integration of synaptic input and the filtering properties of the neuronal structure, resulting in bandpass or highpass filtering properties. Consequently, the response of a neuron to synaptic input depends on the underlying dynamics of membrane excitability. Hodgkin classified neurons according to their spiking behavior upon steady current injection, and the resulting frequency–current relationships (f –I curves) generally can be divided into three distinct classes (Hodgkin, 1948). Numerous subsequent studies have analyzed the relationship between this classification and the output properties of a model neuron (Ermentrout, 1996; Rinzel and Ermentrout, 1998; Gutkin et al., 2003; St-Hilaire and Longtin, 2004; Tateno et al., 2004; Tateno and Robinson, 2006, 2007) and, from a dynamical systems point of view, how this classification relates to the underlying mathematical structure of the model (Izhikevich, 2007; Prescott et al., 2008).
In spite of these theoretical studies, the impact of specific ion channel kinetics on neuronal function remains largely unclear. Drosophila provides a valuable model system for investigating ion channel kinetics and their impact on firing properties. Neurons can be identified individually, and many molecular mechanisms are comparable to those in vertebrate systems. The spiking responses and neuronal morphology in these neurons have been investigated at different developmental stages (Choi et al., 2004) and changes in these properties have been observed during development (Duch and Levine, 2000), after targeted genetic manipulations, and under different pharmacological conditions (Peng and Wu, 2007; Duch et al., 2008; Ryglewski and Duch, 2009). In particular, the use of Drosophila in research is of great interest due to the development of new genetic tools for experimentation. Despite the increasing importance of Drosophila as a model system (reviewed by Baines and Pym, 2006 and Corty et al., 2009), few computational models of its ion channels have been developed. The use of computational modeling techniques can help predict the behavior of membrane dynamics at experimentally inaccessible locations and help connect electrophysiological and other molecular biological findings to neuronal function. Here, we create mathematical models, based on experimental data from Drosophila, in order to conduct a computational study that investigates how changes of different parameters in physiological ranges influence the input-output properties of neurons. We focus on the identified Drosophila motoneuron 5 (MN5) since its morphology, electrophysiology and certain aspects of its behavior during flight have been well-characterized experimentally.
The generation of action potentials, along with their shape and firing patterns, depends in large part on voltage gated sodium (Na+) and potassium (K+) channels. Drosophila has only one confirmed Na+ channel gene DmNav (Miyazaki et al., 1996; Mee et al., 2004), which is subject to alternative splicing. The voltage dependence of the macroscopic currents carried by the different splice variants has been characterized in heterologous expression systems using voltage clamp recordings (Olson et al., 2008; Lin et al., 2009). K+ channels show the greatest diversity among ion channels (Jan et al., 1977; Coetzee et al., 1999), where voltage gated ion channels fall roughly into two categories, the non-inactivating or slowly inactivating delayed rectifier and the rapidly inactivating, transient A-type currents (Hille, 1992). In Drosophila neurons two genes, Shab and Shaw, encode delayed rectifier current conducting channels. Shab and Shaw channels are members of the Kv2 and Kv3 subfamilies, respectively (Wei et al., 1990; Covarrubias et al., 1991; Tsunoda and Salkoff, 1995a,b). Shaw channels demonstrate low voltage sensitivity, suggesting that they operate as leak channels, while Shab channels conduct the majority of delayed rectifier currents (Tsunoda and Salkoff, 1995b). The kinetics of the Shab channel are reported to be comparable to the classical (as described by Hodgkin–Huxley) delayed rectifier K+ channel (Tsunoda and Salkoff, 1995b). Two further genes encode channels conducting A-type currents that show fast inactivation kinetics. Activation and inactivation properties of Drosophila voltage gated K+ channels have been characterized in homologous expression systems as well as in Drosophila neurons (Covarrubias et al., 1991; Islas and Sigworth, 1999; Tsunoda and Salkoff, 1995b; Gasque et al., 2005). Their contributions to firing properties have been studied using pharmacology and genetical manipulations in order to remove the currents (Choi et al., 2004; Gasque et al., 2005; Peng and Wu, 2007; Ryglewski and Duch, 2009; Ping et al., 2011). Comparing different mutant neurons indicates that Shab is required for repetitive spiking, while A-type channels regulate the firing frequency and lower the spike threshold to induce a delay to first spike (Choi et al., 2004; Ping et al., 2011).
However, accessing the specific role of different ion channels is challenged by the considerable amount of variability across individuals. Drosophila MN5 shows considerable animal to animal variation in spiking behaviors upon current injection (Duch et al., 2008), ranging from non-repetitive responses (single spikes or single graded responses) to repetitive spiking with different f –I curves and varying delays to first spike. Data from intracellular recordings in whole cell patch configuration (shown in Figure 1) demonstrate different firing behaviors seen in Drosophila MN5. Some cells exhibit a single AP with no delay to first spike for small Iapp or display only a slight increase in firing frequency with increasing Iapp. This indicates a discontinuous f –I curve and relates to a dynamical system close to a Hopf bifurcation (type II dynamics). Other cells have a comparably lower firing threshold, where smaller Iapp-values induce low frequency spiking. In this case a continuous f –I curve is possible, which can be related to a saddle node on invariant cycle (SNIC) bifurcation (type I dynamics). In rare cases, repetitive spiking occurs after a long delay, where the interspike intervall (ISI) is smaller than the initial delay and higher amplitudes elicit higher frequency spiking without initial delay (see Herrera-Valdez et al., 2013). In dynamical systems this type of behavior is observed if the system is close to a saddle small homoclinic bifurcation (Izhikevich, 2007).
Figure 1. Different firing behaviors from three Drosophila MN5s. Intracellular recordings in whole cell patch configuration carried out by S. Ryglewski in the Duch laboratory reveal substantial animal to animal variation in the responses to current pulses (n = 52). Cell 1 (38 %): A single action potential for a low-amplitude current stimulation, with repetitive spiking that adapts in frequency for slightly increased stimulation. Firing frequency increases with stimulation amplitude. Cell 2 (8 %): A very low-amplitude stimulation produces repetitive spiking with a relatively large duration of single APs. Stimulations with an amplitude that is close to the firing threshold observed in the other cases (0.4 nA) induces a large initial spike with subsequent broader spikes that diminish in amplitude. Increasing the stimulation amplitude further results in dampening oscillations that finalize in a depolarization block. Cell 3 (42 %): Repetitive firing is induced by a stimulation amplitude of 0.4 nA, while the cell remains quiescent with stimulation amplitudes of 0.3 nA. Higher stimulus amplitudes result in only a slight increase in firing frequency (42 %). Experimental methods provided by Ryglewski and Duch (2009).
In addition to animal to animal variability, compensation mechanisms and modulation adjust the membrane properties of neurons, and experimental studies using genetic knock downs can be compromised by these homeostatic regulation mechanisms (Marder, 2011). Here we determine whether this known variability can be captured sufficiently by a minimal model with only two realistic channel types. This provides a foundation for understanding the roles of these channels and for future studies with additional channel types where detailed mathematical analysis is not possible due to the large number of variables. We begin with a model of a patch of excitable membrane where we carefully develop models of Na+ and K+ voltage gated ion channels displaying kinetics based on channels expressed in Drosophila. The Na+ channel is based on kinetics of DmNaV29, which mediates a fast inward current, and the K+ channel is based on kinetics of Shab, which mediates a delayed rectifier outward current. In contrast to the electro-diffusion based model in our previous work (Herrera-Valdez et al., 2013), here we use a conductance-based model, where the parameters of the model are adjusted to better resemble published data. First, using a mathematically reduced two dimensional model, we ask whether changing the density of Shab channels is sufficient to switch between bifurcation types and whether this can cause the different observed response properties. Employing sinusoidal current stimulation, we confirm the occurrence of qualitatively different responses to rhythmic synaptic input. Next, we extend the model to three dimensions and find that adjusting the inactivation kinetics of Na+ promotes a saddle small homoclinic bifurcation leading to a delayed first spike, without changing the kinetics of the Shab channels.
2. Materials and Methods
2.1. Conductance Based Membrane Model
Membrane potential (Vm) dynamics are described with a Hodgkin–Huxley type conductance based model (Hodgkin and Huxley, 1952), where the current balance equation takes the form
Here cm is the membrane capacitance, Iapp the applied current, and Ii is the ionic current mediated by channel type i. The value of cm was determined using the whole cell capacitance (≈1.3nF) as measured in voltage clamp (Ryglewski and Duch, 2009) and the surface area of MN5 (≈10, 000μm2) as determined by geometric reconstruction from confocal image stacks (Vonhoff and Duch, 2010). Ii are modeled according to Ohm's law:
where gi is the voltage dependent conductance and Vi is the reversal potential of channel i. The general formulation for gi includes two gating variables representing activation and inactivation kinetics, p and q, respectively. Each gating variable is associated with a number of independently operating gates j and k, respectively; such that
where is the maximal conductance of channel type i. The open probability of a gate p is described with
where p∞ is the steady state function and τp is the voltage dependent change of the time constant. The description of gate q takes the same form. By assuming that opening or closing of a gate results from charged particles moving across the membrane in response to an electric field, p∞ and τp become
where Jk−1 is the Boltzmann constant, T = 295.15 K the temperature, the half activation voltage, rp the rate of activation, and γp the relative position of the energy barrier for a gating particle in the membrane. zpe is the gating charge, where e ≈ 1.60217733 × 10−19C is the elementary charge and zp reflects the amount and distance the gating particle is moved (Willms et al., 1999; Destexhe and Huguenard, 2000). The parameter values are summarized in Table 1.
2.2. Minimal Ion Channel Model
For the minimal model, one Na+ current, one K+ current, and one leak current determine the intrinsic dynamics of Vm. The voltage dependence of the Na+ current (INa) is based on currents mediated by channels encoded by the splice variant DmNav29, which is among others expressed in adult Drosophila neurons. The K+ current (Isb) is a delayed rectifier current based on kinetics reported for Shab channels, where we neglect the slow inactivation of those channels to aid with mathematical analysis. The leak current, ILeak, accounts for currents with relatively small voltage dependence. Generally, the membrane currents for Na+ and K+ and the leak current take the form
however, we consider several variations on this model in order determine and understand the range of cell dynamic behaviors while fitting the parameters to experimental data.
To begin, the fast activation is assumed to be at steady state, and we assume that the inactivation kinetics can be represented as a function of the activation kinetics of the Shab channel gating (Rinzel, 1985; Av-Ron et al., 1991) in order to reduce the model to two dimensions. Further discussion of these assumptions is provided below. This reduced system takes the form
allowing for phase plane analysis. Voltage clamp recordings are used to constrain the parameters. When the voltage is held constant at Vcl, the solution to Equation (4) is
therefore, the membrane current mediated by a channel denoted by the index i becomes
To model the activation gating (m) of the DmNaV29 channel, we use previously published parameters (Herrera-Valdez et al., 2010). Recall that, initially, the parameters for the inactivation gating (h) are the same as the parameters for the activation gating (b) of the Shab channels. The activation and inactivation curves as reported by Olson et al. (2008) were used to evaluate the channel model. The experimental Na+ activation curve was obtained by clamping the cell at a holding potential of −120 mV and recording the current in response to step potentials from −120 to 60 mV in 5 mV increments. The peak currents for each step potential are divided by (Vm − VNa) to obtain the peak conductance. In experiments, Olson et al. (2008) obtained the steady state inactivation curve for Na+ using pre-steps from −120 to 40 mV followed by a step potential of −5 mV. The pre-step potential is assumed to be sufficiently low so that essentially all inactivation gates (h) are open and all activation gates (m) are closed. The duration of the pre-step is assumed to be long enough that the system is at steady state. The peak current of each trace is divided by the maximal peak current among all steps, which is the one in response to the most hyperpolarized pre-step and assumed to be the current amplitude, given all channels are in the open state. We simulate those experiments using Equation (13) for comparison. Subsequently, we incorporate independent Na+ inactivation by restoring the third gating variable, h using parameters reported for inactivation. Again we simulate the experiments using Equation (13) for comparison of the activation and inactivation curves.
The parameters for Isb were constrained using voltage clamp recordings from Drosophila embryonal cell cultures (Tsunoda and Salkoff, 1995b). Cells were clamped to a holding potential of −50 mV, and the current in response to step potentials from −20 to 50 mV in 10 mV increments was recorded. The published traces were digitized and fitted simultaneously to Equation (13) using a non-linear least-squares optimization algorithm. Different initial conditions were used, that converged all to the same values.
Maximal conductances and can be seen as the combination of the amount of channels in the membrane and the maximal conductance of the corresponding single channels. These values were also fitted but since the number of channels can be very different in different neurons, in what follows we focus on the relative amounts of these values.
Bifurcation diagrams were generated with numerical continuation methods, using PyCont a sub-package of PyDSTool (Clewley et al., 2007), which provides an interface to AUTO (Doedel et al., 2000). Phase response curves (PRCs) were calculated with the adjoint method implemented in PyDSTool. Stable and unstable manifolds for saddle points were obtained using XPPAUT (Ermentrout, 2002) with a time step of 0.001 ms. The phase responses were determined by calculating the system's adjoint (Ermentrout and Kopell, 1991) using PyDSTool. Numerical simulations were performed using Python 2.7. The systems of ordinary differential equations are integrated using odeint with default settings from the scipy package, which uses lsoda from the FORTRAN odepac. For nonstiff problems, AdamsĂŹ method was used, and for stiff problems, a method based on backward differentiation formulas was used. Results are compared to results obtained with XPPAUT using time steps of 0.05 and 0.01 ms to ensure accuracy. Fitting of the electrophysiological data was performed using leastsq from the scipy optimization package.
3.1. Fitting of Ion Channel Models
To obtain models for Shab channels, traces from voltage clamp recordings (Tsunoda and Salkoff, 1995b) are fit with Equation (13) as described above. The number of gates are usually chosen to best represent the experimentally measured ion current dynamics. Since the inactivation of Shab channels is small, we set the number of inactivation gates to zero as shown in the Section 2 above. For the activation gating (b) we use a power of four as used by Hodgkin and Huxley (1952), but also investigate a power of one as in Herrera-Valdez et al. (2013). Both powers provide fits that are in good agreement with the voltage clamp recordings of Shab currents (Figure 2A) although they require different values for other parameters in order to best fit the data. The values of the different parameter sets are summarized in Table 2. In Equation (13) the electromotive force is based on Ohms's law, as it is commonly described in Hodgkin–Huxley type conductance based models and displays a linear current–voltage relationship.
Figure 2. Data fits of Shab and DmNav29 channel parameters. (A) Digitized points from Figure 4A from Tsunoda and Salkoff (1995b) are marked as black circles. Model of Shab channels with parameters obtained from fits with power one (sb1), power four (sb4), hand tuned parameters (sbdef) and fits with power four and fixed zb (sbz). Compare with Figure 7 from Herrera-Valdez et al. (2013). (B) Steady state activation and inactivation curves of DmNav29 according to Olson et al. (2008) (gray) and from models (note that the activation kinetics are the same in all models). Results using parameter sets sb1, sb4, sbdef, and independent inactivation used in the three dimensional model.
For Na+ channels, only the parameters for curves used to fit the peak and tail currents from voltage clamp experiments are available (Olson et al., 2008). Therefore, only zm, zh, , and can be assessed. For the two dimensional reduced model where the inactivation curve for INa depends on the activation of the Shab channels, inactivation is represented best using the term (1−b) with a power of one for the Shab gating variable (b) (Table 2, Figure 2B) because the half activation of a single gate is close to the value reported for DmNav29 (Table 3). The half inactivation values for other DmNav splice variants are between −34.9 and −66.4 mV (Olson et al., 2008; Lin et al., 2009), where none of these reported splice variants has a half inactivation that is close to the value found for Shab channel activation using a power of one. In addition, the slope of the resulting inactivation curve for INa is too flat, indicating that the gating charge zb is too small. The Shab channel kinetics were adjusted in an attempt to improve zb, leading to parameter set sbdef; however, only a slight improvement could be achieved without compromising the Shab kinetics (Figure 2B). Due to these issues, after analyzing the two dimensional model, we extend the model to three dimensions by restoring the variable h for INa inactivation, which allows for different parameters for Na+ inactivation and K+ activation. Adjusting the gating parameter as summarized in Table 2 leads to a good agreement between model and data (Figure 2B) for the three dimensional system.
In the following we use sbdef as shown in Table 2 as a default parameter set to model Shab channels dynamics in order to relate the two dimensional analysis to the three dimensional results. As described below, the results for the three dimensional system agree with the characterization of the reduced two dimensional model.
3.2. General Behavior of the Two Dimensional Membrane Model
In the following we keep constant and vary the dimensionless ratio of and . This serves as a measure for the Shab channel density, assuming that the maximal conductance is the product of the amount of channels in the membrane and the conductance of a single channel in the open state.
With a = 1.1347, the model resembles the basic features of Drosophila MN5 (Figure 3A) in response to steady current injection. Current pulses with amplitudes of the same order of magnitude as used in experiments cause the model to fire repetitively. However, compared to the majority of available recordings, the AP duration is longer and the firing frequency is higher (Figure 3A). It is likely that these features can be adjusted by adding A-type currents to the model. Shaker current is thought to adjust the spike shape (Peng and Wu, 2007), and Shal current is thought to decrease the firing frequency (Tsunoda and Salkoff, 1995a). Yet in some experimental recordings, the AP duration is even longer than in the model (Figure 3B). In contrast to our earlier models (Herrera-Valdez et al., 2013), in this model APs appear on elevated potentials for a wide range of parameters, including the default parameter set sbdef. This is a consequence of the combination of slightly different values for γb, zb, and rb.
Figure 3. Comparison of electrophysiological data and model behavior. (A) Membrane potential in response to squared pulse current injections of different amplitudes in the model and three different cells. (B) Phase plots of electrophysiological recordings and the model in response to squared pulse current injection. Experimental methods provided by Ryglewski and Duch (2009).
In large part, the density of Na+ channels determines the maximal rate of change of the voltage (dVm∕dt) during the rising phase of an action potential. A phase-plane plot showing dVm ∕dt vs. Vm for three different recordings as well as the model is depicted in Figure 3B. The maximal value of dVm∕dt varies among cells and is between 55 and 20 mV/ms.
With a of 15.62 µS (Figure 3B right, black trace) dVm∕dt for the model takes a maximal value of 35 mV/ms, which is between the observed values. Increasing to 20 µS results in a maximal dVm∕dt of about 55 mV/ms during the rising phase, which corresponds to the largest observed values (Figure 3B, right, gray trace). In the following we show results for . However, a model with in combination with a smaller (0.01 µS for example) shows qualitatively similar behavior.
3.3. Shab Channel Density Determines Bifurcation Type
In order to compare the model response to published experimental data, we examine the behavior of Vm over time (Figure 4). Visual inspection reveals that changing has little influence on the spike shape. With , arbitrarily low firing frequencies with long delays to spike can be evoked, indicating type I dynamics. Current pulses with an amplitude of 0.8 nA elicit a single spike and subsequent low amplitude oscillations. Although the transition from rest to spiking occurs via different bifurcations, similar firing patterns are elicited with -values of 0.782 and 1.1347. In both cases a single spike can be elicited with small Iapp, while slightly increased stimulus amplitudes induce repetitive spiking with relatively high frequencies, indicating type II dynamics. Furthermore, the delay of spiking is small compared to observations when . In contrast to what was reported by Herrera-Valdez et al. (2013), the frequently observed pattern of a single spike with small and repetitive spiking with increased stimulus amplitudes in experiments can be generated with a system near a Hopf (Figure 4, bottom) as well as with a system near a saddle node bifurcation (Figure 4, middle). Note however that the range of current amplitudes where a single spike is elicited is narrow.
Figure 4. Firing behavior for different relative Shab densities in the two dimensional model. Responses to steady current injection (left) and f –I curves (right), from top to bottom : 0.6, 0.782, 1.1347 with Rin 70.2, 69.5, 65.4 and 59.4 MΩ. Note that changing the Shab densities affects Rin.
The two parameter bifurcation diagram (Figure 6B) displays the change of the locations of saddle-nodes and Hopf points in the Iapp - plane. Low levels of Shab density result in a saddle node bifurcation, and the saddle-nodes remain until is about 1.2216, where the two branches collide at a cusp point. However, near and Iapp = 0.4087 nA there is a Bogdanov–Takens bifurcation. At this point the saddle node curve meets a Hopf curve, and for higher values of , the fixed point loses stability at a Hopf point. This demonstrates that a non-monotonic I–V-curve does not necessarily indicate a loss of stability via a saddle node bifurcation. At the right branch of the Hopf curve crosses the right branch of the saddle node curve, meaning that the system has two stable fixed points for an increasing range of Iapp.
We use phase plane analysis to further investigate the characteristics of the system with three chosen -values, for which stable periodic solutions can be induced. The model is stimulated with constant current injections of different amplitudes. The amplitudes were chosen to be just below, at and just above the emergence of the limit cycle as determined by the bifurcation diagram of Figure 6. Figure 5 shows the phase plane indicating the fixed points and the stable and unstable manifolds of the neutral saddle, when they exist. Without or with small subthreshold Iapp, the unstable manifold connects the saddle node and the stable fix point to form a heteroclinic trajectory (Figures 5A,B, left). Letting , and injecting a current of 0.296 nA, the unstable saddle point and stable node coalesce and disappear. The trajectory becomes a homoclinic invariant circle and gives rise to a limit cycle with infinite period; that is, stability is lost via a SNIC bifurcation. When Iapp is increased further, low frequency spiking is induced (Figure 5A, right).
Figure 5. Changes of fixed point types and trajectories with increased constant current injection in phase space. 0.6 (A), 0.782 (B), and 1.1347 (C). Trajectories corresponding to initial conditions of the steady state without current injection are color coded from yellow (t = t0) to red (t = tmax). Fixed points, unstable and stable manifolds of the saddle from left to right for in increasing Iapp as indicated.
Increasing to 0.782, a single spike can be elicited with Iapp = 0.39 nA. At this value the unstable manifolds still forms a heteroclinic orbit; however, the shape of the stable manifold changes so that it separates the initial values at rest and the stable fixed point (Figure 5B, middle). Hence a single spike is elicited because the trajectory must take a large excursion in order to approach the equilibrium. At Iapp = 0.395 nA, one unstable manifold converges onto itself, and a limit cycle attractor appears (Figure 5B, right). The trajectory moves around the stable fixed point indicating a saddle big homoclinic bifurcation, where a stable and unstable manifold of a saddle point coincide and form a homoclinic orbit, and the two other branches lie inside the homoclinic orbit. The unstable manifold separates the stable fixed point and the initial conditions at rest, resulting in repetitive spiking. A saddle node bifurcation occurs for Iapp above 0.4 nA.
When increasing to 1.1347, stability is lost via a Hopf bifurcation (Figure 5C). With Iapp = 0.55 nA, a single spike is elicited (Figure 5C, left). A stable limit cycle appears with Iapp = 0.6 nA, while the fixed point remains stable (Figure 5C, middle). The fixed point becomes unstable when increasing Iapp (Figure 5C, right).
We find that changing results in different bifurcation types as the system undergoes a transition from rest to spiking in response to current injection (Figure 6) where low levels of Shab result in a saddle node bifurcation and high levels in a Hopf bifurcation. Figure 6A shows the bifurcation diagram for -values of 0.57, 0.6, 0.782, and 1.347. This diagram shows the local stability of the fixed points (black) of the system and indicates Hopf points (blue) and saddle nodes (green). In all cases there are saddle nodes, but for the stability changes via a subcritical Hopf bifurcation at Iapp = 0.653 nA. A stable limit cycle emerges at Iapp = 0.595 nA, hence a stable fixed point and stable limit cycle coexist. The coexistence of a stable fixed point and a stable limit cycle is also given with . Yet, the stability of the fixed point is lost via a saddle node bifurcation at Iapp = 0.4037 nA, while the stable limit cycle appears at Iapp = 0.3958 nA via a global fold bifurcation of limit cycles (Figure 6A, bottom left panel inset). Decreasing further to 0.6, the stable fixed point looses stability at Iapp = 0.4037 nA and a stable limit cycle emerges (Figure 6A, top right panel inset), indicating a saddle node on invariant cycle bifurcation and hence type I excitability. The saddle node bifurcation is also present with , yet no stable limit cycle could be found. Instead, an unstable limit cycle appears at about Iapp = 0.248 nA and remains unstable until it disappears at the Hopf bifurcation. The trajectories of the system will traverse the phase space toward the more depolarized stable fixed point, if the system is pushed beyond the bifurcation point.
Figure 6. Bifurcation structure for different relative Shab densities. (A) One parameter bifurcation diagrams for : 0.57, 0.6, 0.782, 1.1347 (black to gray lines). Increasing changes the transition from rest into spiking from a saddle node (green points) to a Hopf (blue points) bifurcation. Note that the stable branch continues outside of the shown region for negative values of Iapp. The minimum and maximum of periodic solutions are marked by red to orange lines for increasing . The lines style indicates the stability (solid: stable; dashed: unstable). (B) Two parameter bifurcation diagram showing the course of the Hopf (blue) and saddle node (green) curves in the Iapp x plane.
3.4. Responses to Time Varying Input
These different bifurcation types imply different phase response curves (PRCs) and different responses to periodic forcing. Therefore, the next step is to calculate the PRCs and analyze the model behavior for periodic current stimulation.
We use periodic current injections with the shape of a sine wave with varying periods to assess the model response to periodic forcing and to find the preferred frequency of the model. The frequency sweeps from 0 to 30 Hz linearly over a time interval of 10 s. As expected, systems near a saddle node bifurcation with of 0.6 or 0.782 displays integrator properties, behaving as a low-pass filter of the input current (Figure 7A, top and middle panel). In both cases the amplitude of Vm decreases with stimulus frequency. Current amplitudes above the respective bifurcation points are required to induce APs even at low frequencies. With larger amplitudes repetitive firing is induced at low frequencies and the spike count per cycle declines with increasing frequency. In contrast, with , the model exhibits resonator properties, acting as a band-pass filter of input current. Stimulation with low amplitude swept sine current injection reveals that frequencies around 15 Hz elicit APs, while higher and lower frequencies only cause graded responses (Figure 7A, bottom panel). Upon increasing the stimulus amplitude, the frequency band that yields spiking becomes larger. Sinusoidal current injections confirm that with an amplitude of 0.55 nA a frequency of 15 Hz, but not 14 or 16 Hz results in the generation of APs (Figure 7B). In summary, while the response properties to constant current injection with of 0.782 and 1.1347 are quite similar, the responses to periodic forcing are dramatically different.
Figure 7. Firing behavior in response to periodic forcing for different relative Shab densities in the two dimensional model. (A) Responses to swept sine current injection, with frequencies from 0 to 30 Hz. From top to bottom : 0.6, 0.782, 1.1347. (B) Responses of the model with = 1.1347 to sinusoidal forcing with three different frequencies.
How a rhythmically firing neuron responds to perturbations can be analyzed with the PRC, which provides a measure for the timing of the subsequent AP after a brief stimulus at a specific time during the spiking cycle. If a neuron spikes repetitively with a specific frequency with constant Iapp or without a driving force, the state of the neuron can be expressed by a single phase variable. To obtain the phase, the time of a peak voltage equates to a phase of zero and subsequent times are divided by the spiking period. The PRC gives the phase shift induced by a perturbation as a function of the phase at which the perturbation occurs. Positive values indicate a phase advance, which means the next AP will be sooner, while negative values indicate a phase delay, which means the next AP will be later than without perturbation. PRCs that are mostly positive are called type I, whereas PRCs that have a positive and a negative component are classified as type II. In order to classify PRCs as type I or type II, Tateno and Robinson (2007) introduced a r-value, the ratio of maximal phase delay and advance, of 0.175.
Figure 8 shows the PRCs of the model for the different indicated -values. For each parameterization of three different values for Iapp were used, which starts close to the firing threshold and is increased subsequently. With the PRC can be classified as type I, with r-values of 0.01, 0.04, and 0.15 for Iapp 0.297, 0.3, and 0.37nA. By increasing to 0.782, the PRC has small negative values for small phase values, but the r-values with 0.067, 0.064, and 0.14 for Iapp 0.397, 0.4, and 0.5 still leads to a classification as type I. With a further increase of to 1.1347, the negative values in the early phase increase, for small baseline Iapp of 0.596 and 0.6 nA the PRC is type II, with r-value of 0.192 and 0.179 bigger than 0.175. However, with a baseline Iapp of 0.7 nA the r-value is 0.130. As can be seen by the bifurcation diagram in Figure 6, although the stable fixed point loses stability via a Hopf bifurcation, there is also a saddle node bifurcation near Iapp = 0.7 nA, which might influence the vector field in a way that the model exhibits a type I PRC. For higher baseline Iapp the r-value increases and exhibits a PRC classified as type II (not shown). Also when increasing to 1.5 the PRC is classified type II for all probed stimulus amplitudes that elicit repetitive spiking (not shown).
Figure 8. Phase response curves for different relative Shab densities (indicated by colors) and from left to right for increasing Iapp, as indicated.
The PRC corresponding to all three -values have in common that they are relatively insensitive to perturbations at the beginning and the end of the cycle shortly before and after the peak of the AP. Further, the maximum of the PRC shifts to the right with higher baseline stimulation.
3.5. Small Parameter Changes Induce Delay to First Spike
The only firing profile that could not be reproduced by changing is an onset with a long delay and subsequent higher frequency firing. Usually a saddle small homoclinic bifurcation is responsible for such behavior (Izhikevich, 2007). We explored the parameter space further in order to find a regime where the system undergoes this bifurcation. We set zb to a constant value of 1.5 in order to obtain a better description of the Na+ channel inactivation (see Section 3.1). Then we performed the fitting routine, resulting in the parameter set sbz (see Table 2) with which the desired firing profile can be reproduced (Figure 9C). The bifurcation diagram for reveals that at Iapp ≈ 0.1392 nA the minimum on the periodic orbit collides with the saddle node, and the coexisting stable fixed point is at a smaller Vm (Figure 9A). By increasing , the system exhibits the same bifurcation types as with the default parameters sbdef. With this set of parameters must be much higher in order to produce the various firing patterns than with the previously used parameters.
Figure 9. Saddle small homoclinic bifurcation with changed Shab parameters. (A) Right: Bifurcation diagram shows, that with = 1.35 a stable limit cycle (dark red) coincides with the curve of saddle nodes (black dots), while the stable fixed point vanishes via a saddle node bifurcation at a bigger Iapp. Middle: Increasing (dark to lighter brown and orange for fixed points and periodic solutions) changes the transition from rest into spiking from a saddle node (green points) to a Hopf (blue points) bifurcation. Dashed and solid lines mark stability (solid: stable; dashed unstable). Left: Two parameter bifurcation diagram indicating the location of Hopf points (blue) and saddle nodes (green) as well as the Bogdanov–Takens (square) and cusp (triangle) points with respect to . (B) Changes of fixed point types and trajectories for different Iapp in phase space. From left to right Iapp increases as indicated. Trajectories corresponding to initial conditions of the steady state without current injection are color coded from yellow (t = t0) to red (t = tmax). (C) Left: Responses to constant current injection. The stimulus inducing the response shown in yellow exhibits an additional short current pulse (Is) with 30 ms and 0.01 nA at t = 150 ms. Right: f –I curve. (D) Responses to swept sine current injection, with frequencies from 0 to 30 Hz. (E) PRCs for increasing Iapp, as indicated.
Phase plane analysis confirms that the underlying bifurcation when is indeed a saddle small homoclinic bifurcation (Figure 9B), where one stable and one unstable manifold of the saddle node collide and form a homoclinic orbit; the remaining two manifolds lie outside the homoclinic orbit. At low values of Iapp, one of the unstable manifolds connects to the less depolarized stable fixed point, traversing around the more depolarized stable fixed point as well as around the stable manifold that emerges from the unstable periodic orbit around that fixed point (Figure 9B, left). Increasing Iapp beyond the bifurcation point, the unstable manifold for the saddle node converges onto itself and the stable manifold is now on the outside of the resulting limit cycle (Figure 9B, middle right). The saddle point and the stable node both remain, hence the model is bistable. A constant current pulse does not induce repetitive spiking, because the state of the system is not pushed beyond the stable manifold. Applying an additional short pulse can push the system onto a trajectory that converges onto the periodic orbit (Figure 9B, right and Figure 9C, yellow line). However, the range of Iapp-values where this occurs is extremely narrow, suggesting that it is not likely to be observed experimentally.
As expected for a system near a saddle node bifurcation, the model exhibits integrator properties, as revealed by swept sinusoidal current injection (Figure 9D). The PRC is biphasic (Figure 9E) with an r-value above 0.175, but differs from those shown in Figure 8. A perturbation shortly before, during and after the peak of an AP results in a delay of the subsequent AP. For small Iapp, the region where a perturbation causes the strongest phase advance is shortly after the peak of the AP and much earlier during the spiking cycle at a phase of approximately 0.1. Increasing Iapp to 0.3 nA, the peak of the phase advance shifts to the right, and the PRC has more resemblance to the ones shown in Figure 8.
3.6. Impact of Leak Conductance
An estimate for the input resistance (Rin) is obtained by applying Ohm's law to the steady state voltage response for small Iapp. With the chosen value of 0.036 for , depending on the Shab density, Rin is between 59.4 and 70.2 MΩ (Figure 10A, top left). is the dimensionless ratio of and . This is below the experimentally measured Rin in MN5 of 97±31 MΩ (Duch et al., 2008). However, the amplitude of Iapp for which the stable fixed point loses stability and/or a stable limit cycle emerges for small is below 0.4 nA (Figure 10A, middle left), which is the typical amplitude that induces repetitive spiking in MN5. For higher , on the other hand, the spiking threshold is too high. Rin and the inversely related amplitude of Iapp that elicits APs can be adjusted by changing . We decrease to 0.02 and investigate whether this changes the bifurcations the system exhibits. As before, was set to different values and the fixed points, periodic solutions and bifurcations that occur as Iapp is changed are determined using numerical continuation. To resolve whether a limit cycle occurs via a SNIC or homoclinic bifurcation, we compare the minimal Iapp for which a periodic solution was found with Iapp of the saddle node bifurcation (Figures 10A,B, middle). Equal amplitudes of Iapp indicate a SNIC, otherwise a saddle homoclinic bifurcation is indicated. To further distinguish between a small and big homoclinic bifurcation, the value of Vm at the fixed point and at the minimum of the periodic orbit are compared (Figures 10A,B, bottom). This analysis reveals that with increased , the same bifurcation types can be elicited by changing . The change between bifurcations may occur at different levels; i.e. with stability is lost via Hopf bifurcation when , while it is lost via a saddle node bifurcation when . However, increasing leads to spiking onset in response to smaller Iapp (Figure 10A, middle) and decreasing further ultimately results in an unstable steady state without current injection, with small . Therefore, a loss of stability via a Hopf bifurcation is the only possibility for a model that is quiescent with Iapp = 0 nA. Also, it cannot be excluded that there exists a -value for which the model exhibits a saddle small homoclinic bifurcation resulting in the appearance of a stable limit cycle, but this value was not found.
Figure 10. Small changes of adjust for Rin, but have no influence on bifurcation types the system can exhibit. Different aspects of the model behavior as function of for two different . In panels (A) and (B) the default and changed Shab channel models are used, respectively. Top panels show Rin measured with small current injections. Middle panels show the minimal amplitudes of Iapp for which a stable limit cycle emerges (red pluses) and the fixed point loses stability (blue and green crosses). Bottom panels show Vm of the fixed points near the bifurcation (blue and green crosses) and the minimum and maximum of the limit cycles (red pluses).
The chosen default value for the leak conductance is a compromise in the attempt to match sensitivity to current injection and input resistance with experimental data. The apparent missmatch might be a consequence of the fact that the model incorporates only two channel types. With the changed Shab parameters, the default value of with 0.036 results in higher Rin, but in turn, the bifurcations occur at rather small values of Iapp. Therefore, is increased to 0.05. Figure 10B demonstrates that the change in is not the cause of the saddle small homoclinic bifurcation resulting in the appearance of a stable limit cycle.
3.7. Independent Na+ Channel Inactivation
To further investigate whether the changed bifurcation structure found with the new parameter set sbz is a consequence of altering the activation of the Shab current, the inactivation of the Na+ current or whether it is necessary to alter both channels, we reintroduce an additional variable representing the gating of the Na+ channel inactivation. This allows for a better fit of the Na+ inactivation as demonstrated in Section 3.1. Increasing the gating charge results in a good match to the inactivation curve derived by Olson et al. (2008) (Figure 2B). However, the time constants are not available from these data. Using the same values as before, the model exhibits APs with a peak of -10 mV and the maximal dVm∕dt is below 20 mV/ms. This is adjusted by increasing to 20 µS and decreasing to 0.02.
With the three dimensional formulation, similarly to the two dimensional model, high Shab densities result in a Hopf bifurcation, while low densities result in a saddle node bifurcation (Figure 11). Firing patterns similar to those observed with the two dimensional model are generated. With very low levels, the branch of the minimas of the periodic solution coincides with the unstable branch of saddle points (Figure 11A, upper left panel and inset). The minimum of the periodic orbit is therefore above the stable fixed point, which coexists at the stimulus amplitude where the limit cycle emerges. In the two dimensional model, this was due to a saddle small homoclinic bifurcation (see Section 3.5). As in the two dimensional case, the membrane responds with a long delay to first spike and subsequent higher firing frequency to small square pulse current injections (Figure 12A, first panel), the model acts as low-pass filter (Figure 12B, first panel), and the PRC is biphasic displaying a phase delay shortly before and after an AP. A phase advance occurs at < 0.1 and > 0.9 of the phase; systems with other bifurcation structures display only weak responsiveness to perturbations at this time within the cycle (Figure 12C). Increasing to 0.5, stability is lost via a SNIC bifurcation (Figure 11A, lower left panel and inset). The model exhibits infinitely low frequency and a delay to first spike with the same duration as the following ISI (Figure 12A second panel), low-pass filter (Figure 12B, second panel) properties, and a type I PRC (Figure 12C). With the stable fixed point still loses stability via a saddle node bifurcation; however, a periodic orbit emerges before that point (Figure 11A), lower left panel and inset), with the minimal value of Vm smaller than for the stable fixed point. In the two dimensional model, phase plane analysis revealed that this is due to a saddle big homoclinic bifurcation. With -values of 0.7 and 0.8, a single spike is induced for sufficiently low current stimulation. With higher current injection, a shorter delay to first spike than the following ISI (Figure 12A third and fourth panel) is seen. The PRCs are similar as well (Figure 12C), but the model shows resonator properties only when stability is lost via a Hopf bifurcation with (Figure 12B).
Figure 11. Bifurcation structure for different Shab channel densities in the three dimensional model. (A) Bifurcation diagram shows that increasing changes the transition from rest into spiking from a saddle node bifurcation (blue points) to a Hopf bifurcation (green points). Fixed points are drawn in black to gray and the minimum and maximum of the periodic solutions shown in red to yellow for : 0.3, 0.5, 0.7 and 0.8. Dashed and solid lines mark stability (solid: stable; dashed unstable). (B) Two parameter bifurcation diagram indicating the location of Hopf points (blue) and saddle nodes (green) with respect to .
Figure 12. Firing behavior and PRC for different Shab channel densities in the three dimensional model. Responses to (A) steady and (B) swept sine (frequencies 0–30 Hz) current injection for (from top to bottom) : 0.3, 0.5, 0.7, 0.8; with Rin 94.5, 82.7, 75.1 and 72.1 MΩ. (C) PRCs for increasing Iapp, as indicated.
In summary, with this model, smaller current injection leads to repetitive spiking, and the changes between bifurcation types occur for smaller . Furthermore, increasing the stimulus amplitude results in spikes with smaller amplitude, and the transition out of spiking at high current levels occurs via a subcritical Hopf bifurcation. As in the two dimensional model with the changed parameters for the Shab channels, the three dimensional model displays four types of bifurcations when varying the Shab channel density by a factor of 4. This demonstrates that the additional bifurcation due to changed parameters in the two dimensional model (see Section 3.5) can be achieved by simply adjusting the Na+ channel inactivation kinetics.
In this study, we develop a two dimensional excitable membrane model with currents based on biophysical features of Drosophila Shab and DmNav29 channels. This model was used to assess which features of experimentally observed firing behavior can be produced by varying only the Shab channel density. Using phase plane and bifurcation analysis, the model was investigated and tuned to reproduce experimentally observed firing patterns. The bifurcation types the models undergo as current is injected were assessed as channel densities varied using stability analysis. The two dimensional model was then extended to a three dimensional model in order to include independent Na+ inactivation.
In previous studies, we investigated the influence of A-type channels on Drosophila MN5 firing behavior and found that, in accordance with experimental results, Shal and Shaker can influence the firing patterns and bifurcation structure of a neuron differently (Herrera-Valdez et al., 2009, 2010). It is known that regenerative currents like the persistent Na+ and A-type K+ current support aggregator properties connected to SNIC bifurcations, the type I PRCs, and class I f –I curves. On the other hand, for example the delayed rectifier K+ current evokes resonator properties associated with a Hopf bifurcation, type II PRCs, and the class II f –I curves (Hutcheon and Yarom, 2000; Ermentrout et al., 2001). However, in order to determine for which behaviors the incorporation of K+ channels with voltage dependent inactivation kinetics are necessary, it is useful to know whether those properties can be induced without those channels.
We find that a wide array of firing patterns can be generated with only two currents. Specifically, properties like the long delay to first spike and slow firing frequencies, that were thought to be unique features requiring A-type channels (Choi et al., 2004; Schaefer et al., 2010; Ping et al., 2011), can be induced in our model.
4.1. Minimal Excitable Membrane Model
We developed a minimal model consisting of one regenerative current that amplifies changes, one linear leak current, and one restorative current that counteracts changes of Vm. The regenerative current is a Na+ current; the restorative current is a delayed rectifier K+ current based on channels encoded by Shab. In Drosophila there are only four genes, Shaw, Shab, Shal, and Shaker, that encode voltage gated K+ channels. Shab channels were chosen because it has been shown that they provide delayed rectification in Drosophila and are required to induce repetitive firing, which also can be inferred by the different dynamics of the channels. Shaw channels have a very low voltage sensitivity and mediate a leak current. Shaker and Shal channels mediate transient A-type current, where their steady state activation curves are more hyperpolarized than the one associated with Shab channels, and they express fast inactivation kinetics. This indicates a role in regulating spike initiation and repetitive firing patterns rather than a major contribution to inducing repetitive firing. Additionally, Shaker mutants display broader APs and higher firing frequencies, while Shal has been reported to increase the spiking threshold and to induce a delay to first spike (Choi et al., 2004; Ping et al., 2011).
In the model, was assigned by comparing the wave forms in response to current injection of the in situ measured Vm at the soma of MN5 and the model. We did not aim to reproduce the shape of the APs exactly, since we only incorporate a small subset of channels in our models, and it is likely that the additional channels influence the shape measured in MN5. Instead the maximal dVm∕dt was compared. However, it is not clear whether the depolarization observed at the soma results from passive spread of APs, whether there is an additional spike generating zone close by or at the soma, or whether active properties at the soma boost the AP generated at the spike initiation zone (SIZ). When the SIZ is far away from the soma it is possible that the shape of APs observed at the soma is attenuated and broader than at the SIZ due to the filtering properties of passive structures. The maximal dVm∕dt measured from recordings of MN5 is small compared to values reported in other preparations, where it is about 100 and 300 mV/ms in neurons of the cat visual cortex in vivo and in vitro (Naundorf et al., 2006) or about 200 mV/ms in ganglion cells from the tiger salamander (Fohlmeister and Miller, 1997). This means that the maximal dVm∕dt at the SIZ, and thus may be higher.
4.2. Relation of Channel Density and Firing Patterns
Our results demonstrate that small changes in are sufficient to produce a variety of different firing patterns that relate to those observed experimentally. The differences can be linked to a change in the underlying mathematical structure of the system, yielding qualitatively different transitions between rest and spiking. The differences in observed firing behaviors of MN5 in situ indicate that the neuron might be close to a bifurcation with co-dimension greater than one, where such a change of the bifurcation type with respect to Iapp occurs. However, MNs are required to exhibit reliable outputs during behavior since they are the final processing station, and alternating input-output behavior cannot be adjusted by network properties. Therefore, the observed differences could be due to different modulatory states that tune the neurons' output for specific needs. For example, biogenic amines like octopamine and tyramine influence the take-off likelihood and flight maintenance significantly (Brembs et al., 2007). Further, MN5 exhibits tonic firing during flight, only single APs and no doublets or triplets are observed; during male courtship song, pulse firing is shown. Although this could be due to different input, the different firing patterns may also be supported by differences in the intrinsic excitability and thus the modulatory state.
In the model, low Shab channel concentrations induce behavior thought to be induced by A-type channels. However, at small , the range of stimulus amplitudes for which a stable periodic solution exists is rather narrow. In the three dimensional model and the model with altered channels kinetics, the AP amplitude increases drastically with increased current.
Some experimentally observed features of the firing patterns could not be reproduced well by this model. Further channels may be necessary to support some behavior combinations like integrator properties together with a higher input resistance and a higher firing threshold. Furthermore, adaptation could not be generated, since it requires a third time scale, and therefore a further variable (Guckenheimer et al., 1997).
Bistability occurs in the model for a very specific combination of constant stimulus amplitude and timing of pulse stimulation. Although bistability has not been reported for experimental studies, the model results suggest that it is unlikely this phenomenon would be observed. If observed, the neuron might be rejected and assessed as being not intact. Further, a long delay to first spike with subsequent ISI of similar duration, as typical for a system close to a SNIC bifurcation, was not reported as an observed firing pattern in Herrera-Valdez et al. (2013). However, this firing pattern was recorded at least once in eag and Shaker double knock downs, as shown in the Figure 2C of Duch et al. (2008). It must be noted that the quantification of firing patterns of wild type MN5 as reported in Duch et al. (2008) and Herrera-Valdez et al. (2013) is different, indicating that perhaps the recording conditions changed between these sets of experiments.
4.3. Variability in Ion Channel Density
It has been shown that ion channel expression, as well as the maximal conductances of ion currents in neurons, can be extremely variable. Shab mRNA in MNs of the crustacean stomatogastric nervous system vary two- to four-fold (Schulz et al., 2006, 2007) and four-fold in gastric mill neurons. In Drosophila MN5, voltage clamp recordings reveal that the peak of the total K+ current can vary by a factor of about four (Ryglewski and Duch, 2009). In our models a variation of up to four-fold is sufficient to produce the different firing patterns and bifurcation structures. However, it has also been shown that the expression patterns of certain ion channel combinations are coordinated, where the specific combinations depend on the cell type. In gastric mill neurons Na+ and Shab channels are correlated, but in all other investigated cells, the mRNA abundance of the genes encoding these two channels is independent (Schulz et al., 2007). In a recent review article Marder (2011) pointed out that these findings indicate the need “to measure as many system components as possible within an individual” and that the process of fitting and combining independently obtained and averaged experimental data is unlikely to reproduce realistic model behavior. Especially in the small Drosophila, it is hard to measure a variety of system parameters simultaneously. Similarly, pharmacological manipulations may suffer from low specificity of drugs (Greenwood and Leblanc, 2007). These issues can potentially lead to incorrect conclusions about the role of a certain ion channel, however, the use of computational models can help resolve these issues, by allowing for controlled and independent variation of specific parameters to understand the impact of biological variations and allow for an examination of the impact of parameters that are not experimentally separable.
4.4. Influence of Phasic Input on Output Behavior
During flight, the firing frequency of the MNs is modulated simultaneously (Levine and Wyman, 1973) and studies support the idea that the firing frequency is an inherent property of the neurons (Harcombe and Wyman, 1977), which changes when the required power output is changed (Gordon and Dickinson, 2006). Together with the observation that two MNs never fire in a time window of about 5 ms, it was concluded that all MNs receive common excitatory input and that the MNs interact reciprocally (Harcombe and Wyman, 1977). In other words, the common excitatory input induces periodic spiking of the units, and the reciprocal interactions cause a perturbation of the period depending on the phase when the input is received. An analysis of how a neuron responds to phasic input while firing repetitively is measured with PRCs. The claim of Harcombe and Wyman (1977) and Koenig and Ikeda (1983) that there is a strong interaction if one MN fires simultaneously or shortly before another MN, requires a PRC with large phase shifts at the very beginning and end of the period, as observed in models close to a saddle small homoclinic circle. Those models also display a long delay to first spike with shorter ISIs to constant current injection, a firing pattern that is often observed in neurons at larval stage, but rather seldom in adult MN5s. However, in our model, this pattern corresponds to current injections very close to the spiking threshold. Therefore, it could also mean that this pattern manifests at stimulus amplitudes that fall between the probed current steps.
In the two dimensional model with the default channel kinetics, a delay to first spike could not be reproduced. We have shown with the electro-diffusion model, that this does not require the introduction of a further variable (Herrera-Valdez et al., 2013) rather a change of zb in combination with a low -value elicits this behavior. In the model investigated here, a simple change of zb is not sufficient to produce this profile. However, it is produced by models with moderate changes in more than one parameter and more importantly, after incorporating independent Na+ inactivation dynamics.
4.5. Comparison to Other Excitable Membrane Models
Various minimal conductance based models have been employed in previous studies in order to analyze the electrical behavior of neurons. For example Prescott et al. (2008) employed a two dimensional model based on the Morris–Lecar model (Morris and Lecar, 1981; Rinzel and Ermentrout, 1998) and showed that all three types of excitability could be reproduced. However, here we show that a neuron model based on more realistic ion channel mechanisms also can exhibit all three types of firing behavior. Zeberg et al. (2010) used models in which diverse K+ currents were combined as one recovery variable and found that the variation of channel densities can switch the model from resonator to integrator. This study demonstrates that incorporating only one delayed rectifier K+ channel in the model and varying its density is sufficient to induce this switch.
In an earlier modeling effort, we contributed to the development of an electro-diffusion based model for MN5 membrane potential dynamics (Herrera-Valdez et al., 2013). However, this previous model falls short in reproducing the Na+ channel kinetics, and it is not clear whether the theoretical improvement of the electro-diffusion based model justifies the accompanied increase in complexity. In general, an advantage of the electo-diffusion approach is that the maximal currents can be fitted directly to electrophysiological data. However, electrophysiological data are often preprocessed under the assumption of a linear current–voltage relationship, for example when leak currents are subtracted from the traces. Unlike the conductance-based model that we present here, the previous electro-diffusion based model also fails to generate spikes on elevated potentials as observed experimentally. However, this in situ observation can also emerge when the recording site is far away from the SIZ, which cannot be addressed with a minimal model approach employed by Herrera-Valdez (2012) and here.
In another recent study, a model of Drosophila third-instar larval motoneuron was developed in order to investigate the impact of alternative splicing of Na+ channels (Lin et al., 2012). In that model the Na+ channels have activation kinetics as in the classical Hodgkin–Huxley model formalism, and the time constant is based on electrophysiological recordings in embryonic Drosophila motoneuron (O'Dowd and Aldrich, 1988). The measured Na+ current is possibly a combination of currents mediated by different splice variants of the DmNav29 gene that are specific to embryonic neurons. Here, we use instantaneous activation of Na+ channels for two reasons: first due to the advantage of reducing the model to only two dimensions and second because there are no data for the time course of specific splice variants available.
In conclusion, we find that a wide array of firing patterns like those seen in experimental studies can be generated with only two currents, and we show that a neuron model based on realistic ion channel mechanisms can exhibit all three types of firing behavior. This study demonstrates that incorporating only one delayed rectifier K+ channel in the model and increasing its density is sufficient to induce a switch from type I excitability to type II excitability, with a corresponding switch from integrator to resonator properties. Additionally, properties like the long delay to first spike and slow firing frequencies, that were thought to be unique features requiring A-type channels (Choi et al., 2004; Schaefer et al., 2010; Ping et al., 2011), can be induced in our model. As described above, additional channel types are known to be present and play a role in Drosophila motoneuron dynamics, leading to many possible, more complex alternatives for obtaining similar results. It is plausible that having more channels that in different combinations produce a similar behavior can add to the robustness of the system and provide a way to adjust for natural variability. Certainly, the minimal model described here can guide our understanding of the interactions among channel types in this important model system.
SB and SC were supported in part by the National Science Foundation (NSF IIS 0613404). Support for SB was also provided by the Interdisciplinary Graduate Program in Neuroscience at Arizona State University.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We thank Stefanie Ryglewski and Carsten Duch for providing the electrophysiological recordings from MN5 shown in this paper.
Brembs, B., Christiansen, F., Pflüger, H. J., and Duch, C. (2007). Flight initiation and maintenance deficits in flies with genetically altered biogenic amine levels. J. Neurosci. 27, 11122–11131. doi: 10.1523/JNEUROSCI.2704-07.2007
Choi, J. C., Park, D., and Griffith, L. C. (2004). Electrophysiological and morphological characterization of identified motor neurons in the Drosophila third instar larva central nervous system. J. Neurophysiol. 91, 2353. doi: 10.1152/jn.01115.2003
Clewley, R., Sherwood, W., LaMar, M., and Guckenheimer, J. (2007). Pydstool, A Software Environment for Dynamical Systems Modeling. Available online at: http://pydstool.sourceforge.net
Coetzee, W. A., Amarillo, Y., Chiu, J., Chow, A., Lau, D., McCormack, T., et al. (1999). Molecular diversity of K+ channels. Ann. N. Y. Acad. Sci. 868, 233–255. doi: 10.1111/j.1749-6632.1999.tb11293.x
Destexhe, A., and Huguenard, J. R. (2000). “Which formalism to use for modeling voltag-dependent conductances?,” in Computational Neuroscience: Realistic Modeling for Experimentalists, ed E. DeSchutter (Boca Raton, FL: CRC Press), 129–157.
Doedel, E., Paffenroth, R. C., ChampneyS, A. R., Fairgrieve, T., Kuznetsov, Y. A., Oldeman, B., et al. (2002). AUTO2000: Continuation and bifurcation software for ordinary differential equations. Technical Report, California Institute of Technology.
Duch, C., Vonhoff, F., and Ryglewski, S. (2008). Dendrite elongation and dendritic branching are affected separately by different forms of intrinsic motoneuron excitability. J. Neurophysiol. 100, 2525–2536. doi: 10.1152/jn.90758.2008
Ermentrout, B., Pascal, M., and Gutkin, B. (2001). The effects of spike frequency adaptation and negative feedback on the synchronization of neural oscillators. Neural Comput. 13, 1285–1310. doi: 10.1162/08997660152002861
Gasque, G., Labarca, P., Reynaud, E., and Darszon, A. (2005). Shal and Shaker differential contribution to the K+ currents in the Drosophila mushroom body neurons. J. Neurosci. 25, 2348–2358. doi: 10.1523/JNEUROSCI.4384-04.2005
Gutkin, B., Ermentrout, G. B., and Rudolph, M. (2003). Spike generating dynamics and the conditions for spike-time precision in cortical neurons. J. Comput. Neurosci. 15, 91–103. doi: 10.1023/A:1024426903582
Herrera-Valdez, M. A., Berger, S., Duch, C., and Crook, S. (2009). Predicting changes in neuronal excitability type in response to genetic manipulations of K+-channels. BMC Neurosci. 10(Suppl. 1):P301. doi: 10.1186/1471-2202-10-S1-P301
Herrera-Valdez, M. A., Berger, S., Duch, C., and Crook, S. (2010). Differential contribution of voltage-dependent potassium currents to neuronal excitability. BMC Neurosci. 11:P159. doi: 10.1186/1471-2202-11-S1-P159
Herrera-Valdez, M. A., McKiernan, E. C., Berger, S. D., Ryglewski, S., Duch, C., and Crook, S. (2013). Relating ion channel expression, bifurcation structure, and diverse firing patterns in a model of an identified motor neuron. J. Comput. Neurosci. 34, 211–229. doi: 10.1007/s10827-012-0416-6
Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. doi: 10.1113/jphysiol.1952.sp004764
Koenig, J. H., and Ikeda, K. (1983). Reciprocal excitation between identified flight motor neurons in Drosophila and its effect on pattern generation. J. Comp. Physiol. A 150, 305–317. doi: 10.1007/BF00605020
Lin, W.-H., Günay, C., Marley, R., Prinz, A. A., and Baines, R. A. (2012). Activity-dependent alternative splicing increases persistent sodium current and promotes seizure. J. Neurosci. 32, 7267–7277. doi: 10.1523/JNEUROSCI.6042-11.2012
Lin, W.-H., Wright, D. E., Muraro, N. I., and Baines, R. A. (2009). Alternative splicing in the voltage-gated sodium channel DmNav regulates activation, inactivation, and persistent current. J. Neurophysiol. 102, 1994. doi: 10.1152/jn.00613.2009
Mee, C. J., Pym, E. C., Moffat, K. G., and Baines, R. A. (2004). Regulation of neuronal excitability through pumilio-dependent control of a sodium channel gene. J. Neurosci. 24, 8695. doi: 10.1523/JNEUROSCI.2282-04.2004
Miyazaki, M., Ohyama, K., Dunlap, D. Y., and Matsumura, F. (1996). Cloning and sequencing of the para-type sodium channel gene from susceptible and kdr-resistant German cockroaches (Blattella germanica) and house fly (Musca domestica). Mol. Gen. Genet. 252, 61–68. doi: 10.1007/s004389670007
Olson, R. O., Liu, Z., Nomura, Y., Song, W., and Dong, K. (2008). Molecular and functional characterization of voltage-gated sodium channel variants from Drosophila melanogaster. Insect Biochem. Mol. Biol. 38, 604–610. doi: 10.1016/j.ibmb.2008.01.003
Ping, Y., Waro, G., Licursi, A., Smith, S., Vo-Ba, D.-A., and Tsunoda, S. (2011). Shal/Kv4 channels are required for maintaining excitability during repetitive firing and normal locomotion in Drosophila. PLoS ONE 6:e16043. doi: 10.1371/journal.pone.0016043
Prescott, S. A., De Koninck, Y., and Sejnowski, T. J. (2008). Biophysical basis for three distinct dynamical mechanisms of action potential initiation. PLoS Comput. Biol. 4:e1000198. doi: 10.1371/journal.pcbi.1000198
Rinzel, J., and Ermentrout, B. (1998). “Analysis of neural excitability and oscillations,” in Methods in Neuronal Modeling: From Ions to Networks, eds C. Koch and I. Segev (Cambridge, MA: The MIT Press), 251–292.
Ryglewski, S., and Duch, C. (2009). Shaker and Shal mediate transient calcium-independent potassium current in a Drosophila flight motoneuron. J. Neurophysiol. 102, 3673–3688. doi: 10.1152/jn.00693.2009
Schaefer, J. E., Worrell, J. W., and Levine, R. B. (2010). Role of intrinsic properties in Drosophila motoneuron recruitment during fictive crawling. J. Neurophysiol. 104, 1257. doi: 10.1152/jn.00298.2010
Schulz, D. J., Goaillard, J. M., and Marder, E. (2006). Variable channel expression in identified single and electrically coupled neurons in different animals. Nat. Neurosci. 9, 356–362. doi: 10.1038/nn1639
Schulz, D. J., Goaillard, J. M., and Marder, E. E. (2007). Quantitative expression profiling of identified neurons reveals cell-specific constraints on highly variable levels of gene expression. Proc. Natl. Acad. Sci. U.S.A. 104, 13187. doi: 10.1073/pnas.0705827104
Tateno, T., Harsch, A., and Robinson, H. (2004). Threshold firing frequency-current relationships of neurons in rat somatosensory cortex: type 1 and type 2 dynamics. J. Neurophysiol. 92, 2283. doi: 10.1152/jn.00109.2004
Wei, A., Covarrubias, M., Butler, A., Baker, K., Pak, M., and Salkoff, L. (1990). K+ current diversity is produced by an extended gene family conserved in Drosophila and mouse. Science 248, 599. doi: 10.1126/science.2333511
Willms, A. R., Baro, D. J., Harris-Warrick, R. M., and Guckenheimer, J. (1999). An improved parameter estimation method for hodgkin-huxley models. J. Comp. Neurol. 6, 145–168. doi: 10.1023/A:1008880518515
Keywords: membrane excitability, ion channel kinetics, Drosophila model, neuronal dynamics, bifurcation studies
Citation: Berger SD and Crook SM (2015) Modeling the Influence of Ion Channels on Neuron Dynamics in Drosophila. Front. Comput. Neurosci. 9:139. doi: 10.3389/fncom.2015.00139
Received: 20 May 2015; Accepted: 28 October 2015;
Published: 18 November 2015.
Edited by:Boris Gutkin, École Normale Supérieure, France
Reviewed by:Douglas Zhou, Shanghai Jiao Tong University, China
Baktash Babadi, Harvard University, USA
Copyright © 2015 Berger and Crook. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Sharon M. Crook, firstname.lastname@example.org