Membrane capacitive memory alters spiking in neurons described by the fractional-order Hodgkin-Huxley model.

Excitable cells and cell membranes are often modeled by the simple yet elegant parallel resistor-capacitor circuit. However, studies have shown that the passive properties of membranes may be more appropriately modeled with a non-ideal capacitor, in which the current-voltage relationship is given by a fractional-order derivative. Fractional-order membrane potential dynamics introduce capacitive memory effects, i.e., dynamics are influenced by a weighted sum of the membrane potential prior history. However, it is not clear to what extent fractional-order dynamics may alter the properties of active excitable cells. In this study, we investigate the spiking properties of the neuronal membrane patch, nerve axon, and neural networks described by the fractional-order Hodgkin-Huxley neuron model. We find that in the membrane patch model, as fractional-order decreases, i.e., a greater influence of membrane potential memory, peak sodium and potassium currents are altered, and spike frequency and amplitude are generally reduced. In the nerve axon, the velocity of spike propagation increases as fractional-order decreases, while in a neural network, electrical activity is more likely to cease for smaller fractional-order. Importantly, we demonstrate that the modulation of the peak ionic currents that occurs for reduced fractional-order alone fails to reproduce many of the key alterations in spiking properties, suggesting that membrane capacitive memory and fractional-order membrane potential dynamics are important and necessary to reproduce neuronal electrical activity.

Introduction

The properties of excitable cells and cell membranes have been studied for over a century, dating back to the studies of Weiss [1], Lapicque [2], and Nernst [3] in the early 1900s. These early studies established the now well-known inverse relationship between the strength of a stimulus necessary to elicit a spike (action potential) and the stimulus duration. Theoretical models were developed to explain the empirically-derived strength-duration curves: Representing the passive cell membrane as a parallel resistor-capacitor circuit, assuming ideal capacitive behavior, the membrane voltage V dynamics are given by where C is membrane capacitance, R is membrane resistance, and I represents the applied stimulus current. This simple yet elegant membrane representation is the foundation for essentially all excitable cell models, originating with the Hodgkin-Huxley model [4], in which ionic currents are represented by time- and voltage-dependent conductances, in parallel with the membrane capacitance. Although minimal, the passive membrane model can predict a wide range of sub-threshold and threshold behavior. For example, if the applied current is a current step at time t = 0, i.e., I (t) = I u(t) (where u(t) is the standard Heaviside function), then the solution of this ordinary differential equation, where membrane time constant τ = R C and the membrane is assumed to be initially uncharged (V (0) = 0), can be rearranged to solve for the threshold current I necessary to elicit a spike for a stimulus of duration d, i.e., the strength-duration relationship, where rheobase current I = V /R is the current threshold for an infinite duration stimulus and V is the membrane potential threshold. Eq 3 accurately reproduces the asymptotic behavior for long duration stimuli in the strength-duration curve observed in experiments. However, for short duration stimuli (d ≪ τ), Eq 3 is, to first-order, equivalent to an inverse scaling, and early studies showed that experimental data did not well-fit this simple inverse scaling, d −1. Rather, a general power law relation, d −, was more appropriate, with α ranging from 0.5 to 1 [5]. This was shown as early as 1933, when Cole notes Lapicque’s experimental data fits an average value of α = 0.656, while Ruston’s data fit values of α = 0.76 and 0.86 for warm and cold frog sciatic nerve, respectively [5].

Westerlund and Ekstram showed that Jacques Curie’s 1889 empirical law for current through capacitors and dielectrics [6] can be used to derive the capacitive current-voltage relationship for a non-ideal capacitor [7], where 0 < α < 1, is a fractional-order capacitance with units (amp/volt)⋅sec, and, for now, the exact definition of the fractional-order derivative d V /dt will be left ambiguous but will be clarified shortly.

Magin showed that incorporating this non-ideal capacitor into a fractional-order membrane model can reproduce the scaling of the strength-duration curve at both short and long durations [8]. Specifically, the fractional passive membrane dynamics are given by the fractional-order differential equation (cf. Eq 1) To our knowledge, the physiological source of such non-ideal capacitive behavior is not known, but we speculate may arise due to heterogeneities in the dielectric properties of the membrane, or more generally what is termed “capacitance dispersion” in electroanalytical chemistry [9]. As shown in S1 Text, Eq 6 can be solved via the Laplace transformation, with an analytical solution given by (cf. Eq 2) where , E is the two-parameter Mittag-Leffler function, which generalizes the exponential function, given by and Γ(x) is the gamma function. Note that for α = 1, E 1,1(z) = exp(z), and Eqs 2 and 7 are equivalent, as expected. Eq 7 can be rearranged, as before (cf. Eq 3), and, importantly, yield a strength-duration relation that follows a general power law for short duration stimuli and asymptotic behavior for long duration stimuli (Fig 1A).

Fig 1

Properties of the fractional-order passive membrane.

(A) Strength-duration curves, derived from the fractional passive membrane, are shown as a function of fractional-order α (Eq 9). (B) The magnitude (top) and phase (bottom) of the complex impedance of the fractional-order passive membrane are shown as a function of the normalized frequency ωτ, for different values of α. (C) The normalized membrane potential V /(I R ) response following a current step is shown as a function of normalized time t/τ on a linear (top) and logarithmic (bottom scale), for different values of α.

Despite these early studies which suggest that a non-ideal (or fractional-order) capacitive element may be more appropriate to represent passive membrane dynamics, essentially all excitable cell models assume an ideal capacitor, i.e., α = 1. In recent years, fractional-order dynamics have been shown to provide an improved description of many biological phenomena, including mechanical properties of viscoelastic tissue [10], the tissue-electrode interface [11], pharmacokinetics of drug delivery and absorption [12-14], and anomalous calcium subdiffusion in microdomains [15]. While it is clear that fractional-order membrane potential V dynamics will alter the passive response to sub-threshold stimuli, it is not obvious if fractional-order dynamics will alter the properties of the active neuron and response to super-threshold stimuli, given the bi-directional coupling between V and ionic currents. In this study, we will investigate fractional-order V dynamics in the classical Hodgkin-Huxley model [4] and test the hypothesis that fractional-order V dynamics influence the properties of the active neuron across multiple spatial scales. In this study, we provide a brief background on the mathematics of fractional calculus and elaborate on and illustrate some properties of the fractional-order passive membrane, introduce the fractional-order Hodgkin-Huxley model, and characterize properties of neuronal spikes. We study the properties of sub-threshold and spike propagation in a nerve axon, characterizing the passive and active cable, respectively. Finally, we study the dynamical properties of a network of fractional-order neurons. We conclude with a discussion of our results.

Background on fractional calculus

Definitions of the fractional-order derivative

Our brief background on fractional calculus presented here is by no means complete. For further details, methods, and applications of fractional derivatives and differential equations, we suggest the monographs by Podlubny [16] and Magin [17]. The definition of a fractional-order derivative is not immediately obvious, and indeed, multiple definitions exist. The classical definition of the fractional derivative of order α, known as the Riemann-Liouville fractional derivative, denoted and given by where m−1 < α < m (i.e., m is the first positive integer larger than α), follows from the definition of the fractional integral.

Note that the subscripts 0 and t on indicate limits for the fractional integral and highlight the fact that, in contrast with integer-order derivatives, the fractional-order derivative depends on the previous history of the function and is not strictly a local property. An alternative definition was given by Caputo as where y ((t) is the m -order derivative with respect to time. Note that both definitions involve fractional-order integration and integer-order differentiation, with the order of the two operations different. However, these definitions are not equivalent, except in the case that all initial conditions are equal to 0. This is analogous to the integer-order operators, integration and differentiation, not being commutative and differing by a constant. In this case, the difference is a function of time. The appropriate choice of representation will depend on the situation being modeled. However, the Caputo definition is often preferable for modeling physical processes, since initial conditions of fractional differential equations using the Caputo definition are specified in terms of integer-order derivatives. In contrast, fractional differential equations using the Riemann-Liouville definition requires initial conditions specified in terms of fractional-order derivatives, which typically have vague physical interpretation.

Numerical methods

We consider one more definition for the fractional-order derivative, the Grunwald-Letnikov fractional derivative, which will lead to a simple numerical method for integrating a fractional differential equation [18]. Consider the well-known definition of the first-order derivative: Applying this definition again gives us the second-order derivative, and, again, the third derivative, By induction, the general n -derivative is given by where is the usual notation for the binomial coefficients, and the second equality follows from the gamma function relation Γ(n) = (n−1)!. The use of the gamma function allows one to extend the derivative definition from integer-order values of n to arbitrary fractional-order.

It can be shown that the Grunwald-Letnikov definition is equivalent to the Riemann-Liouville definition [16]. In this study, we will consider the fractional-order Hodgkin-Huxley model initially at rest (V (0) = 0), and thus all three definitions are equivalent. To illustrate the numerical integration scheme, consider a generic fractional differential equation with zero-initial conditions Applying a finite-difference scheme to Eq 13 using Eq 12d, we have where Rearranging Eq 14a, an explicit numerical scheme is given by where the terms can be computed recursively using and . Thus, the terms represent a “weighting” of previous history or memory on the current state of the system, a consequence of the non-locality of the fractional derivative. Values of are shown in S1 Text for different values of α (Fig. A in S1 Text). For α = 1, all terms beyond k = 1 are equal to 0, and Eq 15 is equivalent to the forward Euler method. Importantly, decreases as k decreases, as the current value in the numerical integration depends increasingly less on earlier system states. However, as α decreases, decreases more slowly, i.e., there is a larger dependence on the system history. Note that discretization of the Grunwald-Letnikov definition can be used to solve systems with non-zero initial conditions by including an additional term in the numerical scheme above [18].

Fractional-order neuron model

Fractional-order passive membrane properties

Before investigating the influence of fractional-order α on the properties of the active neuron, it is instructive to consider the response of the fractional-order passive membrane. We first consider the passive response to a sinusoidal current input of varying frequencies, which allows for standard circuit analysis using complex impedances. For the standard passive membrane patch, the resistor complex impedance is strictly real, 𝓩 = R . An ideal capacitor complex impedance is strictly imaginary, , while a non-ideal capacitor complex impedance is given by , where α = 1 refers to an ideal capacitor and α = 0 refers to a pure resistor. Such a non-ideal capacitor is often called a constant phase element. The magnitude and phase of the total impedance for the passive membrane, i.e. the parallel resistor-capacitor circuit, are shown as a function of frequency ω for different values of α in Fig 1B. For all values of α, the total impedance magnitude approaches 1 at low frequencies and decreases at high frequencies, i.e., the membrane acts as a low-pass filter. However, for small α, the drop-off at high frequencies is reduced, i.e., the slope (on a log-log scale) is smaller, reflecting the power law relation present at short time scales. Further, the phase angle observed for high frequencies varies with α, as the phase approaches −90° for α = 1 and −45° for α = 0.5. From Eq 16, it can be shown that the phase at high frequencies is given by . Early studies by Cole and colleagues measured phases corresponding to α values between 0.61 and 0.86, consistent with measurements based on strength-duration data [19-21].

It is worth emphasizing that the order α of the fractional passive membrane does not simply alter the effective time constant, but rather the dynamics at multiple time scales. To illustrate, we consider the response of the fractional passive membrane following a step current input, to characterize the sub-threshold response in a neuron or other excitable cell. The membrane potential V (t) time course, given by Eq 7, is shown as a function of time for different values of α in Fig 1C. For small α, the membrane initially charges faster, compared with the exponential solution when α = 1, following the current step input. Subsequently, V (t) approaches the steady-state voltage more slowly. Thus, in this study, we investigate to what extent the sub-threshold dependence on α will also influence the super-threshold response and the properties of the active neuron and interacting neurons.

Fractional-order membrane patch model

With sufficient numerical tools to proceed, we present here the fractional-order Hodgkin-Huxley (fHH) model of the neuron: where v = V −V represents the membrane potential V relative to the resting potential V , m, h, and n are the sodium current I activation, I inactivation, and potassium current I activation gating variables, respectively, I(t) is a time-dependent applied stimulus, the prime denotes a first-order derivative with respect to time, and d v/dt refers to the fractional derivative in the Caputo sense (since we specify integer-order initial conditions).

Numerical integration

We can use Eq 15 to integrate the entire system (Eq 17) using a specified value for α for integration of Eq 17a and a value of 1 for Eqs 17b–17d. In practice, summation of the previous history becomes computational expensive for long simulations, and thus integration of Eqs 17b–17d was performed via forward Euler coupled with the Grunwald-Letnikov integration scheme in MATLAB (Mathworks, Inc.), using of a time step of Δt = 10−4−10−3 ms. As α decreases, we found that a smaller integration time step Δt is required for numerical stability, which leads to significantly longer simulation times. Maximum time steps and simulation times are shown as functions of α in S1 Text (Fig. B in S1 Text). Standard parameters and initial conditions (Table A in S1 Text), and gating equations are given in S1 Text, and simulation code is provided in S1 Code.

Numerical integration using the standard forward Euler method, i.e., for α = 1, the membrane potential v at time t , v(t ), is a function of state variables at time t , whereas for 0 < α < 1, v(t ) is a function of all of prior history of the membrane potential due to capacitive memory. Using the numerical scheme in Eq 15, we can directly characterize the influence of capacitive memory by defining a voltage memory trace, following the similar approach of Teka et al. [22], that is the weighted sum of previous values of membrane potential, excluding the immediate previous term, where we define v (t 0) = v (t 1) = 0.

Spikes triggered by a brief stimulus pulse

Spike properties

We first characterize the influence of the fractional-order α on the properties of the neuronal spike in response to a brief stimulus pulse (Fig 2). For small α, V is depolarized rapidly, in fact faster than sodium current (I ) activation, such that V subsequently repolarizes until sufficient time has passed that I activation triggers a full depolarization and the spike upstroke (Fig 2A, red and green traces). As a consequence of the rapid membrane depolarization, the spike peak is earlier for small α. However, once triggered, there are only small differences in the properties of the spike and ionic currents. Sodium current I exhibits the largest dependence on α: as α decreases, the peak I magnitude increases (Fig 2B). The I time course is also slightly modified. We found that the increase in I peak magnitude was primarily due to a decrease in I inactivation, i.e., the I inactivation gate at the time of the current peak, h , increases as α decreases.

Fig 2

Properties of the fractional-order Hodgkin-Huxley spike.

(A) The membrane potential V , sodium current I , potassium current I , and voltage memory trace v are shown as a function of time, for different values of fractional-order α. (B). V maximum and minimum (left), I and I peak current magnitude, and h (the sodium inactivation gating variable at the time of peak I current) are shown as a function of α. Spikes are elicited by a brief 0.1-ms duration, 1.5x threshold stimulus.

The magnitude of the voltage memory trace v greatly increases and is non-zero over a longer time period as α decreases, demonstrating a larger influence on membrane potential dynamics as expected. Less prominent effects include a small decrease in the potassium current (I ) peak magnitude, the V maximum (spike peak), and the V hyperpolarization overshoot (an increase in V minimum), such that the spike amplitude is reduced (Fig 2B).

Post-spike refractoriness

To characterize refractoriness following a spike, we measure the minimum time period between stimuli for which two spikes could be triggered by a brief 0.1-ms stimulus (Fig 3). We found that this minimum time decreased as α decreased. However, after accounting for the earlier spike peak (as seen in Fig 2A), the difference in refractoriness as α varied was small: the time between the spike peak and the subsequent stimuli varied less than 1 ms for α between 0.4 and 1. This is likely due to the fact that refractoriness is primarily a consequence of the sodium current I recovery from inactivation, which is minimally influenced by the fractional-order when the stimulus is brief.

Fig 3

Refractoriness in the fractional-order Hodgkin-Huxley model.

The minimum time period between stimuli (top), the time to the spike peak (middle), and their difference (bottom) are shown as a function of fractional-order α. Spikes are elicited by a brief 0.1-ms duration, 1.5x threshold stimulus.

Spikes triggered by a constant stimulus

Repetitive firing in the fractional-order neuron model

It is well-established in the classical Hodgkin-Huxley model (i.e., α = 1), the neuron will repetitively fire or spike in response to a constant applied stimulus I ∈ [I 1, I 2], where I 1 and I 2 represent critical values discussed below [23]. We next investigate to what extent fractional-order α alters the properties of repetitive firing for a constant stimulus. We show the membrane potential V , ionic currents I and I , and voltage memory trace v as functions of time for a constant applied current, for different values of I and α, during a 100-ms duration simulation (Fig 4A).

Fig 4

Repetitive firing in the fractional-order Hodgkin-Huxley model.

(A) The membrane potential V , sodium current I , potassium current I , and voltage memory trace v are shown as a function of time in response to a constant applied current, I = 20 (left), 100 (middle), and 140 (right) μA/cm2, for different values of fractional-order α. (B) The instantaneous spike frequency is shown as a function of the interspike interval (ISI) number for different values of α.

We found a complex time-varying relationship between α, I , and the spike frequency and amplitude. For small values of I (Fig 4, left), spike frequency clearly decreases as α decreases. Spike amplitude and peak I are reduced and peak v is greatly increased, while peak I is slightly increased. The instantaneous spike frequency increases as a function of the interspike interval (ISI) number for all values of α and approaches an asymptotic value after 2–3 intervals (Fig 4B). For a given value of α, the increase between the first and final interval is small, approximately 5 Hz.

For an intermediate value of I (Fig 4, middle), spike frequency decreases to a much smaller extent as α decreases. For all values of α, spike amplitude is reduced, such that the increase in v as α decreases is mitigated. Both peak ionic current magnitudes decrease as α decreases. As before, the instantaneous spike frequency increases as a function of ISI number, approaching an asymptotic value after 4–5 intervals. The increase between the first and final interval is larger, approximately 20 Hz. For a larger value of I (Fig 4, right), for α = 1, small amplitude oscillations or spikes persist, following stimulus onset. However, as α decreases, spike amplitude decreases to 0, such that the neuron no longer fires, and V is held at an elevated level, a phenomenon known as excitation block.

The mechanisms underlying the dependence of spiking properties and excitation block on α are complex and can be generally explained as follows: Following stimulus onset, although peak I increases during the first spike, rapid membrane polarization reduces the time available for I recovery from inactivation following the first triggered spike, which reduces peak I current and peak membrane potential V during the subsequent spike. This, in turn, reduces I activation and peak I current. Reduced I current causes a less hyperpolarized V overshoot, which reduces the magnitude of leak current I , which is depolarizing at this phase of the spike. Reduced I decreases the rate of interspike V depolarization and thus subsequently reduces spike frequency. Capacitive memory further reduces the V depolarization rate, as the voltage memory trace v more heavily weighs the hyperpolarized V values during the preceding overshoot. For larger values of I , over several spikes, the reduced peak V also reduces I activation, which in turn further reduces I peak current. These feedback interactions are self-limited, and spike and ionic current characteristics stabilize after several spikes, resulting in spike frequency and amplitude approaching asymptotic values. By averaging over time periods that include the neuron at rest and the hyperpolarization overshoot, the voltage memory trace v effectively “dampens” changes in the membrane potential. Even though the magnitude of v decreases as I increases, the magnitude of the ionic currents are decreased to a larger extent, such that capacitive memory is more influential for large values of I . For a larger value of I , small fractional-order α may reduce the spike amplitude sufficiently and result in excitation block.

To further determine if the dependence of spiking properties on α is primarily a consequence of the fractional passive membrane dynamics or the subsequent modulation the ionic currents, both passive and active, we run simulations in which we assume first-order V dynamics but also scale the ionic current I , I , and I conductances, g , g , and g , respectively, such that the peak current magnitudes are equivalent to the values for a given α in the fractional-order model, i.e., the asymptotic current peak I and I values in Fig 4A. Peak I magnitude is measured as the maximum depolarizing current, since I is both depolarizing and hyperpolarizing. For example, as α is decreased in this scaled ionic conductance first-order model, g is decreased to account for the smaller peak I magnitude. For I = 20 μA/cm2, g is slightly increased, while for I = 100 μA/cm2, g is decreased. We investigate the influence of scaling the ionic current conductances both individually and combined (Fig 5). Scaling the leak conductance g had minimal influence on spike frequency or amplitude for all values of I , due to its small amplitude relative to the other currents (Fig 5A). For small I , scaling the sodium conductance g had minimal influence (Fig 5B). For intermediate I , scaling g does reduce spike frequency and amplitude, however to a larger extent than observed in the fractional-order neuronal model. For a larger I , scaling g also leads to excitation block, however with a different dependence on α as in the fractional-order model (also see Fig 6B and 6C). For a given value of I , scaling g does not alter spike amplitude and further, increases spike frequency, the opposite effect as observed in the fractional-order model (Fig 5C). When all of the conductances are scaled, the collective influence is that spike frequency increases, not decreases, as the value of α decreases (Fig 5D), demonstrating that simply scaling the magnitude of ionic currents does not reproduce the influence of fractional-order V dynamics. Interestingly, scaling all three conductances does fairly accurately reproduce the reduction in spike amplitude as α decreases.

Fig 5

Repetitive firing in the first-order Hodgkin-Huxley model with scaled ionic current conductances.

The membrane potential V is shown as a function of time in response to a constant applied current, I = 20 (left), 100 (middle), and 140 (right) μA/cm2, for different values of fractional-order α. Sodium, potassium, and leak conductances, g , g , and g , respectively, are scaled, individually (A-C) and combined (D, top), such that peak current measurements are equivalent to values for particular value of α, as described in the text. (D, bottom) The instantaneous spike frequency is shown as a function of the interspike interval (ISI) number for different values of α.

Fig 6

Spike properties in the fractional-order Hodgkin-Huxley model.

(A) Bifurcation diagram of V , sodium current I , and potassium current I , showing steady-state values and limit cycle maximum and minimum, as a function of the applied current I , for different values of α. (B) The critical values denoting I lower and upper limits for spiking (Hopf bifurcations), I 1 and I 2, respectively, are indicated (fractional-order Hodgkin-Huxley model (fHH), solid lines). (C) The spike frequency (top) and amplitude (bottom) are shown as a function of I and α. In B and C, values for I 1, I 2, and spike frequency and amplitude are shown for the first-order model with scaled conductances for comparison (dashed lines, see Fig 5 and main text for more details). In the bottom panel of C, the solid and dashed lines are nearly identical.

In Fig 4, we demonstrate that spike amplitude and frequency approach asymptotic values after several spikes during a constant stimulus. Therefore, we next investigated how these asymptotic values varied for different values of α and I . In the classical Hodgkin-Huxley model, as the applied current magnitude I increases, the resting potential gradually increases until I = I 1 (a subcritical Hopf bifurcation), a critical value at which the neuron fires repetitively. As I increases further, spike frequency increases while the spike amplitude decreases, until I = I 2 (a supercritical Hopf bifurcation), a critical value for excitation block.

The value of α will not influence the steady-state of the system; however α can alter the system stability [24] and thus may influence the values for I 1 and I 2 and the range of I values for which the neuron spikes. In Fig 6A, we show the steady state for a resting (stable) system and maximum and minimum values for a spiking (unstable) system, for V , and ionic currents I and I , as a function of I . As α decreases, the V range is reduced (consistent with the reduced spike amplitude in Fig 4), in conjunction with a reduced range for I and I . The reduction in I is primarily due to a reduction in the sodium activation gate m, while the reduction in I is primarily due to a reduction in the V driving force, with a smaller influence via reduced potassium activation gate n. There is minimal influence of the values on sodium inactivation gate h (Fig. C in S1 Text).

In Fig 6B (solid lines), we show I 1 and I 2 as a function of α. I 1 and I 2 increase and decrease, respectively, as α decreases, such that the I range for which the neuron spikes is reduced. In Fig 6C (solid lines), the asymptotic values of the spiking frequency (top) and spike amplitude (bottom), calculated over the final 30-ms of a 100-ms simulation, are shown as a function of α and I . In the original Hodgkin-Huxley model (α = 1), as I increases, spike frequency increases, while spike amplitude decreases (Fig 6C, black lines). We found this relationship—increasing spike frequency and decreasing spike amplitude, for increasing I —for all values of α (Fig 6C, solid color lines). As observed in Fig 4, for all values of I , spike amplitude decrease as α decreases. For small values of I , spike frequency decreases as α decreases. The reduction in spike frequency is mitigated for larger values of I and, in fact, spike frequency increases as α decreases for I near the critical bifurcation value I 2. We show for comparison (dashed lines in Fig 6B and 6C) measurements of I 1, I 2, and spike frequency and amplitude for the first-order model with all ionic current conductances scaled (Fig 5D). Critical value I 1 is similar in the two models, while I 2 is smaller compared with measurements from the fractional-order model, demonstrating that only accounting for the influence of peak ionic currents underestimates the range of stimuli that result in neuronal spiking. As observed in Fig 5D, the reduction in spike amplitude for decreasing α is in close agreement between the two models, as the dashed and solid lines are nearly identical. However, spike frequency in the scaled conductance first-order model increases as α decreases, in contrast with the fractional-order model.

Fractional-order nerve axon model

Expanding in scale beyond the spatially-clamped membrane patch, we next consider a spatially-extended system, a one-dimensional cable, to investigate whether the fractional-order α influences the properties of electrical propagation in a nerve axon.

Sub-threshold Propagation

As in the previous section, we first investigate the properties of the passive membrane. The fractional-order passive cable equation is given by Eq 6 with the additional of a voltage diffusion term: where g is a longitudinal cable conductance. Eq 19 can be written in standard form, where time constant and space constant . The membrane potential V (x, t) can be determined by the convolution of the applied current I(x, t) and G(x, t), the impulse response, scaled by R . We solve for G(x, t) using an analytical-numerical approach using the Laplace-Fourier transform (see S1 Text), shown in Fig 7A and 7B. At early time points, for small α, G(x, t) is more “spread out” in space, while G(x, t) is less spread out at later time points (Fig 7A). At the site of the impulse (x = 0) and one length constant away (x = λ), the impulse is also more spread out in time for small α (Fig 7B).

Fig 7

Sub-threshold impulse and voltage response in the passive fractional-order cable equation.

(A) The impulse response function G(x/λ, t/τ) is shown as a function of space x, normalized by space constant λ, at times t = 0.05τ and t = τ, where τ is the time constant, on a linear (top) and logarithmic (bottom) scale, for different values of fractional-order α. (B) G(x/λ, t/τ) is shown as a function of normalized time t/τ at location x = 0 (top) and x = λ (bottom). (C) The normalized voltage response to a current step input at the origin x = 0 is shown as a function of normalized time t/τ at locations x = 0 and x = λ. The voltage response in the membrane patch is shown for comparison (dashed lines, Fig 1C). (D) The normalized position of stimulus propagation x/λ is shown as a function of normalized time t/τ (the time at which the normalized voltage response is 0.5) (E) The pseudo-velocity, given by the slope of the stimulus propagation, in units of λ/τ, is shown as a function of α.

As with the membrane patch, we plot the voltage response to a current step input at the origin (x = 0), that is, I(x, t) = I u(t)δ(x), where δ(x) is the Dirac delta function. While it is not possible to characterize a “velocity” of the sub-threshold response, since V decays to 0 as x approaches infinity, we can measure the propagation of V (x, t), normalized by the steady-state voltage at that particular spatial location V (x) (Fig 7C). Membrane polarization is more rapid for small α at x = 0, as in the membrane patch (Fig 1C, shown here for comparison as dashed lines). Downstream of the current input site, at x = λ, V increases at a much earlier time for small α, compared with α = 1, that is, the sub-threshold stimulus propagates down the cable faster for small α. The position of the stimulus propagation is shown as a function of the time when the normalized voltage response reaches 0.5 (Fig 7D). For small α, the stimulus propagation is much faster, and we find that the “pseudo-velocity,” given by the slope of the stimulus propagation position as a function of time, increases as α decreases (Fig 7E). For α = 1, the pseudo-velocity is approximately 2λ/τ, consistent a prior analytical calculation for the classical cable equation [25].

Spike propagation

We next consider the active nerve axon and investigate how fractional-order α influences spike propagation. As in the previous section, we investigate the response of the neuronal model to a brief stimulus pulse and a constant stimulus. Importantly, simulating the response to both a brief and constant stimulus enables the analysis of fractional-order α when capacitive memory may influence spike propagation during a single spike following a resting state and multiple spikes, respectively. The fractional-order Hodgkin-Huxley model can be extended to a cable model by Eq 17a with the addition of a voltage diffusion term: where g is the axon conductance. The gating variable dynamics are similarly functions of space and time. In practice, we integrate this system by converting the fractional-order partial differential equation in Eq 22 to a system of fractional-order ordinary differential equations using the method of lines, with a spatial discretization of Δx = 0.5 mm and utilizing the numerical scheme described above (Eq 15). Simulation code is provided in S1 Code.

Spike propagation trigged by a brief stimulus pulse

We first consider the influence of fractional-order α on spike propagation following a brief stimulus pulse. Representative space-time plots of spike propagation for different values of α are shown in Fig 8A. We find that the spike propagates the length of the cable (1 cm) faster as α decreases. To determine the mechanism underlying faster spike propagation in the cable, we first measure the peak current magnitudes and voltage memory trace along the length of the cable for different values of α (Fig 8B). As in the membrane patch (Fig 2), peak I and v increase and I decreases, respectively, as α decreases. Propagation velocity is measured by calculating the difference between the spike peak time at x = 0.25 and 0.75 cm. We plot the spike propagation velocity as a function of α for different axon conductance levels g (Fig 8C, solid lines). For both values of g, the velocity increases as α decreases, consistent with the passive cable; however, the dependence of the velocity on α is mitigated as g decreases (solid lines, an increase of ∼65% for g = 0.706 μS and ∼88% for g = 7.06 μS, from α = 1 to 0.4).

Fig 8

Spike propagation in the fractional-order Hodgkin-Huxley nerve axon following a brief stimulus pulse.

(A) A space-time plot of the membrane potential V (x, t) is shown for different values of fractional-order α. (B) The peak sodium current I , potassium current I , leak current I , and voltage memory trace v magnitude are shown as a function of position along the cable x, for different values of α. (C) Spike propagation velocity is shown in the fractional-order Hodgkin-Huxley (fHH) nerve axon, as a function of α, for different values of longitudinal conductance g (solid lines). Velocity measurements are also shown (dashed lines) for simulations in which the sodium g and potassium g conductances are scaled such that peak current measurements are equivalent to values for particular value of α. See text for more details. In A and B, axon conductance g = 7.06 μS. Propagating spikes are elicited by a brief 0.1-ms duration, 500-μA/cm2 stimulus at x = 0.

To determine if the dependence of spike propagation velocity on α is primarily a consequence of the properties of the fractional passive cable or the changes in the passive and active ionic currents, as in the previous section, we run simulations in which we assume first-order V dynamics but also scale the ionic current I , I , and I conductances, g , g , and g , respectively, such that the peak current magnitudes are equivalent to the values for a given α and location x (the values in Fig 8B), i.e., we increase g and decrease g and g . In these simulations, spike propagation velocity marginally increases as α decreases (Fig 8C, dashed lines, both increase ∼3% from α = 1 to 0.4) but not nearly to the extent as in the fractional-order nerve axon model. This result, in conjunction with analysis of the sub-threshold stimulus propagation, suggests that the increase in velocity observed for small α is primarily a consequence of the passive cable properties of the axon (Fig 7), and not the modulation of the ionic currents.

Spike propagation trigged by a constant stimulus

We next investigate the influence of fractional-order α on spike propagation during a constant stimulus, in which capacitive memory may alter propagation over successive spikes. In Fig 9A (top), we plot spike propagation velocity as a function of the spike number for different values of α and the applied current amplitude I . In general, spike propagation velocity decreases as a function of the spike number, and in most cases, approaches an asymptotic value after 3–4 spikes. For a small value of I (left), velocity decreases to a small extent as spike number increases, and this decrease is approximately the same for all values of α. We characterize this decrease by calculating the difference between the velocity for a given spike number and between the final (asymptotic) velocity, as a percentage of the final velocity (Fig 9A, bottom). This small change in velocity is due to the fact that, as in the membrane patch, spike frequency is small for small I (Fig 6), such that the time between propagating spikes in the nerve axon is sufficiently long for sodium channel recovery (see Fig. D in S1 Text for further analysis of the relationship between spike frequency and propagation velocity). However, for α = 0.4, propagation fails after a single spike, since the sodium channel recovery is insufficient to maintain the fast propagation velocity over many spikes.

Fig 9

Spike propagation in the fractional-order Hodgkin-Huxley nerve axon during a constant stimulus.

(A) Spike propagation velocity (top) and the change in velocity, as a percentage of the final velocity (bottom), are shown as functions of spike number, for different values of fractional-order α and applied current amplitude I . (B) Velocity measurements are shown for simulations in which the sodium, potassium, and leak conductances, g , g , and g , respectively, are scaled such that peak current measurements are equivalent to values for particular value of α and location x, as described in the text. Axon conductance g = 7.06 μS. Propagating spikes are elicited by a constant stimulus at x = 0.

For an intermediate I (Fig 9, middle, top), velocity decreases to a larger extent as spike number increases, which is a consequence of faster spike frequency and thus less recovery time between spikes. The decrease in the velocity is also larger for smaller α (Fig 9, middle, bottom). For larger I (Fig 9, right, top), we find that propagation fails after 3 and 4 propagating spikes for α = 0.4 and 0.6, respectively, demonstrating that, as in the membrane patch (Fig 6), the range of I values for which the axon can sustain repetitive spiking and propagation decreases as α decreases.

Interestingly, we note that for all values of I , the first spike propagation velocity is equal to the velocity following the brief stimulus pulse (Fig 8), demonstrating that capacitive memory does alter spike propagation progressively over several spikes. If we scale the ionic current conductances as before, with the peak I , I , and I magnitudes measured over the final 30-ms of a 100-ms simulation (Fig 9B), for most parameters, there is only negligible alteration of the spike propagation velocity. Collectively, these results suggest that the influence of fractional-order V dynamics on spike propagation over many spikes is due to fractional-order membrane dynamics and capacitive memory, while the modulation of ionic currents is minimally influential.

Fractional-order neural network model

Finally, we consider a network of N = 50 randomly connected fractional-order Hodgkin-Huxley neurons. We add a synaptic current to the fractional-order Hodgkin-Huxley model (Eq 17a), such that the dynamics of the i neuron are governed by the following equation: where the synaptic current is given by the sum of the excitatory and inhibitory synaptic currents, where S and S are the set of presynaptic neurons with connections to neuron i, with excitatory and inhibitory, respectively, synapses, and s is the gating variable for the postsynaptic conductance, an instantaneous, sigmoidal function of the presynaptic cell potential v with a threshold V [26], that is

Synaptic current parameters are given in S1 Text (Table B in S1 Text), and simulation code is provided in S1 Code. We have previously studied this type of random neural network [27]. Briefly, synaptic connections were determined as follows: The number of presynaptic connections to the i neuron is drawn from a Gaussian distribution with mean μ = 25 and standard deviation σ = 2.5, rounded to the nearest whole number. The presynaptic neuron indices j are chosen at random from integers [1, N]. The type of each synapse, excitatory or inhibitory, is determined at random, such that the probability of an excitatory synapse is on average 0.1. Electrical activity is evoked in the neural network by applying a 40-μA/cm2, 50-ms current in 13 randomly selected neurons at time t = 0. Simulations of 550-ms duration are performed by integrating the system of fractional-order ordinary differential equations utilizing the numerical scheme described above (Eq 15).

We perform simulations of the fractional-order neural network for different values of α, in which the network architecture and synaptic connections are identical. The collective activity of the neural network is represented by the rastergram (Fig 10A) and the the pseudo-electroencephalogram (pEEG) [28], given by Λ(t), the membrane potential averaged over all neurons (Fig 10B). We found that, for α = 1 and 0.9, network activity persists for the entire simulation. However, as α decreases, network activity self-terminates shortly following the cessation of the 50-ms applied stimulus. The time of network quiescence is not strictly a monotonic function of α, as, for example, for this specific network, activity persists longer for α = 0.6 compared with 0.8 (see also Fig. E in S1 Text).

Fig 10

Electrical activity in a fractional-order Hodgkin-Huxley neural network.

(A) Rastergram of spikes in the neural network for different values of fractional-order α. Synaptic connections and network architecture were identical in all simulations. (B) The pseudo-electroencephalogram (pEEG, Eq 25) and (C) firing rate are shown as functions of time, for the simulations in A. Firing rate is measured in a sliding 50-ms window, with 10-ms steps. (D) Interspike interval (ISI) histograms are shown for each simulation. The gray bar denotes the 50-ms applied stimulus, during which a 40-μA/cm2 current was applied to 13 randomly chosen neurons. (E) The mean network activity duration ± standard error of the mean and (F) the fraction of sustained network activity are shown in the fractional-order Hodgkin-Huxley (fHH) neural network, as a function of α (solid lines). Network simulations in which the sodium, potassium, and leak conductances, g , g , and g , respectively, are scaled such that peak current measurements are equivalent to values for particular value of α, as described in the text, are shown for comparison (dashed lines). Values in E and F are calculated for 12 network architectures.

Next, we measure the firing rate per cell as a function of time, by counting the total number of spikes in 50-ms sliding windows, with a 10-ms step (Fig 10C). We found that for α = 0.9, the firing rate is slightly higher compared with α = 1. However, for smaller α, the firing rate is initially higher during the stimulus but then sharply decreases and then becomes 0, when the network is quiescent. Histograms of the interspike intervals (ISIs) for the entire network illustrate that, for small α and networks that are quiescent following stimulus cessation, there is a higher proportion of short ISI values (Fig 10D). Interestingly, this is contrast with the relationship between spike frequency and α in the membrane patch, in which decreasing α decreases spike frequency (Fig 6C). Since the synaptic current is essentially a brief stimulus, shorter ISI values for small α are a consequence of rapid membrane depolarization, as observed in the membrane patch (Fig 2A), that triggers earlier spike peak times and thus shorter timing between excitatory synaptic current stimuli.

Measured for 12 distinct network architectures, the average duration of network electrical activity (Fig 10E, solid line) and the fraction of networks with sustained electrical activity at the end of the simulation (Fig 10F, solid line) both in general decrease with decreasing α. As in previous sections, we run network simulations in which we assume first-order V dynamics and scale ionic current conductances (here, scaling to peak current magnitudes for a constant stimulus of I = 40-μA/cm2 in the membrane patch). In this scaled conductance model, network activity is generally sustained for the entire simulation (Fig 10E and 10F, dashed lines), in agreement with fractional-order V dynamics primarily altering network activity by influencing membrane polarization and not modulating ionic currents.

Finally, we also investigate the influence of fractional-order α on synaptic integration. In Fig 11, we show the excitatory and inhibitory synaptic current, I and I , respectively, averaged over all neurons, for the network analyzed in Fig 10A–10D. By close visual inspection, there appears to be a reduction in the synaptic current magnitudes during the applied stimulus (gray bar), as α decreases. To quantify this change, we calculated the excitatory and inhibitory synaptic charge, Q and Q , the time integral of I and I , respectively, during the 50-ms stimulus, and averaged over 12 network architectures (Fig 11B, solid lines). We found that both Q and Q in general decrease as α decreases, with the the excitatory synaptic charge Q decreasing to a larger extent.

Fig 11

Synaptic activity in a fractional-order Hodgkin-Huxley neural network.

(A) The excitatory and inhibitory synaptic currents, I and I , respectively, averaged over all network neurons, are shown as a function of time, for different values of fractional-order α, for the same network as shown in Fig 10A–10D. The gray bar denotes the 50-ms applied stimulus, during which a 40-μA/cm2 current was applied to 13 randomly chosen neurons. (B) The mean of excitatory and inhibitory current charge magnitude, Q and Q , respectively, ± standard error of the mean, are shown as a function of α, calculated for 12 network architectures (solid lines). Network simulations in which the sodium, potassium, and leak conductances, g , g , and g , respectively, are scaled such that peak current measurements are equivalent to values for particular value of α, as described in the text, are shown for comparison (dashed lines).

Despite shorter ISI values and thus more overall spikes, synaptic current and charge are reduced due to the reduced spike amplitude, which reduces both the synaptic gate s (Eq 24d) and the synaptic current driving force. For small α, at the cessation of the 50-ms applied stimulus, shorter ISI values result in a larger fraction of network neurons in a refractory state, and combined with reduced excitatory synaptic charge, network activity is more likely to self-terminate. In contrast, in the scaled conductance network model, Q and Q both increase to a small extent as α decreases, and these changes negate each other, such that network activity is unaltered (Fig 11B, dashed lines).

Discussion

Summary of main findings

In this study, we investigate the properties of fractional-order neuron models and characterize how the fractional-order α influences the properties of the neuron, the spatially-distributed axon, and a random network of interconnected neurons. In the passive membrane, we show that, for small α, membrane polarization occurs rapidly, and as a consequence, the spike peak is earlier in the neuron model. Further, the bi-directional coupling between V and the ionic currents alters the ionic current time course, and the peak sodium and potassium currents increase and decrease, respectively. During a constant applied stimulus, smaller α decreases spike amplitude, and for small values of I , decreases spike frequency. Further, decreasing α reduces the range of applied stimulus amplitudes, I ∈ [I 1, I 2], that trigger spiking. In the fractional-order passive cable, we show that the impulse response is more spread out in time and space for small α, and as a consequence, sub-threshold stimuli propagate faster in the cable. In the nerve axon, spike propagation velocity also increases as α decreases, with complex dependence on spike frequency for varying amplitudes of a constant stimulus. Finally, we show in a random network of fractional-order neurons that cessation of neural network activity is more likely as α decreases. The fractional-order alters the network activity firing rate; however the dependence on α depends on the specific structure of the network.

Importantly, we found that, while the modulation of ionic currents that occur for small values of fractional-order α can account for some of the changes in spiking properties, such as spike amplitude, peak ionic current changes cannot account for many important spiking properties, including spike frequency and propagation (and their respective time-dependence during a constant stimulus) and the self-termination of neural network activity. The alteration of these spiking properties primarily arises due to fractional-order passive membrane properties, both in the membrane patch and cable, and the resulting capacitive memory effects.

Rationale for the fractional-order neuronal model

The original Hodgkin-Huxley model is natural extension of the passive membrane patch model of the parallel resistor-capacitor circuit, with the active sodium and potassium currents represented in the model by time- and voltage-dependent resistors added in parallel. As in the passive membrane patch, ideal capacitive behavior of the cell membrane and thus first-order V dynamics are assumed in the original Hodgkin-Huxley model. Motivated by the experimental evidence for non-ideal capacitive behavior in the nerve membrane, as discussed in the Introduction, in this study V dynamics are given by a differential equation of fractional order (Eq 17a).

The original Hodgkin-Huxley model assumes that four I activation, three I activation, and one I inactivation gates transition independently and that the kinetics of these gates can be represented by reversible isomerization reactions between “not-activated” (“inactivated”) and “activated” (“not-inactivated”) states (Eqs 17b–17d). Although we know that these assumptions are not strictly correct, e.g., I activation and inactivation is not independent and gating transitions often involve multiple “not-activated” states [29], these complications can be easily addressed by incorporating more detailed gating mechanisms, as many studies have done. However, an assumption that is more challenging to address, implicit in this formulation, is that channel gating is a Markov process, i.e., gating transitions only depend on the current system state and do not depend on prior history. While the Markov assumption is generally accepted in the field of ion channel biophysics, some studies have shown that channel gating kinetics present memory in the form of short- and long-term correlations between open- and closed-dwell times [30-33], and that fractal gating kinetics, i.e., transition rates that scale inversely with the duration in a given state, can reproduce such correlations [30, 34]. Mathematical treatments have begun to address the relationship between these dwell time correlations and fractional-order dynamics in theoretical channel models [35, 36] and more generally the relationship between fractional derivatives and fractal dimension [37]. However, to our knowledge, no studies have investigated the relationship between fractional-order gating dynamics and dwell times in a biophysically-detailed gating model. Further work is needed to investigate this relationship. In this study, to specifically address the influence of non-ideal membrane capacitive behavior and fractional-order V dynamics, we assume Markovian channel gating, such that gating variable dynamics remain unchanged from the original formulation and assumed first-order.

Comparison with prior work

To our knowledge, only two prior studies have investigated the fractional-order Hodgkin-Huxley model [38, 39], and both also showed that the fractional-order influenced the time course of the membrane potential V and gating variables, including an earlier spike peak time, as we show in Fig 2. Both of these studies also assumed that the gating variables are governed by fractional-order dynamics. As discussed above, ion channel gating memory may be present. However, these studies assume the same fractional-order for V and all three gating variables, which we believe is an unjustified modification of the original model given the distinctly different biophysical origin for each source of memory. Further, our analysis is significantly more complete. Here, we characterize the influence of fractional-order on many properties, including spike frequency and amplitude, as well properties of spike propagation in the nerve axon and neural network activity, which are not considered in the aforementioned studies.

A few prior studies have investigated fractional-order V dynamics in other excitable cell models. Shi and Wang demonstrated that the fractional-order Morris-Lecar neuron model can exhibit a wider range of bursting behavior than can be reproduced by the original, first-order model [40]. Jun et al. similarly show that fractional V dynamics in the Hindmarsh-Rose neuronal model alters spiking patterns and also found a larger applied current threshold for repetitive spiking for smaller fractional-order [41], in agreement with our findings in the Hodgkin-Huxley model (Fig 6B and 6C). Teka et al. demonstrated that fractional dynamics and memory in the leaky integrate-and-fire model altered spike adaptation and firing patterns, including spike latency and interspike variability [22].

Momani et al. studied the coupling of fractional-order neurons modeled using the minimal FitzHugh-Nagumo [42]. However, this study was primarily focused on analysis of numerical techniques. Moaddy et al. studied a neuronal model represented by the fractional passive membrane (analyzed in Fig 1) with the addition of a fractional-order inductor [43]. The fractional-order inductor essentially augments the passive membrane model with an additional slow variable, which is comparable to the FitzHugh-Nagumo model. This study was also primarily focused on numerical techniques; however the authors do find that synchronization of coupled (2 or 3) neurons occurs more quickly as the fractional-order decreased, which appears in agreement with our analysis of larger networks (Fig 10), although the number of networks simulated in our study is too small for a definitive conclusion on synchrony.

While this study and the aforementioned modeling studies focus on fractional-order dynamics for the membrane potential and the influence on neuronal spiking properties, fractional-order dynamics for spiking activity has also been observed experimentally. Lundstrom et al. showed that in response to a sinusoidal stimulus, the firing rate of neocortical pyramidal neurons could be well-approximated by a fractional-order derivative of the input [44]. Anastasio showed that the frequency response of the firing rate, relative to eye position, of premotor neurons in the oculomotor system exhibited fractional-order dynamics, with a phase lag less than 90 degrees (see Fig 1B), suggesting fractional-order dynamics are involved in processing of saccade-related activity [45].

Physiological implications

The biophysical basis underlying short- and long-term memory is an important area of neuroscience research [46, 47]. The molecular, cellular, and network level mechanisms underlying memory are complex and inevitably involve the interactions between dynamics across a wide range of spatiotemporal scales. An important consequence of fractional-order dynamics is that the membrane potential V has memory, in the sense that V depends on its entire history from an initial state, albeit with decremental weighting (Fig. A in S1 Text). For fractional-orders in the range considered here, V dynamics theoretically may be influenced by history on the order of minutes to hours (see Fig. A in S1 Text). In this study, the small integration time step required for numerical stability (Fig. B in S1 Text) prohibited parameter studies of minute- to hour-long simulations; thus, although our analysis suggests that several important spiking properties, such as spike frequency and propagation velocity, approach asymptotic values on the order of tens of milliseconds, it is possible that dynamics influenced by memory emerge at this longer time scale. Further, in some neural networks, firing rates approach asymptotic values, while in other networks, firing rate is quite variable in time (Fig. E in S1 Text). Thus, the influence of memory on the longer time scale of minutes to hours may be significant in specific network architectures. While the mechanisms underlying neuronal short-term memory are necessarily more complex than memory associated with the cellular membrane dielectric behavior, our study shows that this memory behavior contributes to, or at least modulates, the electrical activity underlying short-term memory. Further work is necessary to thoroughly address this question.

Practical considerations for neuronal modeling

In this study, we demonstrate that accounting for fractional-order dynamics of the membrane potential alters many properties at the level of the neuron, axon, and network. While we have shown that accounting for fractional-order dynamics may be important, it is not our intention to argue that studies using the “standard” excitable cell model, assuming ideal (first-order) capacitive behavior, are inaccurate. In the neuronal membrane patch following a brief stimulus, the fractional-order α primarily influences measures quantitatively, e.g., peak current magnitudes and refractoriness, but did not alter the overall qualitative behavior of the neuron. These properties may be similarly reproduced by altering model parameters.

However, at longer time scales, on the order of tens to hundreds of milliseconds, capacitive memory and fractional-order membrane properties are more significant, such that scaling ionic current conductances does not reproduce important spiking properties, such as spike frequency. Further, at the level of the nerve axon and neural network, we found that fractional-order V dynamics can significantly alter emergent properties of the system that are not reproduced by appropriately rescaling ionic current conductances. For example, modulating peak ionic currents does not reproduce the significant increase in propagation velocity in the fractional-order nerve axon model (Fig 8), nor the self-termination of network electrical activity in a random network (Fig 10).

Thus, it is a difficult challenge to determine for which models and settings it is necessary to represent V with fractional-order dynamics. Our study suggests that it is important particularly at longer time scales and when accounting for neural membranes that are interconnected via propagation in a nerve axon or synaptic connections in a random network. However, this conclusion may be model-specific (see Limitations section). The use of fractional differential equations to model excitable cells and biological systems in general is presently limited by the complex mathematics involved and limitations on numerical methods. While numerical integration of fractional differential equations is typically more computational expensive compared with integer-order differential equations, advanced integration methods have been recently developed [48-50]. As these methods become more widely used, characterized, and optimized, we expect that modeling of fractional-order excitable cell systems and networks will become more computational feasible, which will allow for further determination of models in which accounting for fractional-order dynamics is important.

Limitations

The prior studies of fractional-order dynamics, mentioned above, typically focused on minimal neuron models, such as the FitzHugh-Nagumo and Morris-Lecar models, which while valuable, are limited in their biophysical detail. The Hodgkin-Huxley model of the squid giant axon neuron is a classical biophysical model of an excitable cell that is well-characterized, and thus it was a reasonable model to use as a starting point to investigate the influence of fractional-order V dynamics on neuronal properties. However, more biophysically-detailed neuronal models relevant to mammalian physiology have been described, incorporating more detailed sodium, potassium, and calcium currents [51, 52]. It may also be important to account for intracellular calcium signaling, specifically stochastic dynamics that can influence electrical activity [53, 54], including spiking frequencies and synaptic transmission.

The properties of the neural networks considered in this study were random; while this type of network provides general insight, the dynamics of activity in networks with a physiological architecture may differ, particularly in the context of disease states such as epileptic seizures [55]. Further work is needed investigate fractional-order dynamics in these more physiologically-realistic neuronal models and networks. Nonetheless, our study provides valuable insight and serves as a reference for investigating fractional-order dynamics in models of neuronal activity.

Supporting Analysis, Figures, and Tables.

Supporting information includes model equations, initial conditions, and parameters, analytical/numerical solutions of passive membrane patch and cable, memory weighting terms, numerical stability analysis, gating variables during spiking, spike frequency and propagation velocity analysis, and neural network firing rate analysis.

(PDF)

Click here for additional data file.

MATLAB membrane patch, nerve axon, and neural network simulation code.

(ZIP)

Click here for additional data file.