Abstract
Modeling is a useful tool for investigating various biophysical characteristics of neurons. Recent simulation studies of propagating action potentials (spike conduction) along axons include the investigation of neuronal activity evoked by electrical stimulation from implantable prosthetic devices. In contrast to point-neuron simulations, where a large variety of models are readily available, Hodgkin–Huxley-type conductance-based models have been almost the only option for simulating axonal spike conduction, as simpler models cannot faithfully replicate the waveforms of propagating spikes. Since the amount of available physiological data, especially in humans, is usually limited, calibration, and justification of the large number of parameters of a complex model is generally difficult. In addition, not all simulation studies of axons require detailed descriptions of nonlinear ionic dynamics. In this study, we construct a simple model of spike generation and conduction based on the exponential integrate-and-fire model, which can simulate the rapid growth of the membrane potential at spike initiation. In terms of the number of parameters and equations, this model is much more compact than conventional models, but can still reliably simulate spike conduction along myelinated and unmyelinated axons that are stimulated intracellularly or extracellularly. Our simulations of auditory nerve fibers with this new model suggest that, because of the difference in intrinsic membrane properties, the axonal spike conduction of high-frequency nerve fibers is faster than that of low-frequency fibers. The simple model developed in this study can serve as a computationally efficient alternative to more complex models for future studies, including simulations of neuroprosthetic devices.
Significance Statement
Conduction of electrical impulses (action potentials) along the axon is essential for information transfer between neurons. Simulation studies of propagating action potentials, which earlier focused on the biophysical mechanisms of conduction, have progressed to investigations of pathologic malfunctions of nerves and electrical stimulations via prostheses. In contrast to dimensionless, single-neuron modeling, for which a number of different approaches are available, simulation of nerve conduction generally requires a complex model of ionic conductances to reproduce propagating action potentials. In this study, we present a simplified phenomenological model of axonal conduction with increased computational efficiency and a reduced number of parameters. This simple model can be used as an alternative to conventional models, especially for applications including prosthetic simulations of nerve conduction.
Introduction
Since Louis Lapicque (Lapicque, 1907) first introduced its underlying idea of parallel capacitor and leak resistor in combination with threshold crossing, the integrate-and-fire (IF) model has served as a useful tool to simulate spiking activity of neuronal membranes. Even after the detailed ionic dynamics underlying spike generation were discovered and modeled with a more complex conductance-based description (Hodgkin and Huxley, 1952), the IF model and its variations are still frequently used as a convenient alternative to Hodgkin–Huxley (HH)-type models, especially when computational efficiency and mathematical transparency are required. Applications of the IF model include large-scale simulations of a neuronal network, and rigorous analysis of neuronal spiking responses driven by random synaptic inputs (for representative examples, see Koch, 1999; Gerstner et al., 2014).
Within the family of IF-type models, various nonlinear versions have been created (Gerstner et al., 2014). For instance, the exponential IF (EIF) model was introduced to describe the exponential growth of the membrane potential at spike initiation (Fourcaud-Trocmé et al., 2003). Its subthreshold response properties to stimulus current injections remain largely unchanged from the Wang–Buzsáki (WB) model, which itself is a modification of the HH model (Wang and Buzsáki, 1996). In the category of single-compartment models, a modified version of the EIF model was shown to replicate the rapid initiation of action potentials even better than more detailed HH-type models (Brette, 2015). Although the EIF model was originally created to simulate the spiking activity of cortical neurons, its variations are now used to simulate a wider range of cells [auditory nerve (AN) fibers: Rutherford et al., 2012; Joshi et al., 2017; visuomotor system: Morén et al., 2013; cerebellar Purkinje cell: Ostojic et al., 2015; optic nerves: Arancibia-Cárcamo et al., 2017].
In the domain of single-compartment models, several levels of abstraction are possible: from biophysical conductance-based descriptions equipped with a variety of ion channels, via IF-type models with intermediate complexity, to more phenomenological “black-box” approaches that focus solely on the input-output functions (Herz et al., 2006). This gradient of biological plausibility and computational efficiency enables a user to select an appropriate single-compartment model depending on the specific purpose of modeling (Ashida et al., 2017). In contrast, for multicompartment neuronal modeling, in which multiple nonlinear excitable units are connected with each other, HH-type models are normally the only options, since the abrupt reset of the membrane potential in an IF-type model is generally incompatible with spatially propagating electrical activity over the modeled membrane. Earlier studies using multicompartment spiking neuron models simulated the conduction of action potentials along the axon. For myelinated axons, for example, each node of Ranvier was modeled as a HH-type excitable compartment that was interconnected with an axial resistance (FitzHugh, 1962; Brill et al., 1977; Moore et al., 1978). For unmyelinated axons, the HH model was combined with the cable equation to account for the spatial extension of the axon (Cooley and Dodge, 1966). Compartmental models of spike conduction were later applied to simulate, for example, pathologic changes of axons (Coggan et al., 2010; Brown and Hamann, 2014) and the interaction between nerves and prosthetic devices (for a review, see Rattay et al., 2003).
Despite the general success of HH-type models in reproducing axonal spike conduction, not all simulation studies actually require the detailed descriptions of ion channel dynamics. Moreover, neurophysiological data from humans, in particular for single-cell properties, are usually sparse, making it difficult to calibrate or justify the parameters of a model used for prosthetic nerve simulations. Recent prosthetic modeling aims to simulate tens of nodes in thousands of nerves distributed three-dimensionally (Nogueira et al., 2016), which requires the efficient model representation of excitable units. In this study, we propose a simple model of action-potential propagation along the axon based on the EIF model of spike generation. The model has much fewer parameters than the HH model but still faithfully reproduces axonal spike conduction. As an example application, we fit the model to known physiological data from ANs. Our simulated conduction velocities match the experimentally measured range in AN fibers, confirming the applicability of the model. We expect that the model introduced in this study will serve as a simpler replacement for the HH model especially when computational performance and structural simplicity are preferred over the biophysical details of spike generation.
Materials and Methods
Overview of the three models
In this article, we compare three types of spiking membrane models: (1) the WB model, a variation of the HH model, having nonlinear sodium and potassium conductances (Wang and Buzsáki, 1996); (2) the original version of the EIF model (Fourcaud-Trocmé et al., 2003), which we refer to as the standard EIF (sEIF) model; and (3) a modified version of the EIF model, called the bounded EIF (bEIF) model, with a “ceiling” for the spike-generating, depolarizing current. The sEIF model was originally created (and fitted) to study the spike generation of the WB model (Fourcaud-Trocmé et al., 2003), and we here introduce the bEIF model as a modification of the sEIF model to account for axonal spike conduction.
All these models share an equation of the form:
(1)
where Cm is the membrane capacitance, GL is the time- and voltage-independent (linear) leak conductance, EL is the leak reversal potential, Iinj is the intracellularly injected current, and ψ(V) is a nonlinear function of the membrane potential V responsible for spike generation. In the WB model, ψ(V) is a sum of sodium INa and potassium IK currents:
(2a)
The activation/inactivation kinetics of these currents are described by three additional differential equations (Table 1). In the sEIF model, ψ(V) is equal to the depolarizing current Idep:
(2b)which is an exponential function of the membrane potential V describing the exponential growth of the membrane potential at spike initiation (Table 2). The after-spike repolarization in the sEIF model is simply a reset of V to the resting membrane potential. In the bEIF model, ψ(V) is a sum of depolarizing Idep and repolarizing Irep currents:
(2c)
Equations and parameters for the single-compartment WB model
Equations and parameters for the single-compartment sEIF model
The depolarizing current Idep represents the exponentially growing, inward current for spike initiation, while the repolarizing current Irep corresponds to the outward current responsible for after-spike repolarization (Table 3). More detailed descriptions of each model are provided below. As in the previous work, in which the sEIF model was first introduced (Fourcaud-Trocmé et al., 2003), we use the WB model as a reference and compare its responses with those of the EIF models.
Equations and parameters for the single-compartment bEIF model
Wang-Buzsáki model
The WB model is a set of nonlinear differential equations that describe the dynamics of the membrane potential V(t), the activation m(t) and inactivation h(t) of sodium channels, and the activation n(t) of potassium channels (Table 1). While in the original work of Wang and Buzsáki (1996), the sodium activation was assumed to be instantaneous, here we adopted a voltage-dependent time constant for sodium activation as in the original HH model (Hodgkin and Huxley, 1952). The resulting differences between instantaneous and time-delayed sodium activations are generally minor and limited to high-frequency sinusoidal input currents (Fourcaud-Trocmé et al., 2003). The membrane parameters we used (Table 1) were taken from Wang and Buzsáki (1996).
The single-compartment WB model has seven membrane parameters: one membrane capacitance (Cm), three conductances (GL, GK, GNa), and three reversal potentials (EL, EK, ENa). In addition, each rate function for channel activation and inactivation (αm, βm, αh, βh, αn, βn) requires three parameters (amplitude, reference voltage, and slope factor). In total, the WB model needs 25 parameters to be calibrated to fit physiological data.
Standard exponential integrate-and-fire model
The sEIF model, which is a nonlinear modification of the leaky (linear) IF model, phenomenologically describes the exponentially increasing sodium inward current at spike initiation (Fourcaud-Trocmé et al., 2003). Its spike-generating depolarizing current,
(3)is characterized by the (soft) threshold VT and the slope factor KT, which determine the excitability of the model neuron. Once the membrane potential V crosses the spike-detection threshold Vspike, it is reset to and held at the resting potential Vreset for the refractory period τref.
The single-compartment sEIF model has eight parameters (Table 2). Since the spiking voltage of the sEIF model quickly diverges to infinity in a finite amount of time, the spike-detection threshold Vspike does not play a major role in determining the response property of the model (Touboul, 2009), thus reducing the effective number of unconstrained parameters to seven. As in the original sEIF study (Fourcaud-Trocmé et al., 2003), the parameters of the sEIF model in our study were selected so that the initial part of its spike waveform, including the subthreshold response (Fig. 1A), its spiking threshold, and its frequency-current (f-I) relationship resembled those of the WB model (Fig. 1B).
Response properties of single-compartment models. A, Spike responses of the sEIF (green), bEIF (blue), and WB (gray) models driven by a step current input. B, f-I curves of the models. Step currents of varied amplitudes were injected and the numbers of spikes in 1000 ms were calculated. C, Voltage dependence of the exponential growth factors of the sEIF and bEIF models. D, Depolarizing and repolarizing spike currents of the bEIF model. The horizontal corresponds to the expanded time in A. E, Spike shapes of the bEIF and WB models. F, Membrane currents of the bEIF and WB models. In E, F, traces are aligned such that time 0 corresponds to the peak timings of the action potentials shown in A.
Bounded exponential integrate-and-fire model
As we will see in the Results, the behavior of the sEIF membrane potential diverging to infinity (Fig. 1A) is incompatible with spike propagation along the axon. Hence, we replaced the exponential growth term exp((V-VT)/KT) in the sEIF model (Eq. 3) with AT/(1+AT exp(-(V-VT)/KT)) to set a ceiling AT for the inward current (Fig. 1C). Thus, in this modified model, named the bEIF here, the spike-generating current is written as:
(4)and is bounded as: max(Idep) = GLKTAT. This modification resulted in a slower (and probably more realistic) voltage increase near the peak of an action potential while keeping the sub- and near-threshold responses almost identical to the sEIF model (Fig. 1A).
In contrast to the instantaneous potential reset of the sEIF model, the bEIF model has an additional repolarizing current Irep to mimic the downward trajectory of the membrane potential after each spike (Fig 1A; for the equations, see Table 3). This current is initiated when the potential V reaches the preset starting voltage Vrep, rapidly overwhelming the depolarization current (Fig. 1D) to bring the membrane potential back to the resting level (Fig. 1A). We used an alpha function for the repolarizing conductance, since it allows fast and exact calculation at each time step (Rotter and Diesmann, 1999). Simulated spike shapes were similar between the single-compartment WB and bEIF models, except for the slightly narrower spike width and the lack of after-spike hyperpolarization in the bEIF model (Fig. 1E). The total membrane currents during a spike are also comparable between these two models, both in amplitude and time course (Fig. 1F).
The bEIF model has nine parameters (Table 3): three for subthreshold responses (Cm, GL, and EL), another three for spiking (VT, KT, and AT), and the remaining three for repolarization (Vrep, τrep, and Arep). In this study, the ceiling factor AT and three repolarization parameters of the bEIF model were adjusted to mimic the spike shape (Fig. 1A) and f-I curve (Fig. 1B) of the WB model, while the other five parameters were unchanged from the sEIF model. It should be noted that, unlike the sEIF model, the bEIF model does not have an explicit refractory period as a model parameter, because the repolarizing current Irep, which rapidly overcomes the spike current Idep, effectively suppresses spike generation for a certain time period. The length of this “dead time” is determined by the shape (amplitude and time scale) of the repolarizing current.
Simulating myelinated axons
Myelinated axons were modeled as a series of excitable units interconnected with an axial resistance (Table 4). The excitability of each nodal compartment is described either by the WB model or the bEIF model. For simplicity, we considered the ideal situation, in which the myelinated internodes are perfectly insulated (i.e., with negligible capacitance and transmembrane conductance; McNeal, 1976; Keener and Sneyd, 2009), although simulations suggested imperfect insulation might affect both excitability and conduction (McIntyre et al., 2002; Young et al., 2013). Furthermore, we also simply assumed that all ion channels (of the WB model) are located at the nodes of Ranvier, despite the accumulating evidence of non-uniform distribution of ion channels at and around the node (Hossain et al., 2005; Yi et al., 2010; Kim and Rutherford, 2016; for reviews, see Salzer et al., 2008; Debanne et al., 2011; Freeman et al., 2016). The default parameter values of our myelinated axon model are shown in Table 4. For each excitable node, we used the same WB (Table 1) or bEIF (Table 3) description as for the single-compartment model.
Equations and parameters for myelinated axon models
We simulated 141 nodes along a one-dimensional (non-branching) axon. Stimulus currents (amplitude Iinj = 100 pA, duration T = 1 ms) were injected intracellularly into the node #20 to evoke action potentials. In simulations where the axonal diameter D (μm) was changed, the current amplitude was linearly adjusted with the diameter as Iinj = (D/2) × 100 pA to securely initiate spikes. To estimate the conduction velocity, we measured the travel time between nodes #40 and #90 by calculating the difference of the times at which the membrane potential reached its peak at these nodes. Then we divided the distance between the two nodes by the travel time to obtain the conduction velocity.
Simulating unmyelinated axons
Unmyelinated axons were simulated as a series of excitable elements combined with the cable model (Koch, 1999). While a partial differential equation is used for the mathematical formulation, the simulated axon is actually divided into compartments when the propagation of an action potential is numerically calculated (Cooley and Dodge, 1966). Table 5 summarizes the equations and default parameters of the unmyelinated axon model. Similarly to the myelinated axon model, we used either the WB model (Table 1) or bEIF model (Table 3) for each compartment to simulate the spikes traveling along unmyelinated axons.
Equations and parameters for unmyelinated axon models
We simulated 301 compartments along a one-dimensional (non-branching) axon. The length of each compartment was set to 20 μm, which is sufficiently small compared to the length constant λ = 1.1 mm of this model axon (calculated as λ2 = D/(4GLRax) = 1.25 mm2). Stimulus currents (amplitude Iinj = 10 nA, duration T = 1 ms) were injected intracellularly into the compartment #50 to evoke action potentials. In simulations where the axonal diameter D (μm) was changed, the current amplitude was linearly adjusted with the diameter as Iinj = (D/10) × 10 nA, to securely initiate spikes. To estimate the conduction velocity, we measured the travel time between nodes #100 and #200 by calculating the difference of the times at which the membrane potential reached its peak at these nodes. Then we divided the distance between the two nodes by the travel time to obtain the conduction velocity.
Simulating extracellular stimulation
Extracellular stimulation of myelinated nerves can be simulated similarly to intracellular stimulation by introducing an additional variable of extracellular potential Uex (Table 6), which is inversely proportional to the distance r from the current source (i.e., dimensionless, point electrode) as: Uex = ρexIex/(4πr), with ρex being the extracellular resistivity and Iex the amount of injected current (Abzug et al., 1974; McNeal, 1976; Rattay, 1986). The intracellular axonal current Iaxon is determined by the gradient of intracellular potential Uin, which is the sum of the membrane potential V and the extracellular potential Uex (Table 6). The equations for the active and passive membrane currents that depend on the membrane potential are unchanged from the case of intracellular stimulation (Table 4).
Equations and parameters for extracellular stimulation of myelinated axon models
We simulated 141 nodes along a one-dimensional, straight axon. To evoke action potentials, stimulus currents (amplitude Iex = -1 mA, duration T = 0.1 ms) were injected into the extracellular space located 1 mm away from the node #20. We tested both the WB and the bEIF model.
Simulating AN axons
As an application of the bEIF model, we simulated spike conduction along the central axon of the mammalian AN. The model equations are the same as those for the myelinated bEIF axon (Table 4). Here we consider AN fibers that are tuned to either low frequency (located in the apex of the cochlea) or high frequency (located in the base of the cochlea). Table 7 lists the parameters for the low- and high-frequency AN models, and Table 8 summarizes relevant anatomic and physiological data used for calibrating the models. Since no physiological data were available for AN axons, we used values measured in cell bodies of spiral ganglion (AN) neurons.
Parameters for low- and high-frequency AN models
Anatomic and physiologic data used for calibrating the AN models
From the reported capacitance of 10 pF (Table 8) and a standard capacitance density of 1.0 μF/cm2, the effective surface area of the cell body is estimated as 1000 μm2. Using the leak conductance densities of 0.2 mS/cm2 (low frequency AN model) and 0.4 mS/cm2 (high-frequency AN model), the membrane resistance of a 1000 μm2 membrane patch is calculated as 500 MΩ (low frequency) and 250 MΩ (high frequency), matching the measured physiologic values of mammalian spiral ganglion neurons (Table 8). The adopted leak reversal potential of -65.3 mV was unchanged from the single-compartment bEIF model (Table 3), as it matched the measured resting potential (Table 8). The (soft) threshold VT for the spike-generating current of the EIF model is not directly related to measured spike thresholds, because of the fundamental differences in their definitions. We selected the value of VT to roughly mimic the recorded spike waveforms of ANs (Adamson et al., 2002).
Reported internodal lengths of cat AN axons (Liberman and Oliver, 1984) were larger in the high-frequency region than in the low-frequency region, while the diameters of axons were similar between these regions (summarized in Table 8). As there were no systematic measurements available for the length of the node of Ranvier, we simply assumed that both low- and high-frequency AN fibers share the same nodal length of 2 μm. Previous simulations showed that the effect of nodal length of conduction velocity is relatively minor (Arancibia-Cárcamo et al., 2017). We simulated 40 nodes along a one-dimensional (non-branching) axon. The total length of the modeled axon roughly corresponds to the length of the cat AN fiber innervating the cochlear nuclei (Ryugo and Rouiller, 1988). Stimulus currents (amplitude Iinj = 60 pA, duration T = 1 ms) were injected intracellularly to the node #1 to evoke action potentials. The conduction velocity of a propagating spike was calculated as the distance between nodes #10 and #30 divided by the travel time between these nodes.
Simulation environment
Numerical integration of the model equations was performed with the explicit (forward) Euler method, in combination with the Crank-Nicolson method for axonal current propagation (Moore et al., 1978). The time step we used was fixed to 4 μs, unless otherwise stated. To obtain an f-I curve for the single-compartment models, we injected a step current of varied amplitudes with a duration of 1000 ms to calculate the output spike rate of the model. Code was implemented with MATLAB R2015b (MathWorks).
To evaluate the computational costs of calculating axonal spike conductions with different models, we computed the integration time of voltage traces of an axon with 141 nodes, using identical configurations of the modeled axon to those described above in Simulating myelinated axons. Each trial was 400 ms long, and we repeated the computation 50 times to obtain an average integration time. In addition to MATLAB, numerical algorithms were also implemented in D (Alexandrescu, 2010), which is compiled into native machine code and is expected to run faster. Simulations were conducted on a desktop computer (Dell 1398 Precision T1700) with a 64-bit Windows 7 Professional Operating System, Intel Xeon CPU E3-1270 1399 v3 (4 core, 3.5 GHz) and 16 GB memory.
Code accessibility
MATLAB implementation of the models is publicly available online at https://github.com/pinkbox-models.
Results
Responses the bEIF model
The main goal of this study was to construct a simple model of spike conduction along the axon. To this end, we first modified the sEIF model by limiting the exponentially growing inward current (Fig. 1C) and introducing a repolarizing current after spikes (Fig. 1D; see Materials and Methods). With these modifications, the resulting single-compartment bEIF model showed spike waveforms that were more similar to those of the HH-type WB model than the sEIF model (Fig. 1A,E), while keeping its f-I relationship largely comparable to those of both the WB and sEIF models (Fig. 1B).
Spike conduction in myelinated axons
We simulated spike conduction along a myelinated axon by connecting excitable compartments with an axial resistance (Fig. 2A; Materials and Methods). The voltage dynamics of each compartment was simulated by either the WB or bEIF model. As in earlier studies with the HH model (FitzHugh, 1962; Brill et al., 1977; Rattay et al., 2003), stable propagation of an action potential was observed for the WB model (Fig. 2B). With the bEIF model, simulated spike propagation was also stable (Fig. 2C) and the estimated conduction velocity was comparable to that of the WB model. The relative insensitivity of the conduction velocity to the detailed spike-generating mechanisms was reported in an earlier study that compared the conductance-based HH model with the permeability-based Frankenhaeuser–Huxley model (Moore et al., 1978).
Response properties of myelinated axon models. A, Schematic drawing (top) and modeled electric circuit (bottom) of a myelinated axon with a diameter D, nodal length Ln and internodal length Li. Each nodal compartment has a capacitance cm and (voltage-dependent) resistance rm and is interconnected with neighboring nodes with an axial resistance ra. See Table 4 for the default parameter values. In panels B, C, G, H, voltage responses at every 10th node (i.e., each 2 mm apart) are shown. Currents were injected intracellularly into the node #20 (gray traces). B, Spike conduction along the modeled myelinated axon simulated with the WB model. C, Spike conduction along the modeled myelinated axon simulated with the bEIF model. D, Shapes of conducted spikes of the bEIF and WB models. The peaks of both traces are aligned at time 0. E, Dependence of simulated conduction velocity u (m/s) on the axon diameter D (μm). The dotted curve shows a square root fit by 
. F, Dependence of conduction velocity u (m/s) on the internodal length Li (μm). The dotted curve shows a square root fit by 
. G, Spike conduction along the modeled myelinated axon simulated with the sEIF model. Inset, Expanded traces around the spike generation. In panel G, we used a time step of 0.2 μs to faithfully simulate the rapidly increasing membrane potentials. H, Spike conduction along the modeled myelinated axon simulated with the bEIF model, with an instantaneous potential reset instead of the repolarizing current. I, Dependence of conduction velocity u on the ceiling value AT of the bEIF model (blue). Conduction velocity of the WB model (5.7 m/s) is also shown as a reference (thicker gray line). When necessary (typically for large and small values of AT), the starting voltage for repolarization current Vrep was readjusted (down to -20 mV from the default value of +15 mV using a step of 5 mV) to make sure the membrane potential returned to rest after spiking.
In contrast to the single-compartment setting (Fig. 1E), the simulated propagating spike waveform was wider for the bEIF model than for the WB model (Fig. 2D). This is because the repolarizing conductance of the bEIF model is static (i.e., independent of the additional axial current), whereas the ionic conductances in the WB model are more dynamically regulated by the membrane potential. The overall dependences of the conduction velocity on the internodal length (Fig. 2E) and the axonal diameter (Fig. 2F) were very similar between the WB and bEIF model, both well fitted by a square root curve. Assuming that the internodal length and the axonal diameter vary in proportion to each other (Hursh, 1939), the combined effect of these square root relationships results in the conduction velocity changing linearly with the size of the myelinated axon (Waxman, 1980).
Both the ceiling for the depolarizing conductance (Fig. 1C) and the after-spike repolarizing conductance in the bEIF model (Fig. 1D) are necessary for simulating stable spike propagation. Because of the exponential dependence of the depolarizing current on the membrane potential (Eq. 3), the membrane potential of the sEIF model quickly blows up to infinity once a spike is initiated (Touboul, 2009). Due to this instantaneous divergence of the membrane potential, the simulated action potential did not properly propagate along the axon when the sEIF model was used instead of the bEIF model (Fig. 2G). The ceiling of the depolarizing conductance in the bEIF model slows down the voltage change near the peak of the action potential, resulting in a more realistic spike shape and stable spike conduction than the sEIF model.
When the repolarization current of the bEIF model was replaced by an abrupt potential reset, propagation of electrical activity was observed, but the simulated waveforms were not uniform across the axonal compartments (Fig. 2H). This is because the potential reset (as in the sEIF model) is equivalent to injecting an enormous negative current within a small time step, which leads to discontinuous changes of the membrane potential. In contrast, the spike-mimicking repolarizing current in the bEIF model does not cause such discontinuous, unstable changes. As suggested by the simulation results for the sEIF model (Fig. 2G), the level of the depolarization ceiling factor (value of AT; Table 3; Fig. 1C) in the bEIF model affects the conduction velocity (Fig. 2I).
Spike conduction in unmyelinated axons
By connecting excitable units (Fig. 3A), spike conduction along an unmyelinated axon can also be simulated (Fig. 3B, for the WB model, Fig. 3C, for the bEIF model). As in myelinated axons, the simulated waveform was slightly wider for the bEIF model than for the WB model. Both models, however, showed similar conduction velocities. Furthermore, the dependence of the simulated conduction velocity on the axon diameter was largely comparable between these two models (Fig. 3D), displaying a square root relationship typical for unmyelinated axons (Waxman, 1980). These simulation results demonstrate that the bEIF model can be used for simulating propagation of action potentials in both myelinated and unmyelinated axons.
Response properties of unmyelinated axon models. A, Schematic drawing (top) and modeled electric circuit (bottom) of an unmyelinated axon with a diameter d. Each nodal compartment of length l has a capacitance cm and (voltage-dependent) resistance rm and is interconnected with neighboring nodes with an axial resistance ra. See Table 5 for the default parameter values. In panels B, C, voltage responses at every 25th node (i.e., each 0.5 mm apart) are shown. Currents were injected intracellularly into the node #50 (gray traces). D, Dependence of simulated conduction velocity u (m/s) on the axon diameter D (μm). The dotted curve shows a square root fit by 
.
Extracellular stimulation
Neuroprosthetic devices usually use extracellular stimulation (Rattay et al., 2003). To test the applicability of the model to prosthetic stimulation, we simulated the responses of the modeled axon to extracellular current injection (Fig. 4A). In contrast to the case of intracellular stimulation, where the extracellular space was assumed to be isopotential, extracellular injection of current produces a gradient of extracellular potential, which is the source of the intracellular axial current along the modeled axon. For both the WB (Fig. 4B) and the bEIF (Fig. 4C) models, an extracellularly injected negative current induced positive responses at the closest node (Fig. 4B,C, small arrows) that lead to spike initiation. Depending on the relative location between the node and the electrode, the response of each node to the extracellular current is either depolarization or hyperpolarization (Fig. 4D), as was demonstrated in earlier modeling studies (Rattay, 1986; Warman et al., 1992). The overall responses, including spike generation and conduction, were largely similar between the two models (Fig. 4B,C), confirming that the bEIF model can be used not only with intracellular current injection, but also with extracellular stimulation.
Responses of myelinated axon models to extracellular stimulation. A, Schematic drawing of a myelinated axon stimulated with extracellular current injection. The extracellular voltage at each node is determined by the distance between the node and the stimulus electrode (see Materials and Methods for the equations). B, Spike conduction along the modeled myelinated axon simulated with the WB model. C, Spike conduction along the modeled myelinated axon simulated with the bEIF model. In panels B, C, voltage responses at every 10th node (i.e., each 2 mm apart) are shown. Gray arrows indicate the intracellular responses caused by the extracellular negative current injection. D, Location-dependent voltage responses to extracellular stimulation. Simulated membrane potentials at the offset of extracellular current stimulation (-1 mA, 0.1 ms) are plotted as a function of the distance from the node #20 (gray traces in B, C), which is the closest node to the stimulus electrode, with a separation of 1 mm (see inset for a schematic drawing).
Application to ANs
As an application of the bEIF model, we simulated the spike initiation and conduction of AN fibers. First, we constructed a single-compartment model of low- and high-frequency ANs (see Materials and Methods for details). Previous physiological measurements showed tonotopic (frequency-dependent) variations of membrane properties in spiral ganglion neurons (AN cells; see Rusznák and Szú́cs, 2009; Davis and Crozier, 2015 for reviews). The difference in input resistance between low- and high-frequency ANs was simply represented as the difference in the leak conductance density GL of the bEIF model (Tables 7, 8), while other physiological parameters were identical between low- and high-frequency AN models. This modification led to a delayed spike initiation for the low-frequency model (Fig. 5A), although the overall voltage trajectory after scaling the time axis was largely indistinguishable between these two models (Fig. 5B). These simulation results are supported by physiological measurements reporting that low- and high-frequency neurons in mice had similar spiking thresholds but that the spike response latency was larger for low-frequency neurons (Adamson et al., 2002). As shown in Figure 1B, the bEIF (and WB) model expresses a Type I spiking behavior (i.e., zero spiking frequency at threshold), contrasting to the standard HH model that is Type II (non-zero spiking frequency at threshold). This Type I response property of the bEIF model enabled us to simulate the observed delayed spike generation in AN cells, which is generally inconceivable with Type II models.
Response properties of AN axon models. See Table 7 for the default parameter values. A, Spike responses of the low-frequency (red) and high-frequency (blue) single-compartment AN models driven by step current inputs. B, Same traces as in A but with a rescaled time axis. C, Spike conduction along the modeled low-frequency myelinated AN axon. D, Spike conduction along the modeled high-frequency myelinated AN axon. In panels C, D, voltage responses at every five nodes are shown.
Next, we adopted the same parameter sets to simulate spike conduction along the central part of myelinated AN fibers (Materials and Methods). The axonal diameter and internodal length were determined from previous anatomic measurements in cats (Liberman and Oliver, 1984). The simulated propagating spike waveform was wider for the low-frequency model (Fig. 5C) than for the high-frequency model (Fig. 5D), reflecting the difference in response latency (Fig. 5A). Calculated conduction velocities were 9.1 and 14.3 m/s for low- and high-frequency models, respectively. These values correspond to the measured velocities in cats (11.6 ± 1.6 m/s, Nguyen et al., 1999; ∼10 m/s, Miller et al., 2004). Our simulation results predict that high-frequency AN fibers should have higher conduction velocity than low-frequency fibers, because of the shorter response latency in the former. Testing this prediction will be a subject of future physiological studies in the field.
Computational time
To compare the computational performances of the models, we calculated the average integration time for a modeled axon (Table 9). In agreement with the reduced number of equations and parameters, the bEIF model was several times faster than the WB model. In addition, the implementation with a compiled language led to a computation several times faster than the MATLAB code, while the computational advantage of the bEIF model was consistent between these implementations. These results roughly correspond to previous reports that compared the computational costs between HH-type and IF-type models (Destexhe, 1997; Ashida et al., 2015) and between C and MATLAB implementations (Goodman and Brette, 2008).
Computational time for calculating axonal spike conduction
Discussion
In this study, we introduced a simple model of nerve spike conduction based on the EIF model (Fig. 1). In comparison to the conventional HH-type model, our bEIF model has much fewer parameters (9 vs 25) and better computational performance (Table 9), but still retains fundamental functions for reproducing action potential propagation along the modeled axon (Figs. 2–4). Application of the model to ANs replicated measured conduction velocity in cats, and predicted that the velocity varies along the tonotopic axis (Fig. 5).
Advantages of simple models
Simple phenomenological models can serve as a practical substitute for complex, conductance-based models, especially when descriptions of detailed ion channels dynamics are not required, when computational simplicity and mathematical transparency are desired, or when only insufficient empirical data are available for constraining a complex model. A simple neuron model was recently adopted, for instance, for a real-time simulation combined with an artificial fingertip (Oddo et al., 2016).
Fundamental lack of relevant biophysical data is a frequent impediment to the development of neural models for prosthetic simulation. To systematically tune the parameters of a neuron model (either IF or HH), measurements of intracellular membrane potentials are generally required (Badel et al., 2008; Rossant et al., 2010; Meliza et al., 2014). Because of ethical and technical limitations, however, measured data from human nerves are usually sparse (if not totally unavailable). Furthermore, a number of studies demonstrated that anatomical and physiological properties of human neurons may differ considerably from those of other animals (Felix, 2002; Angelino and Brenner, 2007; Eyal et al., 2016; Zhang et al., 2017). This makes it even more difficult to extrapolate non-human data to humans (see discussion below for related limitations in prosthetic simulations). Known as the “curse of dimensionality” (Almog and Korngreen, 2016), fitting model parameters of a nonlinear system becomes increasingly troublesome with an increase in the number of unconstrained parameters. Moreover, the geometry of “good” parameter sets may be highly skewed in the high-dimensional parameter space (Marder and Taylor, 2011), leading to a general difficulty in justifying the selection of parameters of complex models with a limited amount of data. The reduced number of parameters in the bEIF model may help mitigate these difficulties, at least when compared to complex HH-type models.
Because of its mathematical simplicity, the IF model and its variations have been used widely in theoretical and computational neuroscience (Koch, 1999; Gerstner et al., 2014, and references therein). Bifurcation analyses, for example, allow direct examination of spiking mechanisms of the EIF model (Touboul and Brette, 2008). Occasionally, IF-type models have also been used in multi-compartment simulations (Rospars and Lánský, 1993; Clopath et al., 2007; Saparov and Schwemmer, 2015; Aspart et al., 2016). However, these models normally have only one spike initiation site to avoid the problem with the instantaneous potential reset (see related discussion by Cessac and Viéville, 2008). In the bEIF model, the potential reset, which led to unstable waveforms of conducting spikes (Fig. 2H), was replaced with repolarizing conductance to replicate the downward trajectory of the spike waveform. Introduction of spike-mimicking current to IF-type models had already been suggested in prior studies (Ashida et al., 2015; Saparov and Schwemmer, 2015).
In early mathematical analyses of spike propagation, FitzHugh–Nagumo-type models were preferred, because they have fewer variables and are thus much easier to analyze than the HH model (Rinzel and Keller, 1973). The FitzHugh–Nagumo model, however, has major drawbacks: its fast activation variable (usually written as V) does not directly correspond to the membrane potential of a real neuron, and the parameters of the model have no clear biological interpretations. Our EIF model-based approach may be useful in alleviating these problems, as its parameters and variables have more intuitive biophysical meanings (e.g., membrane potential, conductances, spike-generating currents, etc.) while keeping a similar level of mathematical complexity as the FitzHugh–Nagumo model.
Disadvantages and limitations
Previous studies have revealed a number of anatomical and physiological specializations in axons, which are nevertheless not always considered in existing axon models, including ours. For example, Nav1, Kv3, and Kv7 channels are clustered at the nodes of Ranvier, while Kv1 channels are distributed at juxtaparanodes under the myelin sheath (Debanne et al., 2011; Freeman et al., 2016; Kim and Rutherford, 2016). With some rare exceptions (McIntyre et al., 2002; Brown and Hamann, 2014), however, models of myelinated axons do not take the detailed distributions of ion channels into account. In our simulations (Fig. 2), spike-generating ionic currents were simply decomposed into depolarizing and repolarizing components without considering their actual ionic compositions. Moreover, histological studies suggested that ion channels are distributed unevenly along the actual AN fiber (Hossain et al., 2005; Yi et al., 2010; Kim and Rutherford, 2016). Hence our naive assumption that the axonal properties match the somatic properties (used in Fig. 5) is likely to be violated. Further refinement of the model would require detailed physiological characterization along each AN fiber and across the tonotopic axis.
The simulated voltage of the bEIF model does not undershoot after an action potential, since its repolarizing current is driven by the leak reversal potential EL. Introducing a different reversal potential (such as EK) could make the simulated spike waveform more realistic and closer to that of the WB model (Fig. 1A), but at the cost of adding another unconstrained parameter to tune. To calculate the repolarization conductance, we used an alpha function solely because of its simplicity, which nevertheless might have to be revised with a different function when a fine tuning of the depolarization phase is important. Moreover, spike shapes can significantly differ between the cell body and the axon (Kole et al., 2007). Further modifications and tuning of the model currents would thus be necessary to improve the physiological plausibility of simulated spikes propagating along the axon.
Possible expansions and applications
To better account for the nonlinear membrane dynamics of a real neuron, a number of modifications of IF-type models have been proposed. Examples include bursting with T-type calcium current (Smith et al., 2000) or with persistent sodium current (Breen et al., 2003); spike-rate accommodation with slowly adapting current (Brette and Gerstner, 2005; Barranca et al., 2014) or with an after-hyperpolarization current (Zhou and Colburn, 2010; Barranca et al., 2014); subthreshold nonlinearity with low-voltage-activated potassium current (Svirskis and Rinzel, 2003; Ashida et al., 2015); and adaptation and stochastic fluctuation of the threshold (Takanen et al., 2016). Similar modifications can be incorporated into the bEIF model. We note, however, that the bEIF model would not fully replace HH-type models, but these two model types should rather complement each other. A user can and should choose an appropriate model according to the intended goals of modeling (Ashida et al., 2017): HH-type models for simulating the detailed ionic dynamics, and IF-type models for more phenomenological, handy description of neuronal spiking behavior.
Random opening of ion channels is suggested to affect neuronal coding properties (Ashida and Kubo, 2010; White et al., 2000; Negm and Bruce, 2014; Moezzi et al., 2016), including nerve conduction (Faisal and Laughlin, 2007). The utility of channel noise in cochlear implants has also been suggested (Rubinstein et al., 1999; White et al., 2000). To fully account for the ion channel stochasticity, Markov channel models in combination with an HH-type membrane equation would be required (Goldwyn and Shea-Brown, 2011). In practice, adding an adequate amount of artificial noise in a model can, at least in part, mimic the stochastic activity of electrically stimulated nerves (Joshi et al., 2017). Similarly, the addition of a noise term to the bEIF model would be necessary for simulating non-deterministic responses of the modeled nerve.
Earlier simulations of myelinated (FitzHugh, 1962; Brill et al., 1977; Moore et al., 1978) or unmyelinated (Cooley and Dodge, 1966) axons focused primarily on the biophysical mechanisms of nerve conduction. More recent modeling approaches aimed to gain clinical and engineering implications. Examination of degraded spike conductions caused by demyelination is one such example (Coggan et al., 2010; Brown and Hamann, 2014; Resnick et al., 2018). Moreover, driven by the rapid development of implantable devices that use electrical pulses to restore the functions of peripheral and central nerves (for reviews, see Masani and Popovic, 2011; Eiber et al., 2013), modeling approaches to simulate the activity patterns of electrically stimulated nerves have become an important tool for evaluating and predicting the performance of these prostheses (Rattay et al., 2003; also see Craver, 2010 for related philosophical considerations).
Recent prosthetic simulations using conductance-based (or other related) models include electrical stimulation of the retina (Barriga-Rivera et al., 2017), ANs (Briaire and Frijns, 2005; Negm and Bruce, 2014; O'Brien and Rubinstein, 2016; Joshi et al., 2017; Nogueira and Ashida, 2018), motor nerves (ElBasiouny and Mushahwar, 2007), and spinal and other peripheral nerves (Ladenbauer et al., 2010; Raspopovic et al., 2011; Capogrosso et al., 2013; Kent and Grill, 2013). A modeling approach in the cochlear implant study, for example, simulated several tens of thousands of excitable nodes distributed in three-dimensional space to predict the aggregated electrical response of the tissue (Nogueira et al., 2016). Such large-scale modeling studies generally require efficient phenomenological descriptions of neuronal spiking activity. Furthermore, as discussed above, simulations of human nerves often suffer from the lack of relevant physiological data. Our bEIF model thus offers a computationally efficient alternative to more complex (e.g., HH-type) models to be used in future prosthetic simulations.
Supplementary 1
Acknowledgments
Acknowledgements: We thank the two anonymous reviewers for their constructive comments and suggestions. Language services were provided by www.stels-ol.de.
Footnotes
The authors declare no competing financial interests.
This research was supported by the Deutsche Forschungsgemeinschaft(DFG) Cluster of Excellence EXC 1077/1 “Hearing4all” (G.A. and W.N.).
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
References
In this issue
