Power spectral analysis of voltage-gated channels in neurons

This article develops a fundamental insight into the behavior of neuronal membranes, focusing on their responses to stimuli measured with power spectra in the frequency domain. It explores the use of linear and nonlinear (quadratic sinusoidal analysis) approaches to characterize neuronal function. It further delves into the random theory of internal noise of biological neurons and the use of stochastic Markov models to investigate these fluctuations. The text also discusses the origin of conductance noise and compares different power spectra for interpreting this noise. Importantly, it introduces a novel sequential chemical state model, named p₂, which is more general than the Hodgkin–Huxley formulation, so that the probability for an ion channel to be open does not imply exponentiation. In particular, it is demonstrated that the p₂ (without exponentiation) and n⁴ (with exponentiation) models can produce similar neuronal responses. A striking relationship is also shown between fluctuation and quadratic power spectra, suggesting that voltage-dependent random mechanisms can have a significant impact on deterministic nonlinear responses, themselves known to have a crucial role in the generation of action potentials in biological neural networks.

1 Introduction

The purpose of this article is to elucidate the fundamental behavior of excitable neuronal membranes by using recent methods in the frequency domain. Since the historical work of the French mathematician and physicist Fourier (1822) in The Analytical Theory of Heat, it is known that many kinds of signals admit a dual representation either as a real valued function u(t) of the time variable t or as a complex valued function û(ω) of the frequency variable ω where û is the Fourier transform of u. Using this approach, individual excitable cells usually studied by their responses to constant stimuli can also be investigated by their responses to multi-sinusoidal stimuli over a broad frequency range.

With the current clamp technique, the reactions of individual cells to current stimuli generally show threshold impulses of the membrane potential that provide the means for signal transmission throughout the nervous system. Precisely because of this threshold property, it is difficult to determine the ionic currents underlying those action potential impulses. A major technical advance in measuring neuronal properties occurred with the advent of the voltage clamp technique (Bear et al., 2016), which was first invented by Cole and Marmont (1949) and later exploited by Hodgkin and Huxley (1952) who developed the Hodgkin–Huxley (HH) equations that have become the gold standard for most neuronal models in real time simulations. With the voltage clamp, a retroactive electronic device controls the membrane potential such that the neuronal properties can be quantitatively determined from the measured current elicited by a change in the potential.

Both current clamp and voltage clamp techniques can be used in the frequency domain to characterize neuronal function. A typical approach considers linear analysis by calculating impedance and admittance, as described by Mauro (1970) for the squid giant axon. However, neuronal behavior is fundamentally nonlinear due to the voltage dependence of most ionic channels (for instance potassium or sodium). Quadratic sinusoidal analysis (QSA) is a recent method developed by Magnani and Moore (2011) that provides a fundamental insight of the linear and quadratic neuronal behaviors using matrix calculus in the frequency domain. Concrete applications have been done with neurons involved in the oculo-vestibular integrator (Magnani et al., 2013) as well as with neurons of the medial entorhinal cortex which are part of the grid cell network (Magnani et al., 2014). These linear and nonlinear approaches in the frequency domain are much more efficient and concise than time domain methods for extracting stationary and dynamic features from neurons by slightly perturbing them with multi-sinusoidal signals around a steady state.

Although such a smooth deterministic description by frequency waves is able to capture fundamental properties of the neuronal function, biological neurons are significantly perturbed by internal noise. Among the different kinds of noise sources described by Stevens (1972), conductance fluctuations reflect the stochastic nature of ionic channels at the microscopic level. For this reason, stochastic Markov models have been used in this article to investigate the intrinsic fluctuations and their relationships to the complicated nonlinear behavior of neurons. This approach extends QSA analysis with stochastic Markov simulations and compares the power spectra of linear, quadratic and stochastic neuronal processes.

Some of the earliest measurements of membrane voltage noise were done on the node of Ranvier by Verveen and Derksen (1968) who suggested that the opening and closing of ionic channels lead to voltage fluctuations. Later, voltage clamp measurements on squid axons suggested that current channel noise was filtered by the axonal membrane that can be measured as a voltage power spectrum. Extensive measurements have been made of spontaneous fluctuations under both current and voltage clamp conditions along with linear impedance analysis to assess the basis of the measured noise power spectrum (Fishman et al., 1983; Conti and Wanke, 1975; Conti, 1984; Poussart, 1969; Poussart and Ganguly, 1977).

An early and perceptive paper by Stevens (1972) showed that the Hodgkin–Huxley equations themselves provide two different interpretations of the origin of conductance noise, namely the opening and closing of whole conductance channel whose opening is controlled by multiple gating particles (n, m or h), or alternatively, fluctuations of individual gating particles. The first case involves probability and correlation functions having multiple terms. For the potassium conductance where $g_{\text{K}}=g_{\text{K}}^{\text{max}}{n}^4$, this leads to conductance noise power spectra consisting of four Lorentzian terms. In the second case, the Hodgkin–Huxley equations are linearized and the potassium conductance power spectrum consists of a single term related to the potassium channel time constant. In the first case, the nonlinear properties of the channel (n⁴) would be involved in the origin of the spontaneous fluctuation, while in the second case, the linearized impedance (n) would likely be a good predictor of the voltage noise.

The measured squid axon spontaneous noise and the nonlinear power spectra of the Hodgkin–Huxley model both show resonance in the voltage measurements and non-resonating Lorentzian functions under voltage clamp. There is strong experimental and theoretical evidence that the spontaneous conductance noise cannot be predicted by the linear response from the same axon (Fishman et al., 1983; Poussart, 1969).

With the advent of measurements of single ionic channels in neuronal membranes, many of the detailed properties of different ionic channels as well as their macromolecular basis have dominated excitable membrane research. These findings have stimulated the development of stochastic ionic conductance models, again using the Hodgkin–Huxley equations as a fundamental basis. The gold standard method to simulate the stochastic behavior appears to be Markov models which provide fluctuation noise power spectra identical to the first Hodgkin–Huxley interpretation by Stevens discussed above. Numerous papers have derived stochastic fluctuation power spectra based on this nonlinear character (O'Donnell and van Rossum, 2014; Goldwyn and Shea-Brown, 2011, for instance). In addition to squid axon, measurements from nodes of Ranvier (Conti and Wanke, 1975; Conti, 1984; Sigworth, 1980; Elinder et al., 2001) and other preparations are consistent with the nonlinear origin of ionic power spectra. The interpretation, namely that certain nonlinear properties can show a probabilistic or stochastic character, suggests that they are involved in the fundamental origins of spontaneous channel noise in neurons. Since this clearly indicates that the stochastic behavior of excitable membranes is nonlinear, it is useful to quantitatively compare the power spectrum content of the nonlinear QSA responses with the corresponding simulations of Markov models. This will be done rigorously from the Hodgkin–Huxley equations described in the methods.

In the methods section, the Hodgkin–Huxley model is briefly reviewed. Linear analysis in the frequency domain is introduced based on the work of Mauro (1970). The quadratic sinusoidal analysis (QSA) is introduced in two ways, first the basic theory as described by Magnani and Moore (2011) and second, a new algorithm of frequency averaging to calculate power spectra. Fluctuation simulations are done with Markov models based on the work of Goldwyn et al. (2011) and Goldwyn and Shea-Brown (2011). The theoretical expressions for fluctuation power spectra are derived from the work of O'Donnell and van Rossum (2014).

In the results section, the exponentiation of the potassium n⁴ model is reduced to the minimal degree of nonlinearity n². Then, a novel sequential chemical state model named p₂ is introduced as a generalization of the n² model without exponentiation. Linear analysis, fluctuation simulations and theoretical power spectra are applied to the p₂ model by adapting the previous methods. Linear and quadratic functions are compared to neuronal fluctuations by adapting the previous methods. The p₂ and n⁴ models are compared and it is demonstrated that they can produce similar neuronal responses. Remarkably, it is illustrated how the fluctuation-dissipation theorem can be violated at depolarized membrane potentials. The time constants of the p₂ and n⁴ models seem more consistent with stochastic analysis than linear analysis.

Finally, the discussion section deals with the origin of fluctuations in the nonlinear neuronal responses. Surprisingly, simulations show that, for certain stimulus amplitudes and membrane surfaces, random voltage-dependent fluctuations can significantly modify deterministic nonlinear responses.

2 Methods

2.1 Hodgkin–Huxley model

The standard Hodgkin–Huxley model was originally proposed by Hodgkin and Huxley (1952) to describe the initiation and propagation of action potentials in the squid giant axon based on voltage clamp experiments. In this model, the current across the lipid bilayer is defined by

$C_m\frac{dV}{dt}=I-I_{\text{L}}-I_{\text{K}}-I_{\text{Na}}$

where C_m is the membrane capacitance, V is the membrane potential, I is the total membrane current, I_L is the leak current, I_K is the current through potassium ion channels and I_Na is the current through sodium ion channels.

The ionic currents are expressed with conductances

$\{\begin{array}{l}{I}_{\text{L}}=g_{\text{L}}(V-V_{\text{L}})\\ I_{\text{K}}=g_{\text{K}}{n}^4(V-V_{\text{K}})\\ I_{\text{Na}}=g_{\text{Na}}{m}^{3}h(V-V_{\text{Na}})\end{array}$

where g_L is the leak conductance, g_K is the maximum potassium conductance and g_Na is the maximum sodium conductance. The constants V_L, V_K and V_Na are reversal potentials for leak, potassium ion channel and sodium ion channel respectively. The gating variables n, m and h represent potassium channel activation, sodium channel activation and sodium channel inactivation respectively, their values are constrained between 0 and 1.

The kinetics of gating variables satisfies the following first order differential equations determined by pairs of rate constants α_i(V) and β_i(V)

$\{\begin{array}{l}\frac{dn}{dt}={\alpha }_n(V)(1-n)-{\beta }_n(V)n\\ \frac{dm}{dt}={\alpha }_m(V)(1-m)-{\beta }_m(V)m\\ \frac{dh}{dt}={\alpha }_h(V)(1-h)-{\beta }_h(V)h\end{array}$

At the steady state corresponding to a membrane potential level V₀, the time derivatives vanish with constants

$\{\begin{array}{l}{n}_{\infty }(V_0)={\alpha }_n(V_0){\tau }_n(V_0)\\ m_{\infty }(V_0)={\alpha }_m(V_0){\tau }_m(V_0)\\ h_{\infty }(V_0)={\alpha }_h(V_0){\tau }_h(V_0)\end{array}$

where

${\tau }_n={({\alpha }_n+{\beta }_n)}^{-1},{\tau }_m={({\alpha }_m+{\beta }_m)}^{-1},{\tau }_h={({\alpha }_h+{\beta }_h)}^{-1}$

2.2 Linear analysis

Cole (1941) was among the first to use the frequency domain to investigate neurons and suggested that the linear response of the axon membrane could be modeled with equivalent circuit elements such as inductances, capacitances and resistances. Hodgkin and Huxley also showed that their findings were consistent with the potassium conductance being described as an inductive reactance. Later, a mathematical equivalence between nonlinear conductance Hodgkin–Huxley models and electrical circuits was developed by Mauro (1970) and experimentally confirmed on the squid axon.

In papers published by Mauro (1970) and Fishman et al. (1977), the linearization of the Hodgkin–Huxley equations in the frequency domain is obtained for small perturbations at steady state. More precisely, let X = (V, n, m, h) be a dynamic state, X₀ = (V₀, n₀, m₀, h₀) a steady state and $\frac{dX}{dt}=F(X)$ the system of differential equations. Clearly, n₀ = n_∞(V₀), m₀ = m_∞(V₀) and h₀ = h_∞(V₀) where V₀ can be controlled directly in voltage clamp or indirectly in current clamp. Then, differential calculus linearizes the system for a small perturbation steady state

$\delta F=F(X)-F(X_0)~\delta X^{T}DF(X_0)$

It follows linearization of the kinetic equations, for example for the potassium variable

$\delta \frac{dn}{dt}=\frac{d{\alpha }_n}{dV}\delta V-{\alpha }_n\delta n-n_0\frac{d{\alpha }_n}{dV}\delta V-{\beta }_n\delta n-n_0\frac{d{\beta }_n}{dV}\delta V$

This leads to the first order linear differential equation

$\frac{d\delta n}{dt}+({\alpha }_n+{\beta }_n)\delta n=\left[\frac{d{\alpha }_n}{dV}-n_0\frac{d({\alpha }_n+{\beta }_n)}{dV}\right]\delta V$

Applying the Fourier transform

$i\omega \widehat{\delta n}+({\alpha }_n+{\beta }_n)\widehat{\delta n}=\left[\frac{d{\alpha }_n}{dV}-n_0\frac{d({\alpha }_n+{\beta }_n)}{dV}\right]\widehat{\delta V}$

This can be solved in the frequency domain as

$\widehat{\delta n}=\frac{\frac{d{\alpha }_n}{dV}-n_0\frac{d({\alpha }_n+{\beta }_n)}{dV}}{i\omega +{\alpha }_n+{\beta }_n}\widehat{\delta V}$

and similarly for $\widehat{\delta m}$ and $\widehat{\delta h}$. In order to simplify notations, this fraction will be denoted as $\widehat{D_n}$, $\widehat{D_m}$ and $\widehat{D_h}$ respectively.

Similarly, linearization of the ionic currents is given by

$\{\begin{array}{l}\delta I_{\text{L}}=g_{\text{L}}\delta V\\ \delta I_{\text{K}}=4g_{\text{K}}{n}_0^3(V_0-V_{\text{K}})\delta n+g_{\text{K}}{n}_0^4\delta V\\ \delta I_{\text{Na}}=3g_{\text{Na}}{m}_0^2{h}_0(V_0-V_{\text{Na}})\delta m+g_{\text{Na}}{m}_0^3(V_0-V_{\text{Na}})\delta h\\ +g_{\text{Na}}{m}_0^3{h}_0\delta V\end{array}$

which can be solved in the frequency domain as

$\{\begin{array}{l}\widehat{\delta I_{\text{L}}}=g_{\text{L}}\widehat{\delta V}\\ \widehat{\delta I_{\text{K}}}=g_{\text{K}}\left[4n_0^3(V_0-V_{\text{K}})\widehat{D_n}+n_0^4\right]\widehat{\delta V}\\ \widehat{\delta I_{\text{Na}}}=g_{\text{Na}}[3m_0^2{h}_0(V_0-V_{\text{Na}})\widehat{D_m}+m_0^3(V_0-V_{\text{Na}})\widehat{D_h}\\ +m_0^3{h}_0]\widehat{\delta V}\end{array}$

Similarly, linearization of the current across the lipid bilayer is given by

$C_m\frac{d\delta V}{dt}=\delta I-\delta I_{\text{L}}-\delta I_{\text{K}}-\delta I_{\text{Na}}$

or equivalently in the frequency domain

$C_{m}i\omega \widehat{\delta V}=\widehat{\delta I}-\widehat{\delta I_{\text{L}}}-\widehat{\delta I_{\text{K}}}-\widehat{\delta I_{\text{Na}}}$

The linearized response current $\widehat{\delta I}$ to the linearized stimulus voltage $\widehat{\delta V}$ is characterized by the admittance

$\hat{Y}=\frac{\widehat{\delta I}}{\widehat{\delta V}}$

$\begin{array}{l}\hat{Y}=C_{m}i\omega \\ +g_{\text{L}}\\ +g_{\text{K}}\left[4n_0^3(V_0-V_{\text{K}})\widehat{D_n}+n_0^4\right]\\ +g_{\text{Na}}\left[3m_0^2{h}_0(V_0-V_{\text{Na}})\widehat{D_m}+m_0^3(V_0-V_{\text{Na}})\widehat{D_h}+m_0^3{h}_0\right]& (1)\end{array}$

Equivalently, the impedance is defined as

$\begin{array}{ll}\hat{Z}=\frac{1}{\hat{Y}}& (2)\end{array}$

Note that the frequency f in Hz is interchangeable with angular frequency ω = 2πf, both are used in this article depending on the context.

One of the most striking features of neurons and their models is the presence of an impedance resonance at certain membrane potentials V₀. This is illustrated for the voltage clamped Hodgkin–Huxley model especially by the two amplitude peaks at 5 and 25.2 mV depolarizations shown in Figure 1A. The plots are superimposed magnitudes |Ẑ(f, V₀)| of the Equation 2, which is identical to a small signal sinusoidal stimulus of the full nonlinear Hodgkin–Huxley equations. The simulation results are similar to actual measurements on squid axons independent of the electrode properties.

Figure 1

Figure 1. Impedance functions for the Hodgkin–Huxley model. Different membrane potential displacements V₀ were applied from the resting potential of zero. Curves from top to bottom at lowest frequencies: –5, 0, 5, 10.2, 15.2, 25.2 in mV. Abscissa: frequency f in kHz (logarithmic scale). Ordinate: magnitude |Ẑ(f, V₀)| in Ω (logarithmic scale). Parameters of the simulation are given in Table 1. (A) Simulations with potassium and sodium conductances. (B) Simulations with potassium conductance only (g_Na = 0).

Although impedance resonance is a linear property of the Hodgkin–Huxley equations, it is due to the voltage dependence of the steady state values of the ionic conductances. In particular, $\frac{d{n}_{\infty }}{dv}>0$ for the potassium conductance [i.e. n_∞(V₀) does increase with depolarizing levels V₀]. In Equation 1, the potassium term in g_K was shown to be equivalent to an inductive reactance by Mauro (1970). Similarly the sodium term g_Na can be described by other circuit elements. Thus, the Hodgkin–Huxley model or any excitable cell can be analyzed by a piecewise linear analysis at different membrane potentials as shown in Figure 1A. Data collected in this manner over a range of membrane potentials allows one to determine the voltage dependence of the active conductances in addition to passive properties, and in turn construct a system of nonlinear differential equations for a particular model. Thus, a further advantage of frequency domain measurements is that model discrimination using parameter estimation can be more accurate if both real time and impedance results are used, as shown by Murphey et al. (1995).

Since neurons are composed of a minimum of two conductances, inward sodium or calcium currents and various outward potassium and other currents, it is useful to consider their individual contributions to the frequency domain behavior. This paper is focused on the potassium conductance as a model of any of the individual conductances, thus g_Na = 0. Figure 1B illustrates one broad impedance resonance maximum for the Hodgkin–Huxley potassium conductance alone, which shifts to higher frequencies with increasing depolarizations. In contrast, the impedances of Figure 1A have sharper resonances and two peaks (5 and 25.2 mV) clearly due to the presence of the sodium conductance in conjunction with potassium. Thus, one useful aspect of this analysis is that the number of resonance peaks can give an indication of the minimum number of active conductance processes present.

2.3 Quadratric sinusoidal analysis

Biological neurons and their models are fundamentally nonlinear systems that sometimes significantly contradict the linear superposition principle. System identification methods such as Volterra and Wiener series have been widely used to characterize nonlinear neuronal functions. In the pioneering work of Marmarelis and Naka (1972), Wiener theory was applied to predict the nonlinear behavior of a neuron chain in the catfish retina. Unfortunately, it is practically difficult and time-consuming to calculate the kernel coefficients of the Volterra and Wiener series, as these methods generally require extensive data analysis, analogous to averaging a small signal from extraneous noise.

In general, the use of random broad band stimuli similar to a white noise is time consuming, even for linear methods requiring averaging over experiments. However, significant reduction of experimental time is possible if a pseudo random stimulus containing a limited number of frequencies is applied to the preparation (Fishman et al., 1977; Poussart, 1969). Responses to such stimuli have an excellent signal to noise ratio and little or no averaging is necessary when measuring the impedance of a single cell.

In the case of nonlinear systems, multi-sinusoidal stimuli can be used to precisely measure deterministic linear and quadratic responses, provided that generated output frequencies do not overlap at first and second orders. For example, input frequencies f₁ = 1, f₂ = 2, f₃ = 3, f₄ = 4 (in Hz) generate unwanted overlaps because f₁ + f₂ = f₃ or f₁ + f₄ = f₂ + f₃. Multi-sinusoidal stimuli with selected frequencies to avoid overlap at first and second orders were used by Magnani and Moore (2011) to study subthreshold neuronal responses, as well as previously by Victor and Shapley (1980) to study cat retinal ganglion cells.

Such non-overlapping multi-sinusoidal stimuli must have sufficiently large amplitudes to elicit both linear and quadratic responses, while remaining sufficiently small to avoid higher order contamination. The quadratic sinusoidal analysis, termed QSA, was introduced by Magnani and Moore (2011) to provide a flexible way to capture the linear and quadratic neuronal functions at subthreshold membrane potentials, as well as to compare biological experiments with theoretical models (Magnani et al., 2013, 2014). In particular, the subthreshold membrane potential just below the threshold is fundamentally a nonlinear process which is critically involved in action potential generation.

It is well known that a sinusoidal signal can be expressed as a sum of complex exponentials in Fourier analysis, for example :

$cos(\omega t)=\frac{1}{2}{e}^{i\omega t}+\frac{1}{2}{e}^{-i\omega t}$

where i² = −1 and ω is the angular frequency. More generally, a multi-sinusoidal signal of N frequencies can be expressed as a sum of complex exponentials

$x(t)={\displaystyle \sum _{k\in \text{Γ}}}{x}_k{e}^{i{\omega }_{k}t}$

where Γ = {−N, ⋯ , −1, +1, ⋯ , +N} enumerates integers between −N and +N (zero excluded), x_k are complex Fourier coefficients such that $x_{-k}=\overline{x_k}$ and ω_k are angular frequencies such that ω_−k = −ω_k. The duration T (s) of the experiment determines the lowest frequency $\frac{1}{T}$ (Hertz). The angular frequencies are defined by integer multiples of the lowest frequency, that is to say ${\omega }_k=2\pi \frac{n_k}{T}$ where n_k are integers such that n_−k = −n_k. Each Fourier coefficient is a complex number that can be decomposed as $x_k=\left|x_k\right|e^{i{\theta }_k}$ where |x_k| is the amplitude (nonnegative real number) and θ_k is the phase (between 0 and 2π). It is good practice to randomize the phases of multi-sinusoidal stimuli to avoid biases in neuronal responses.

Such a multi-sinusoidal stimulus x(t) can be applied to neuronal cells in voltage clamp and current clamp experiments. Assuming that quality criteria for the QSA method are satisfied (Magnani and Moore, 2011), the output signal measured by the electrode can be decomposed as

$y(t)=y_0+y_1(t)+y_2(t)$

where y₀ is the DC component, y₁(t) is the linear component and y₂(t) is the quadratic component. More specifically

$\begin{array}{ll}y(t)=y_0+{\displaystyle \sum _{k\in \text{Γ}}}{l}_k{x}_k{e}^{i{\omega }_{k}t}+{\displaystyle \sum _{i,j\in \text{Γ}}}{b}_{i,j}{x}_i{e}^{i{\omega }_{i}t}{x}_j{e}^{i{\omega }_{j}t}& (3)\end{array}$

where l_k and b_{i, j} are complex numbers characterizing the neuronal response to the stimulus, by similarity with impedance or admittance. By convention, b_{−k, k} = 0 for all k ∈ Γ so that all DC components are encoded in y₀.

The QSA theory is based on a vector representation of multi-sinusoidal signals in an orthonormal vector basis {_{ek}k ∈ Γ} representing the complex exponentials ${\left\{e^{i{\omega }_{k}t}\right\}}_{k\in \Gamma }$. In this way, the multi-sinusoidal stimulus x(t) is encoded as a time independent vector

$\text{x}={\displaystyle \sum _{k\in \text{Γ}}}{x}_k{\text{e}}_k$

The corresponding time dependent vectors are defined by

${\text{x}}_t={\displaystyle \sum _{k\in \text{Γ}}}{x}_k{e}^{i{\omega }_{k}t}{\text{e}}_k$

Putting the coefficients l_k in a row matrix L and the coefficients b_{i, j} in a square matrix B, Equation 3 can be reformulated with linear algebra as

$\begin{array}{ll}y(t)=y_0+\text{L}{\text{x}}_t+{\text{x}}_t^T\text{B}{\text{x}}_t& (4)\end{array}$

By noticing that

$\begin{array}{ll}{\displaystyle \sum _{i,j\in \text{Γ}}}{b}_{i,j}{x}_i{e}^{i{\omega }_{i}t}{x}_j{e}^{i{\omega }_{j}t}& ={\displaystyle \sum _{i,j\in \text{Γ}}}{b}_{-i,j}\overline{x_i{e}^{i{\omega }_{i}t}}{x}_j{e}^{i{\omega }_{j}t}\end{array}$

we are led to define the QSA matrix Q by

$Q_{i,j}=B_{-i,j}$

so that

$y(t)=y_0+\text{L}{\text{x}}_t+{\overline{\text{x}}}_t^T\text{Q}{\text{x}}_t$

Remarkably, the QSA matrix Q is Hermitian (Magnani and Moore, 2011), which means that

${\overline{\text{Q}}}^T=\text{Q}$

This algebraic approach has been widely used in modern physics, especially in quantum physics where Hermitian operators represent physical observables. In particular, the eigenvalues of the QSA matrix are real numbers which form a spectrum that can be interpreted as a signature of the quadratic neuronal function.

The most glaring property of the quadratic neuronal response is the generation of second order frequencies that are not present in the stimulus. More precisely, Equation 3 shows that the linear components $l_k{x}_k{e}^{i{\omega }_{k}t}$ contain only the stimulus frequencies ω_k, while the quadratic components $b_{i,j}{x}_i{e}^{i{\omega }_{i}t}{x}_j{e}^{i{\omega }_{j}t}=b_{i,j}{x}_i{x}_j{e}^{i({\omega }_i+{\omega }_j)t}$ contain new frequencies ω_i + ω_j that are not present in the stimulus in the absence of overlap. The row matrix L and the QSA matrix Q can be interpreted as linear and quadratic filters respectively, by extension of the impedance or admittance concept.

2.4 Multi-sinusoidal power spectra

The concept of power spectrum is useful in statistical analysis, both for signal processing and for stochastic processes, so it is natural to use it in this paper. The power spectrum of a multi-sinusoidal signal u(t) is given by |û(ω)|² where û is the Fourier transform of u and ω is the frequency variable. The power spectrum can be averaged over a set of measurements.

For a single measurement of the output y(t) in response to a stimulus x(t) with fundamental frequencies |ω_k|, the multi-sinusoidal power spectra cover the positive frequencies of linear and quadratic analyses:

• Linear power spectrum at fundamental frequencies:

$S_L(|{\omega }_k|)=|{\hat{y}}_1(|{\omega }_k|){|}^2$

• Quadratic power spectrum at frequency doubling:

$S_D(2|{\omega }_k|)=|{\hat{y}}_2(2|{\omega }_k|){|}^2$

• Quadratic power spectrum at frequency sums for |ω_i| ≠ |ω_j|:

$S_P(|{\omega }_i|+|{\omega }_j|)=|{\hat{y}}_2(|{\omega }_i|+|{\omega }_j|){|}^2$

• Quadratic power spectrum at frequency differences for |ω_i| ≠ |ω_j|:

$S_M(||{\omega }_i|-|{\omega }_j||)=|{\hat{y}}_2(||{\omega }_i|-|{\omega }_j||){|}^2$

The quadratic power spectra are indexed by second order frequencies, which are not stimulus frequencies because they were chosen without overlap. Thus, it would be convenient to have also a quadratic power spectrum indexed by fundamental frequencies as an alternative representation. To this end, a function similar to the “R summation function” of Magnani et al. (2013) is introduced in a different way below, computing the mean squared quadratic output by matrix columns:

$S_R({\omega }_j)=\frac{1}{2N}{\displaystyle \sum _{i\in \text{Γ}}}|Q_{i,j}\overline{x_i}{x}_j{|}^2$

Although these multi-sinusoidal power spectra reflect exact neuronal responses, they are not accurate because they are defined over a small set of non-overlapping frequencies. Therefore, it is necessary to average multiple measurements to increase the accuracy of the power spectra.

To perform this averaging, a set of M stimuli x^(m)(t) is generated with sets of non-overlapping random frequencies Ω^(m) for m = 1, …, M. Each set of first order frequencies Ω^(m) determines a set of second order frequencies Ξ^(m). The global sets of first and second order frequencies are defined by merging the individual sets, respectively

$\{\begin{array}{l}\Omega ={\displaystyle \bigcup _{m=1}^M}{\Omega }^{(m)}\\ \Xi ={\displaystyle \bigcup _{m=1}^M}{\Xi }^{(m)}\end{array}$

Importantly, different individual sets of the same order can share frequencies, so it is necessary to count redundancies

$\{\begin{array}{l}{N}_{\Omega }(\omega )=\#\left\{m\mid \omega \in {\Omega }^{(m)}\right\}\\ N_{\Xi }(\xi )=\#\left\{m\mid \xi \in {\Xi }^{(m)}\right\}\end{array}$

where # indicates the number of elements in a set.

Fourier analysis of the measured linear responses $y_1^{(m)}(t)$ yields output sets ${\hat{y}}_1^{(m)}(\omega )$ defined at the first order frequencies ω ∈ Ω^(m) and zero elsewhere. This allows to calculate the averaged linear power spectrum for ω ∈ Ω

$S_L(\omega )=\frac{1}{N_{\Omega }(\omega )}{\displaystyle \sum _{m=1}^M}|{\hat{y}}_1^{(m)}(\omega ){|}^2$

where ω ∈ Ω implies ${\mathcal{N}}_{\Omega }(\omega )\ge 1$ so the denominator is not zero.

Fourier analysis of the measured quadratic responses $y_2^{(m)}(t)$ yields output sets ${\hat{y}}_2^{(m)}(\xi )$ defined at the second order frequencies ξ ∈ Ξ^(m) and zero elsewhere. This allows to calculate the averaged quadratic power spectrum for ξ ∈ Ξ

$S_2(\xi )=\frac{1}{N_{\Xi }(\xi )}{\displaystyle \sum _{m=1}^M}|{\hat{y}}_2^{(m)}(\xi ){|}^2$

where ξ ∈ Ξ implies $\mathcal{N}(\xi )\ge 1$ so the denominator is not zero. The partial quadratic power spectra S_D, S_P, S_M can be calculated using the same method.

The power spectra S_R are calculated on first order frequencies using the same method for ω ∈ Ω

$S_R(\omega )=\frac{1}{N_{\Omega }(\omega )}{\displaystyle \sum _{m=1}^M}{S}_R^{(m)}(\omega )$

In particular, ω ∈ Ω implies ${\mathcal{N}}_{\Omega }(\omega )\ge 1$ so the denominator is not zero.

It should be noted that although different sets Ω^(m) and Ξ^(m) may share frequencies, such post-result frequency redundancy is not frequency overlap, with each individual measurement y^(m)(t) being the response to a stimulus without overlap.

2.5 Stochastic Markov simulations

Although the Hodgkin–Huxley model is empirical, it has led to various biophysical interpretations. The potassium ion channel is generally considered to have four independent identical subunits, each characterized by a two-state process

$0\stackrel{\alpha }{{\displaystyle \underset{\beta }{\rightleftharpoons }}}1$

where α_n(V) and β_n(V) are voltage-dependent transition rates given in the Hodgkin–Huxley model. All four subunits must be open for the channel to be open. If n denotes the gating variable for subunit activation, the total potassium conductance is proportional to n⁴.

This four-subunits interpretation suggests molecular-scale conformational changes that exceed the accuracy of the Hodgkin–Huxley model. However, this concept provides a theoretical interpretation of the fundamental nonlinearity of the neuron that is interesting for exploring equivalent forms of the Hodgkin–Huxley model.

The original Hodgkin–Huxley model is nonlinear in two ways: one is the use of an exponentiation for the gating variable, such as n⁴; the other is the voltage dependence of the rate constants. In the first case, an exponentiation is used to describe the delay in current response observed after a step change in membrane potential. In the second case, the use of the voltage clamp at a constant membrane potential V₀ has the effect of linearizing the differential equation $\frac{dn}{dt}={\alpha }_n(V_0)(1-n)-{\beta }_n(V_0)n$ since the rate constants α_n(V₀) and β_n(V₀) become constants.

The two-state process can be extended to multiple sequential states typical of higher order chemical relaxation models. More precisely, the four independent identical subunits of the potassium ion channel can be modeled as a five-state Markov chain

$\begin{array}{ll}0\stackrel{4{\alpha }_n}{{\displaystyle \underset{{\beta }_n}{\rightleftharpoons }}}1\stackrel{3{\alpha }_n}{{\displaystyle \underset{2{\beta }_n}{\rightleftharpoons }}}2\stackrel{2{\alpha }_n}{{\displaystyle \underset{3{\beta }_n}{\rightleftharpoons }}}3\stackrel{{\alpha }_n}{{\displaystyle \underset{4{\beta }_n}{\rightleftharpoons }}}4& (5)\end{array}$

where each state 0, 1, 2, 3, 4 corresponds to the number of open subunits at a given time. Each channel has the probability p_k(t) to be in state k. In particular, 1 = p₀ + p₁ + p₂ + p₃ + p₄. The channel is open when all four subunits are open, which corresponds to state 4 with probability p₄. When the system is in the state k, there are k open subunits and (4 − k) closed subunits. For k < 4, the transition k → k + 1 opens one of the (4 − k) closed subunits, which corresponds to the rate (4 − k) α _n. Similarly, for k > 0, the transition k → k − 1 closes one of the k open subunits, which corresponds to the rate kβ_n.

More generally, arbitrary rate constants could be chosen other than integer multiples of α_n and β_n. Numerous models of this type with more than two independent rate constants have been compared to the original Hodgkin–Huxley model, e.g. Vandenberg and Bezanilla (1991).

Stochastic Markov simulations were based on Gillespie (1977) and Goldwyn et al. (2011); Goldwyn and Shea-Brown (2011) considering a statistical population of N potassium ion channels and a time interval t ∈ [0;T] discretized by Δt. Let N_k(t) be the number of channels in state k at time t. In particular, N = N₀ + N₁ + N₂ + N₃ + N₄. In the limit, for a large population, the proportion of channels in state k tends to the corresponding probability

${\displaystyle \underset{N\to \infty }{lim}}\frac{N_k(t)}{N}=p_k(t)$

Let ΔN_k(t) be the variation of the number of channels in state k during Δt. Let ΔN_i→j(t) be the number of channels switching from state i to state j during Δt.

To ensure validity, it will always be assumed that p_k, N_k, ΔN_k, ΔN_i→j are zero for indices outside the range {0, 1, 2, 3, 4}.

The proportion of channels switching from state k to state k ± 1 during Δt is determined by the transition rates

$\{\begin{array}{l}\text{Δ}{N}_{k\to k+1}=(4-k){\alpha }_n{N}_k\text{Δ}t\\ \text{Δ}{N}_{k\to k-1}=k{\beta }_n{N}_k\text{Δ}t\end{array}$

Furthermore, the variation of the number of channels in state k during Δt corresponds to the number of transitions to state k minus the number of transitions from state k. More precisely

$\text{Δ}{N}_k=\text{Δ}{N}_{k-1\to k}+\text{Δ}{N}_{k+1\to k}-\text{Δ}{N}_{k\to k-1}-\text{Δ}{N}_{k\to k+1}$

Replacing expressions

$\begin{array}{r}\text{Δ}{N}_k=(5-k){\alpha }_n{N}_{k-1}\text{Δ}t+(k+1){\beta }_n{N}_{k+1}\text{Δ}t-k{\beta }_n{N}_k\text{Δ}t\\ -(4-k){\alpha }_n{N}_k\text{Δ}t\end{array}$

Dividing by N and taking the limit, a system of differential equations is obtained

$\frac{d{p}_k}{dt}=(5-k){\alpha }_n{p}_{k-1}+(k+1){\beta }_n{p}_{k+1}-k{\beta }_n{p}_k-(4-k){\alpha }_n{p}_k$

More explicitly, this gives the master equation

$\{\begin{array}{l}\frac{d{p}_0}{dt}=-4{\alpha }_n{p}_0+{\beta }_n{p}_1\\ \frac{d{p}_1}{dt}=4{\alpha }_n{p}_0-(3{\alpha }_n+{\beta }_n)p_1+2{\beta }_n{p}_2\\ \frac{d{p}_2}{dt}=3{\alpha }_n{p}_1-2({\alpha }_n+{\beta }_n)p_2+3{\beta }_n{p}_3\\ \frac{d{p}_3}{dt}=2{\alpha }_n{p}_2-({\alpha }_n+3{\beta }_n)p_3+4{\beta }_n{p}_4\\ \frac{d{p}_4}{dt}={\alpha }_n{p}_3-4{\beta }_n{p}_4\end{array}$

The conductance of the population of channels is calculated as g_Kf where g_K is the conductance of an individual channel and f is the proportion of open channels (Goldwyn et al., 2011). In the limit, for a large population, the proportion of open channels tends to the probability p₄ that a channel is open, which corresponds to the deterministic conductance g_Kp₄.

Although the master equation does not use exponentiation, the equivalence with the n⁴ model is discussed by Dayan and Abbott (2001). Indeed, in the Hodgkin–Huxley model, n is a gating variable representing the probability that a subunit is open and 1 − n the probability that it is closed. In state k, k of the four subunits are open and the 4 − k others are closed, thus $p_k=\left(\begin{array}{c}4\\ k\end{array}\right)n^k{(1-n)}^{4-k}$. In particular, in state 4, all four subunits are open, thus $p_4=n^4$ which is consistent with the exponentiation of the Hodgkin–Huxley model.

Therefore, at a fundamental level, nonlinearities generated by neuronal responses are likely to reflect fluctuations between internal states for which the master equation is controlled by transition rates involving energy for the movement of charges, the associated Boltzmann factor is discussed by Dayan and Abbott (2001). The master equation provides a deterministic description of the probabilities of these fluctuations.

Although Markov simulations are widely used to describe ion channel noise in the Hodgkin–Huxley model, alternative methods have been proposed, offering different trade-offs between accuracy and complexity. The Langevin approach, introduced by Fox and Lu (1994) and further refined by Fox (1997), incorporates stochasticity by adding Gaussian white noise terms to the equations, making it suitable for systems with a large number of channels. Goldwyn and Shea-Brown (2011) critically evaluated such stochastic formulations, highlighting the limitations of certain approximations in accurately representing neuronal behavior. Linaro et al. (2011) developed an improved diffusion approximation that enhances the accuracy of simulations, providing a more reliable alternative for neuronal modeling. Furthermore, Baravalle et al. (2017) presented a path integral approach to the Hodgkin–Huxley model, using methods from theoretical physics to analyze stochastic dynamics in which noise processes may be influenced by neural network interactions and feedforward correlations.

The advantage of Markov models lies in their ability to explicitly represent individual states of ion channels and the stochastic transitions between them, closely reflecting microscopic behavior. Unlike traditional Hodgkin–Huxley models, they inherently capture stochastic fluctuations in channel opening and closing, as transitions are governed by probabilistic processes. This enables Markov models to account for complex behaviors, including multiple transition pathways, nonlinear kinetics, and independent inactivation or recovery processes that are challenging to address with conventional methods. Notably, Markov models establish a connection between stochastic processes and nonlinear multiplicative effects, such as $p_4=n^4$, while Langevin approaches are based on the addition of noise. Furthermore, hierarchical Markov models (Siekmann et al., 2016) allow for the inclusion of modal gating of ion channels, capturing both transitions between modes and stochastic dynamics within modes. Overall, Markov models provide a scalable and versatile method for simulating neuronal phenomena, widely recognized as a standard in the field.

2.6 Markov power spectra

Stochastic analysis of Markov models provides power spectra of conductance noise identical to the first interpretation of Stevens (1972) involving probability and correlation functions. Many papers have derived noise power spectra based on the nonlinear nature of the potassium channel (n⁴), especially O'Donnell and van Rossum (2014) which is used in this article. In addition to the squid axon, measurements on the nodes of Ranvier (Conti, 1984) are consistent with the nonlinear origin of the conductance noise power spectra. Such an interpretation that certain nonlinear properties can induce a probabilistic or stochastic character, suggests that they are involved in the fundamental origin of spontaneous fluctuations in neurons.

The conductance noise can be predicted from the Hodgkin–Huxley equations for potassium conductance as described by O'Donnell and van Rossum (2014). Their approach is explored here for further adaptation in this article.

Let X_t be the random variable measuring if the channel is open (X_t = 1) or closed (X_t = 0) at time t. Let p(t) be the probability that the channel is open at time t, namely p(t) = p(X_t = 1). Let p_∞ be the steady state probability, which coincides with the average of p(t) over time. Actually, $p=p_4=n^4$ using the previous Markov probability notations, but the p notation is used here to be more general.

By definition, the autocorrelation r(t₁, t₂) characterizes the similarity of X_t between times t₁ and t₂. Denoting by E the expected value of a random variable

$r(t_1,t_2)=\text{E}\left[X_{t_1}{X}_{t_2}\right]$

Since X_t only takes the values 0 or 1, the autocorrelation is given by

$\begin{array}{l}r(t_1,t_2)={\displaystyle \sum _{i,j\in \left\{0,1\right\}}}i·j·p(X_{t_1}=i,X_{t_2}=j)\\ =p(X_{t_1}=1,X_{t_2}=1)\\ =p(X_{t_2}=1\mid X_{t_1}=1)p(X_{t_1}=1)\end{array}$

or more concisely

$r(t_1,t_2)=p(t_1)p(t_2\mid t_1)$

The autocorrelation can also be reformulated betwen times t₀ and t₀ + s to make the time lag s explicit

$r(t_0,t_0+s)=p(t_0)p(t_0+s\mid t_0)$

Here, p(t₀ + s ∣ t₀) is the conditional probability that the channel is open at time t₀ + s provided that it is open at time t₀. Importantly, if the channel is open at time t₀ then p(t₀) = 1. In this special case, the Markov state (p₀, p₁, p₂, p₃, p₄) = (0, 0, 0, 0, 1) is unique due to the constraint p₀ + p₁ + p₂ + p₃ + p₄ = 1. The same reasoning is valid for other sequential models with fewer or more Markov states. Thus, the time evolution of p(t₀ + s ∣ t₀) from the unique state (0, 0, 0, 0, 1) between t₀ and t₀ + s is identical to the time evolution of p(s ∣ 0) from the unique state (0, 0, 0, 0, 1) between 0 and s, the two trajectories of the dynamical system are identical because they start from the same unique state (0, 0, 0, 0, 1) and have the same duration s. Therefore

$r(t_0,t_0+s)=p(t_0)p(s\mid 0)$

Averaging over the time origin t₀ gives the autocorrelation as a function of the time lag

$\begin{array}{l}r(s)={\langle r(t_0,t_0+s)\rangle }_{t_0}\\ ={\langle p(t_0)p(s\mid 0)\rangle }_{t_0}\\ ={\langle p(t_0)\rangle }_{t_0}p(s\mid 0)\\ =p_{\infty }p(s\mid 0)\end{array}$

By definition, the autocovariance C(t₁, t₂) characterizes the similarity of X_t − E[X_t] between times t₁ and t₂. Namely

$\begin{array}{l}C(t_1,t_2)=\text{E}\left[(X_{t_1}-\text{E}\left[X_{t_1}\right])(X_{t_2}-\text{E}\left[X_{t_2}\right])\right]\\ =\text{E}\left[(X_{t_1}-p_{\infty })(X_{t_2}-p_{\infty })\right]\\ =\text{E}\left[X_{t_1}{X}_{t_2}\right]-p_{\infty }\text{E}\left[X_{t_1}\right]-p_{\infty }\text{E}\left[X_{t_2}\right]+p_{\infty }^2\\ =r(t_1,t_2)-p_{\infty }^2\end{array}$

The autocovariance can also be reformulated betwen times t₀ and t₀ + s to make the time lag s explicit

$C(t_0,t_0+s)=r(t_0,t_0+s)-p_{\infty }^2$

Averaging over the time origin t₀ gives the autocovariance as a function of the time lag

$\begin{array}{l}C(s)={\langle C(t_0,t_0+s)\rangle }_{t_0}\\ =r(s)-p_{\infty }^2\\ =p_{\infty }p(s\mid 0)-p_{\infty }^2\end{array}$

When the channel is open, the single-channel current is given by Ohm's law, namely i_K = γ_K(V − V_K) where γ_K is the single-channel conductance, V is the voltage imposed by voltage clamp and V_K is the reversal potential. Then, the fluctuating single-channel current is defined as i_KX_t. In particular, the fluctuating single-channel current is zero when the channel is closed and is given by Ohm's law when the channel is open. For a population of N_K channels, let $X_t^n$ be the random variables for each channel n = 1…N_K. The total fluctuating current is

$I_{\text{K}}(t)={\displaystyle \sum _{n=1}^{N_{\text{K}}}}{i}_{\text{K}}{X}_t^n$

By definition, the autocovariance C_IK(t₁, t₂) of the total fluctuating current between times t₁ and t₂ is given by

$\begin{array}{l}{C}_{I_{\text{K}}}\left(t_1,t_2\right)=\text{E}\left[\left(I_{\text{K}}\left(t_1\right)-\text{E}\left[I_{\text{K}}\left(t_1\right)\right]\right)\left(I_{\text{K}}\left(t_2\right)-\text{E}\left[I_{\text{K}}\left(t_2\right)\right]\right)\right]\\ =\text{E}[\left({\displaystyle \sum _{n=1}^{N_{\text{K}}}}{i}_{\text{K}}{X}_{t_1}^n-\text{E}\left[{\displaystyle \sum _{n=1}^{N_{\text{K}}}}{i}_{\text{K}}{X}_{t_1}^n\right]\right)({\displaystyle \sum _{m=1}^{N_{\text{K}}}}{i}_{\text{K}}{X}_{t_2}^m\\ -\text{E}\left[{\displaystyle \sum _{m=1}^{N_{\text{K}}}}{i}_{\text{K}}{X}_{t_2}^m\right]\left)\right]\\ =i_{\text{K}}^2\text{E}[\left({\displaystyle \sum _{n=1}^{N_{\text{K}}}}{X}_{t_1}^n-{\displaystyle \sum _{n=1}^{N_{\text{K}}}}\text{E}\left[X_{t_1}^n\right]\right)({\displaystyle \sum _{m=1}^{N_{\text{K}}}}{X}_{t_2}^m\\ {\displaystyle -\sum _{m=1}^{N_{\text{K}}}\text{E}[X_{t_2}^m]\left)\right]}\\ =i_{\text{K}}^2\text{E}\left[\left({\displaystyle \sum _{n=1}^{N_{\text{K}}}}\left(X_{t_1}^n-\text{E}\left[X_{t_1}^n\right]\right)\right)\left({\displaystyle \sum _{m=1}^{N_{\text{K}}}}\left(X_{t_2}^m-\text{E}\left[X_{t_2}^m\right]\right)\right)\right]\\ =i_{\text{K}}^2{\displaystyle \sum _{n=1}^{N_{\text{K}}}}{\displaystyle \sum _{m=1}^{N_{\text{K}}}}\text{E}\left[\left(X_{t_1}^n-p_{\infty }\right)\left(X_{t_2}^m-p_{\infty }\right)\right]\end{array}$

It can be noticed that

$\begin{array}{l}\text{E}\left[(X_{t_1}^n-p_{\infty })(X_{t_2}^m-p_{\infty })\right]=\text{E}\left[X_{t_1}^n{X}_{t_2}^m\right]-p_{\infty }\text{E}\left[X_{t_1}^n\right]-p_{\infty }\text{E}\left[X_{t_2}^m\right]+p_{\infty }^2\\ =\text{E}\left[X_{t_1}^n{X}_{t_2}^m\right]-p_{\infty }^2\end{array}$

If n ≠ m, the random variables $X_t^n$ and $X_t^m$ are independent, then

$\text{E}\left[(X_{t_1}^n-p_{\infty })(X_{t_2}^m-p_{\infty })\right]=\text{E}\left[X_{t_1}^n\right]\text{E}\left[X_{t_2}^m\right]-p_{\infty }^2=p_{\infty }{p}_{\infty }-p_{\infty }^2=0$

If n = m, the random variables $X_t^n$ and $X_t^m$ are the same as X_t, then

$\text{E}\left[(X_{t_1}^n-p_{\infty })(X_{t_2}^m-p_{\infty })\right]=\text{E}\left[X_{t_1}{X}_{t_2}\right]-p_{\infty }^2=r(t_1,t_2)-p_{\infty }^2=C(t_1,t_2)$

These remarks on the indices n and m allow to simplify the double sum in the autocovariance of the total fluctuating current so that

$\begin{array}{l}{C}_{I_{\text{K}}}(t_1,t_2)=i_{\text{K}}^2{\displaystyle \sum _{n=1}^{N_{\text{K}}}}C(t_1,t_2)\\ =N_{\text{K}}{i}_{\text{K}}^{2}C(t_1,t_2)\end{array}$

The autocovariance can also be reformulated betwen times t₀ and t₀ + s to make the time lag s explicit

$C_{I_{\text{K}}}(t_0,t_0+s)=N_{\text{K}}{i}_{\text{K}}^{2}C(t_0,t_0+s)$

Averaging over the time origin t₀ gives the autocovariance as a function of the time lag

$\begin{array}{l}{C}_{I_{\text{K}}}(s)={\langle N_{\text{K}}{i}_{\text{K}}^{2}C(t_0,t_0+s)\rangle }_{t_0}\\ =N_{\text{K}}{i}_{\text{K}}^2{\langle C(t_0,t_0+s)\rangle }_{t_0}\\ =N_{\text{K}}{i}_{\text{K}}^{2}C(s)\end{array}$

The Wiener–Khinchin theorem provides the power spectrum $S_{I_{\text{K}}}^{\pm }$ as the Fourier transform of the autocovariance C_IK, using the non-unitary Fourier transform with angular frequencies

$S_{I_{\text{K}}}^{\pm }(\omega )={\displaystyle {\int }_{-\infty }^{\infty }}{C}_{I_{\text{K}}}(s)e^{-i\omega s}ds$

The autocovariance is an even function of the lag because under stationarity, the average over the time origin t₀ is independent of a time lag

$\begin{array}{l}C(s)={\langle C(t_0,t_0+s)\rangle }_{t_0}\\ ={\langle C(t_0-s,t_0)\rangle }_{t_0-s}\\ ={\langle C(t_0-s,t_0)\rangle }_{t_0}\\ ={\langle C(t_0,t_0-s)\rangle }_{t_0}\\ =C(-s)\end{array}$

Thus, the power spectrum can be decomposed into two parts

$\begin{array}{l}{S}_{I_{\text{K}}}^{\pm }(\omega )={\displaystyle {\int }_{-\infty }^0}{C}_{I_{\text{K}}}(s)e^{-i\omega s}ds+{\displaystyle {\int }_0^{+\infty }}{C}_{I_{\text{K}}}(s)e^{-i\omega s}ds\\ ={\displaystyle {\int }_0^{+\infty }}{C}_{I_{\text{K}}}(-s)e^{i\omega s}ds+{\displaystyle {\int }_0^{+\infty }}{C}_{I_{\text{K}}}(s)e^{-i\omega s}ds\\ ={\displaystyle {\int }_0^{+\infty }}{C}_{I_{\text{K}}}(s)\left[e^{i\omega s}+e^{-i\omega s}\right]ds\\ =2{\displaystyle {\int }_0^{+\infty }}{C}_{I_{\text{K}}}(s)cos(\omega s)ds\end{array}$

To follow the convention of O'Donnell and van Rossum (2014), the power spectrum S_IK will be considered for positive frequencies only, i.e.

$S_{I_{\text{K}}}(\omega )=2S_{I_{\text{K}}}^{\pm }(\omega )=4\Re {\displaystyle {\int }_0^{+\infty }}{C}_{I_{\text{K}}}(s)e^{-i\omega s}ds$

In order to calculate the Markov power spectra for the fluctuating potassium current, it is necessary to determine the autocovariance $C(s)=p_{\infty }p(s\mid 0)-p_{\infty }^2$ based on the two components p(s∣0) and p_∞. The kinetics of the potassium gating variable n is given by

$\frac{dn}{dt}={\alpha }_n(1-n)-{\beta }_{n}n$

or equivalently

$\frac{dn}{dt}=\frac{n_{\infty }-n}{{\tau }_n}$

The general solution can be formulated with an exponential decay and an arbitrary initial condition n₀

$n(t)=n_{\infty }+(n_0-n_{\infty })e^{-t/{\tau }_n}$

The solution n(t∣0) is interpreted as the conditional probability that the gate is open at time t provided that it is open at time 0, which is implemented by the initial condition n₀ = 1

$n(t\mid 0)=n_{\infty }+(1-n_{\infty })e^{-t/{\tau }_n}$

The potassium ion channel having four independent identical subunits, this provides the probability that the channel is open at time t

$p(t)=n^4(t)$

At steady state

$p_{\infty }=n_{\infty }^4$

The conditional probability that the channel is open at time t provided that the channel is open at time t = 0 is given by

$\begin{array}{l}p(t\mid 0)=n^4(t\mid 0)\\ ={\left[n_{\infty }+(1-n_{\infty })e^{-t/{\tau }_n}\right]}^4\\ ={\displaystyle \sum _{q=0}^4\left(\genfrac{}{}{0.0pt}{}{4}{q}\right)}{n}_{\infty }^{4-q}{(1-n_{\infty })}^q{e}^{-qt/{\tau }_n}\end{array}$

The autocovariance does follow

$\begin{array}{l}{C}_{I_{\text{K}}}(t)=N_{\text{K}}{i}_{\text{K}}^{2}C(t)\\ =N_{\text{K}}{i}_{\text{K}}^2\left[p_{\infty }p(t\mid 0)-p_{\infty }^2\right]\\ =N_{\text{K}}{i}_{\text{K}}^2\left[n_{\infty }^{4}p(t\mid 0)-{(n_{\infty }^4)}^2\right]\\ =N_{\text{K}}{i}_{\text{K}}^2{n}_{\infty }^4\left[p(t\mid 0)-n_{\infty }^4\right]\\ =N_{\text{K}}{i}_{\text{K}}^2{n}_{\infty }^4{\displaystyle \sum _{q=1}^4\left(\genfrac{}{}{0.0pt}{}{4}{q}\right)}{n}_{\infty }^{4-q}{(1-n_{\infty })}^q{e}^{-qt/{\tau }_n}\end{array}$

where the term corresponding to q = 0 in the sum has been canceled with $-n_{\infty }^4$.

The power spectrum S_IK is computed by integrating each term of the sum in C_IK

$\begin{array}{l}{S}_q(\omega )={\displaystyle {\int }_0^{+\infty }\left(\genfrac{}{}{0.0pt}{}{4}{q}\right)}{n}_{\infty }^{4-q}{(1-n_{\infty })}^q{e}^{-qt/{\tau }_n}{e}^{-i\omega t}dt\\ ={\displaystyle \left(\genfrac{}{}{0.0pt}{}{4}{q}\right)}{n}_{\infty }^{4-q}{(1-n_{\infty })}^q{\displaystyle {\int }_0^{+\infty }}{e}^{-qt/{\tau }_n}{e}^{-i\omega t}dt\end{array}$

Using the non-unitary Fourier transform $\mathcal{F}\left[e^{-a\left|t\right|}\right]=\frac{2a}{a^2+{\omega }^2}$ with angular frequencies and divided by 2 because of the positive half of the time [0; + ∞]

$\begin{array}{l}{S}_q(\omega )={\displaystyle \left(\genfrac{}{}{0.0pt}{}{4}{q}\right)}{n}_{\infty }^{4-q}{(1-n_{\infty })}^q\frac{q/{\tau }_n}{q^2/{\tau }_n^2+{\omega }^2}\\ ={\displaystyle \left(\genfrac{}{}{0.0pt}{}{4}{q}\right)}{n}_{\infty }^{4-q}{(1-n_{\infty })}^q\frac{q{\tau }_n}{q^2+{\omega }^2{\tau }_n^2}\end{array}$

Applying S_q to each q

$\{\begin{array}{l}{S}_1(\omega )=4n_{\infty }^3(1-n_{\infty })\frac{{\tau }_n}{1+{\omega }^2{\tau }_n^2}\\ S_2(\omega )=6n_{\infty }^2{(1-n_{\infty })}^2\frac{2{\tau }_n}{4+{\omega }^2{\tau }_n^2}\\ S_3(\omega )=4n_{\infty }{(1-n_{\infty })}^3\frac{3{\tau }_n}{9+{\omega }^2{\tau }_n^2}\\ S_4(\omega )={(1-n_{\infty })}^4\frac{4{\tau }_n}{16+{\omega }^2{\tau }_n^2}\end{array}$

the power spectrum is obtained

$\begin{array}{ll}{S}_{I_{\text{K}}}(\omega )=4N_{\text{K}}{i}_{\text{K}}^2{n}_{\infty }^4\left[S_1(\omega )+S_2(\omega )+S_3(\omega )+S_4(\omega )\right]& (6)\end{array}$

There are four Lorentzians with corner frequencies ${\omega }_1=\frac{1}{{\tau }_n}$, ${\omega }_2=\frac{2}{{\tau }_n}$, ${\omega }_3=\frac{3}{{\tau }_n}$, ${\omega }_4=\frac{4}{{\tau }_n}$ respectively. In particular, jω_q correspond to the poles.

The parameters of Table 1 were reused with the potassium channel density ${\rho }_{\text{K}}=18/\text{μ}{\text{m}}^2$, the membrane area $A_{\text{K}}=500\text{μ}{\text{m}}^2$, the number of potassium channels N_K = ρ_KA_K and the potassium single-channel conductance γ_K = g_K/N_K.

Table 1

Table 1. Parameters of the simulation of the Hodgkin–Huxley equations used to compute impedances given by Equation 2 at different membrane potential displacements V₀.

Figure 2A illustrates a voltage clamp Markov simulation (Equation 5) for a 5 mV depolarization of the spontaneous potassium current fluctuations superimposed on the theoretically predicted power spectrum S_IK, the four Lorentzians, the four corner frequencies and the admittance. The simulation was done with the Markov model and QSA multi-sinusoidal stimuli of amplitude 0.25 mV. At this potential, the predicted power spectrum S_IK (magenta curve) have relative equivalent contributions from all the Lorentzian functions S₁, S₂, S₃, S₄ and certainly not just those present in the lowest corner frequency, ${\omega }_4=\frac{1}{{\tau }_n}$. The voltage responses at the stimulus frequencies are accurately described by the squared admittance |Ŷ|² (green curve) after adjustment of the vertical offset. The square of the admittance is obtained with an input of constant amplitude, which makes it possible to compare it to other curves independent of the input. Figure 2B shows that simulation done without stimulus gave identical power spectra.

Figure 2

Figure 2. Markov simulations of the n⁴ model. The power spectra approximated by 128 iterations are represented by the blue scatterplot. The predicted power spectrum S_IK (magenta curve) accurately fits the blue scatterplot. The four Lorentzians compose the predicted power spectrum with S₁ (brown curve), S₂ (orange curve), S₃ (yellow curve), S₄ (blue curve). The corresponding corner frequencies ω₁, ω₂, ω₃, ω₄ are represented by vertical lines. (A) Simulations for a 5 mV depolarization with QSA multi-sinusoidal stimulus of amplitude 0.25 mV and frequencies (2, 3, 10, 21, 35, 50, 76, 104, 134, 143, 223, 239, 285, 388, 405, 515, 564, 636, 815, 892, 982) Hertz. The squared admittance |Ŷ|² (green curve) matches the linear Markov responses after adjustment of the vertical offset. (B) Simulations for a 5 mV depolarization without stimulus. The fluctuation power spectra are the same as those in (A) but there is no admittance. (C) Simulations for a 55 mV depolarization without stimulus.

2.7 Validity

Although this article focuses on potassium channels, the methods are compatible with the full Hodgkin–Huxley model, including sodium channels.

Linear analysis allows the analytical linearization of conductance-based models around a steady state (Mauro, 1970) and their frequency-domain representation using the Fourier transform. This method efficiently captures the behavior of the system under small perturbations.

Quadratic sinusoidal analysis (QSA) investigates nonlinear responses in the frequency domain around a steady state. It supports multi-sinusoidal measurements in both conductance-based models and experimental data (Magnani and Moore, 2011). This method has also been applied beyond Hodgkin–Huxley, including neurons of the prepositus hypoglossi nucleus (Magnani et al., 2013) and neurons of the medial entorhinal cortex (Magnani et al., 2014).

Multi-sinusoidal power spectra extend QSA by averaging multiple measurements while preserving its general applicability. Although demonstrated here for potassium channels, this method is also applicable to other channels, such as sodium.

Stochastic Markov simulations are here considered for potassium channels to focus on the interplay between spontaneous fluctuations and nonlinear responses in a single ion channel type. However, this Markov modeling method is also well-documented for sodium channels (Destexhe and Rudolph-Lilith, 2012; Goldwyn et al., 2011; Goldwyn and Shea-Brown, 2011; O'Donnell and van Rossum, 2014) and continues to evolve in the literature (e.g., Ramlow and Lindner, 2024).

Markov power spectra, computed here for potassium channels, are a method that can be adapted to sodium channels by replacing n⁴ with m³h, as described by O'Donnell and van Rossum (2014).

By focusing on potassium channels, this article isolates specific dynamics. Extending to sodium channels in future work would provide complementary insights into neuronal behavior, particularly in understanding the integrated roles of multiple ion channel types.

3 Results

3.1 The n² model

Since chemical relaxation models generally obey the fluctuation-dissipation theorem (FDT), it is useful to compare their linear and nonlinear behaviors using both averaged quadratic (QSA) and noise (Markov) power spectra. This will be done with the potassium conductance of Hodgkin–Huxley equations, namely n⁴, and compared with more general relaxation models having fewer sequential states. A similar comparison could be done with the sodium conductance.

A variant of the Hodgkin–Huxley equations, the Frankenhaeuser and Huxley (1964) equations for myelinated nerve, uses a n² model for the potassium conductance, which can be simulated as a three-state Markov chain

$0\stackrel{2{\alpha }_n}{{\displaystyle \underset{{\beta }_n}{\rightleftharpoons }}}1\stackrel{{\alpha }_n}{{\displaystyle \underset{2{\beta }_n}{\rightleftharpoons }}}2$

Thus, the n² model has only two independent rate constants, α_n and β_n. More generally, such a n² model is a particular sequential kinetic model 0⇌1⇌2, arbitrarily called p₂ in this paper, having the four specific rate constants above with only two independent rate constants. However, a general p₂ model, which has four independent rate constants, will be shown to be a reasonable model for a potassium conductance. In particular, the rate constants of the p₂ model can be selected to show similar behavior to the Hodgkin–Huxley n⁴ model. In this case, the p₂ model is a model without exponentiation, which is not identical to the n² model, but approximates the Hodgkin–Huxley n⁴ model.

3.2 The p₂ model

The general p₂ model has four independent rate constants, two forward k₁, k₃, and two backward k₄, k₂, as follows:

$\begin{array}{ll}{p}_0\stackrel{k_1}{{\displaystyle \underset{k_2}{\rightleftharpoons }}}{p}_1\stackrel{k_3}{{\displaystyle \underset{k_4}{\rightleftharpoons }}}{p}_2& (7)\end{array}$

The n² model is a special case of the p₂ model with k₁ = 2α_n, k₃ = α_n, k₄ = 2β_n and k₂ = β_n.

More generally, it may be convenient to rewrite k₁, k₂, k₃, k₄ in terms of α_n and β_n such that k₁ = Aα_n, k₃ = α_n, k₄ = Bβ_n and k₂ = β_n, where A and B are arbitrary factors. In this way, the n² model is a special case of the p₂ model with A = B = 2.

In this article, the p₂ model simulations use A = 0.35 and B = 4. The other parameters are the same as the n⁴ model simulations given in Table 1.

The master equation is deduced from the three-state Markov chain as follows

$\begin{array}{ll}\{\begin{array}{l}\frac{d{p}_0}{dt}=-k_1{p}_0+k_2{p}_1\\ \frac{d{p}_1}{dt}=k_1{p}_0-(k_2+k_3)p_1+k_4{p}_2\\ \frac{d{p}_2}{dt}=k_3{p}_1-k_4{p}_2\end{array}& (8)\end{array}$

Since p₀ + p₁ + p₂ = 1, the system can be reduced to the variables p₁ and p₂

$\begin{array}{ll}\{\begin{array}{l}\frac{d{p}_1}{dt}=-(k_1+k_2+k_3)p_1+(-k_1+k_4)p_2+k_1\\ \frac{d{p}_2}{dt}=k_3{p}_1-k_4{p}_2\end{array}& (9)\end{array}$

This can be rewritten in matrix form

$\begin{array}{ll}\frac{d}{dt}\left(\begin{array}{l}{p}_1\\ p_2\end{array}\right)=\left(\begin{array}{ll}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right)\left(\begin{array}{l}{p}_1\\ p_2\end{array}\right)+\left(\begin{array}{c}-{\alpha }_{22}-{\alpha }_{12}\\ 0\end{array}\right)& (10)\end{array}$

where

$\{\begin{array}{l}{\alpha }_{11}=-(k_1+k_2+k_3)\\ {\alpha }_{12}=-k_1+k_4\\ {\alpha }_{21}=k_3\\ {\alpha }_{22}=-k_4\end{array}$

At steady state

$\begin{array}{l}\left(\begin{array}{c}{p}_1^{\infty }\\ p_2^{\infty }\end{array}\right)={\left(\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right)}^{-1}\left(\begin{array}{c}{\alpha }_{22}+{\alpha }_{12}\\ 0\end{array}\right)\\ =\frac{1}{{\alpha }_{11}{\alpha }_{22}-{\alpha }_{12}{\alpha }_{21}}\left(\begin{array}{cc}{\alpha }_{22}& -{\alpha }_{12}\\ -{\alpha }_{21}& {\alpha }_{11}\end{array}\right)\left(\begin{array}{c}{\alpha }_{22}+{\alpha }_{12}\\ 0\end{array}\right)& (11)\end{array}$

$\begin{array}{ll}{p}_1^{\infty }=\frac{{\alpha }_{22}({\alpha }_{22}+{\alpha }_{12})}{{\alpha }_{11}{\alpha }_{22}-{\alpha }_{12}{\alpha }_{21}},p_2^{\infty }=\frac{-{\alpha }_{21}({\alpha }_{22}+{\alpha }_{12})}{{\alpha }_{11}{\alpha }_{22}-{\alpha }_{12}{\alpha }_{21}}& (12)\end{array}$

3.3 Linear analysis of the p₂ model

The linear admittance of the p₂ model can be obtained using a method similar to that previously described for the Hodgkin–Huxley model by Mauro (1970). The matrix form given in Equation 10 represents the p₂ model of Equation 9 with k₁, k₂, k₃, k₄ replaced by α₁₁, α₁₂, α₂₁, α₂₂ as follows

$\begin{array}{ll}\{\begin{array}{l}\frac{d{p}_1}{dt}={\alpha }_{11}{p}_1+{\alpha }_{12}{p}_2-({\alpha }_{22}+{\alpha }_{12}),\\ \frac{d{p}_2}{dt}={\alpha }_{21}{p}_1+{\alpha }_{22}{p}_2\end{array}& (13)\end{array}$

The linearization for a small perturbation at steady state is given by

$\{\begin{array}{l}\delta \frac{d{p}_1}{dt}=\frac{d{\alpha }_{11}}{dV}\delta V{p}_1^{\infty }+{\alpha }_{11}\delta p_1+\frac{d{\alpha }_{12}}{dV}\delta V{p}_2^{\infty }\\ +{\alpha }_{12}\delta p_2-\frac{d({\alpha }_{22}+{\alpha }_{12})}{dV}\delta V\\ \delta \frac{d{p}_2}{dt}=\frac{d{\alpha }_{21}}{dV}\delta V{p}_1^{\infty }+{\alpha }_{21}\delta p_1+\frac{d{\alpha }_{22}}{dV}\delta V{p}_2^{\infty }+{\alpha }_{22}\delta p_2\end{array}$

or equivalently

$\{\begin{array}{l}\delta \frac{d{p}_1}{dt}-{\alpha }_{11}\delta p_1-{\alpha }_{12}\delta p_2=[\frac{d{\alpha }_{11}}{dV}{p}_1^{\infty }+\frac{d{\alpha }_{12}}{dV}{p}_2^{\infty }\\ {\displaystyle -\frac{d({\alpha }_{22}+{\alpha }_{12})}{dV}]\delta V}\\ \delta \frac{d{p}_2}{dt}-{\alpha }_{21}\delta p_1-{\alpha }_{22}\delta p_2=\left[\frac{d{\alpha }_{21}}{dV}{p}_1^{\infty }+\frac{d{\alpha }_{22}}{dV}{p}_2^{\infty }\right]\delta V\end{array}$

Applying the Fourier transform with angular frequency ω

$\{\begin{array}{l}(i\omega -{\alpha }_{11})\widehat{\delta p_1}-{\alpha }_{12}\widehat{\delta p_2}=[\frac{d{\alpha }_{11}}{dV}{p}_1^{\infty }+\frac{d{\alpha }_{12}}{dV}{p}_2^{\infty }\\ {\displaystyle -\frac{d({\alpha }_{22}+{\alpha }_{12})}{dV}]\widehat{\delta V}}\\ (i\omega -{\alpha }_{22})\widehat{\delta p_2}-{\alpha }_{21}\widehat{\delta p_1}=\left[\frac{d{\alpha }_{21}}{dV}{p}_1^{\infty }+\frac{d{\alpha }_{22}}{dV}{p}_2^{\infty }\right]\widehat{\delta V}\end{array}$

Then $\widehat{\delta p_2}$ can be deduced from

$\{\begin{array}{l}{\alpha }_{21}(i\omega -{\alpha }_{11})\widehat{\delta p_1}-{\alpha }_{21}{\alpha }_{12}\widehat{\delta p_2}={\alpha }_{21}\left[\frac{d{\alpha }_{11}}{dV}{p}_1^{\infty }+\frac{d{\alpha }_{12}}{dV}{p}_2^{\infty }-\frac{d({\alpha }_{22}+{\alpha }_{12})}{dV}\right]\widehat{\delta V}\\ (i\omega -{\alpha }_{11})(i\omega -{\alpha }_{22})\widehat{\delta p_2}-(i\omega -{\alpha }_{11}){\alpha }_{21}\widehat{\delta p_1}=(i\omega -{\alpha }_{11})\left[\frac{d{\alpha }_{21}}{dV}{p}_1^{\infty }+\frac{d{\alpha }_{22}}{dV}{p}_2^{\infty }\right]\widehat{\delta V}\end{array}$

By summation

$\begin{array}{l}\begin{array}{l}\left[(i\omega -{\alpha }_{11})(i\omega -{\alpha }_{22})-{\alpha }_{21}{\alpha }_{12}\right]\widehat{\delta p_2}\\ ={\alpha }_{21}\left[\frac{d{\alpha }_{11}}{dV}{p}_1^{\infty }+\frac{d{\alpha }_{12}}{dV}{p}_2^{\infty }-\frac{d({\alpha }_{22}+{\alpha }_{12})}{dV}\right]\widehat{\delta V}\\ +(i\omega -{\alpha }_{11})\left[\frac{d{\alpha }_{21}}{dV}{p}_1^{\infty }+\frac{d{\alpha }_{22}}{dV}{p}_2^{\infty }\right]\widehat{\delta V}\end{array}\end{array}$

This implies the rate of variation

$\begin{array}{l}\frac{\widehat{\delta p_2}}{\widehat{\delta V}}=\frac{\begin{array}{l}{\alpha }_{21}\left(\frac{d{\alpha }_{11}}{dV}{p}_1^{\infty }+\frac{d{\alpha }_{12}}{dV}{p}_2^{\infty }-\frac{d\left({\alpha }_{22}+{\alpha }_{12}\right)}{dV}\right)\\ +\left(i\omega -{\alpha }_{11}\right)\left(\frac{d{\alpha }_{21}}{dV}{p}_1^{\infty }+\frac{d{\alpha }_{22}}{dV}{p}_2^{\infty }\right)\end{array}}{\left(i\omega -{\alpha }_{11}\right)\left(i\omega -{\alpha }_{22}\right)-{\alpha }_{21}{\alpha }_{12}}\end{array}$

The current across the lipid bilayer is given by

$I=C_m\frac{dV}{dt}+g_{\text{L}}(V-V_{\text{L}})+g_{\text{K}}{p}_2(V-V_{\text{K}})$

It is linearized by

$\delta I=C_m\frac{d(\delta V)}{dt}+g_{\text{L}}\delta V+g_{\text{K}}\left[\delta p_2(V_0-V_{\text{K}})+p_2^{\infty }\delta V\right]$

Applying the Fourier transform with angular frequency ω

$\widehat{\delta I}=C_{m}i\omega \widehat{\delta V}+g_{\text{L}}\widehat{\delta V}+g_{\text{K}}\left[\widehat{\delta p_2}(V_0-V_{\text{K}})+p_2^{\infty }\widehat{\delta V}\right]$

This provides the admittance of the p₂ model

$\hat{Y}=\frac{\widehat{\delta I}}{\widehat{\delta V}}=C_{m}i\omega +g_{\text{L}}+g_{\text{K}}\left(\frac{\widehat{\delta p_2}}{\widehat{\delta V}}(V_0-V_{\text{K}})+p_2^{\infty }\right)$

3.4 Markov power spectra of the p₂ model

From Equations 10, 11, the differential equation can be written in the homogeneous form

$\frac{d}{dt}\left(\begin{array}{l}{p}_1-p_1^{\infty }\\ p_2-p_2^{\infty }\end{array}\right)=\left(\begin{array}{ll}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right)\left(\begin{array}{l}{p}_1-p_1^{\infty }\\ p_2-p_2^{\infty }\end{array}\right)$

More concisely

$\frac{d\text{p}}{dt}=\text{M}\text{p}$

where $\text{p}=\left(\begin{array}{c}{p}_1-p_1^{\infty }\\ p_2-p_2^{\infty }\end{array}\right)$ and $\text{M}=\left(\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right)$.

The general solution is a linear combination

$\begin{array}{ll}\text{p}=c_1{e}^{{\text{λ}}_{1}t}{\text{v}}_1+c_2{e}^{{\text{λ}}_{2}t}{\text{v}}_2& (14)\end{array}$

where λ₁, λ₂ are eigenvalues of M, v₁, v₂ are eigenvectors of M and c₁, c₂ are constants. Indeed

$\begin{array}{l}\text{M}\text{p}=c_1{e}^{{\text{λ}}_{1}t}\text{M}{\text{v}}_1+c_2{e}^{{\text{λ}}_{2}t}\text{M}{\text{v}}_2\\ =c_1{e}^{{\text{λ}}_{1}t}{\text{λ}}_1{\text{v}}_1+c_2{e}^{{\text{λ}}_{2}t}{\text{λ}}_2{\text{v}}_2\\ =\frac{d\text{p}}{dt}\end{array}$

The trace T = α₁₁ + α₂₂ and determinant D = α₁₁α₂₂ − α₁₂α₂₁ of M determine the eigenvalues

$\begin{array}{c}{\text{λ}}_{1,2}=\frac{T\pm \sqrt{T^2-4D}}{2}\\ {\text{λ}}_{1,2}=\frac{{\alpha }_{11}+{\alpha }_{22}\pm \sqrt{{({\alpha }_{11}+{\alpha }_{22})}^2-4({\alpha }_{11}{\alpha }_{22}-{\alpha }_{12}{\alpha }_{21})}}{2}\end{array}$

Since p represents probabilities, λ₁ and λ₂ must be negative to have exponential decreases rather than exponential increases. This determines two time constants for the p₂ model

${\tau }_1=-\frac{1}{{\text{λ}}_1},{\tau }_2=-\frac{1}{{\text{λ}}_2}$

The eigenvectors v₁ and v₂ are obtained by solving the equation Mv = λv and substituting eigenvalues λ = λ₁ and λ = λ₂

${\text{v}}_1=\left(\begin{array}{c}\frac{{\text{λ}}_1-{\alpha }_{22}}{{\alpha }_{21}}\\ 1\end{array}\right),{\text{v}}_2=\left(\begin{array}{c}\frac{{\text{λ}}_2-{\alpha }_{22}}{{\alpha }_{21}}\\ 1\end{array}\right)$

Inserting v₁ and v₂ into Equation 14

$\begin{array}{c}\left(\begin{array}{c}{p}_1-p_1^{\infty }\\ p_2-p_2^{\infty }\end{array}\right)=c_1{e}^{{\text{λ}}_{1}t}\left(\begin{array}{c}\frac{{\text{λ}}_1-{\alpha }_{22}}{{\alpha }_{21}}\\ 1\end{array}\right)+c_2{e}^{{\text{λ}}_{2}t}\left(\begin{array}{c}\frac{{\text{λ}}_2-{\alpha }_{22}}{{\alpha }_{21}}\\ 1\end{array}\right)\end{array}$

The time t = 0 determines the coefficients c₁ and c₂ as follows

$\begin{array}{c}{\displaystyle \{\begin{array}{l}{\displaystyle p_1\left(0\right)-p_1^{\infty }=c_1\frac{{\text{λ}}_1-{\alpha }_{22}}{{\alpha }_{21}}+c_2\frac{{\text{λ}}_2-{\alpha }_{22}}{{\alpha }_{21}}}\\ {\displaystyle p_2\left(0\right)-p_2^{\infty }=c_1+c_2}\end{array}}\\ {\displaystyle \{\begin{array}{l}{\displaystyle p_1\left(0\right)-p_1^{\infty }=c_1\frac{{\text{λ}}_1-{\alpha }_{22}}{{\alpha }_{21}}+\left(p_2\left(0\right)-p_2^{\infty }-c_1\right)\frac{{\text{λ}}_2-{\alpha }_{22}}{{\alpha }_{21}}}\\ {\displaystyle c_2=p_2\left(0\right)-p_2^{\infty }-c_1}\end{array}}\\ {\displaystyle \{\begin{array}{l}{\displaystyle c_1\left(\frac{{\text{λ}}_2-{\text{λ}}_1}{{\alpha }_{21}}\right)=p_2\left(0\right)\frac{{\text{λ}}_2-{\alpha }_{22}}{{\alpha }_{21}}-p_2^{\infty }\frac{{\text{λ}}_2-{\alpha }_{22}}{{\alpha }_{21}}-p_1\left(0\right)+p_1^{\infty }}\\ {\displaystyle c_2=p_2\left(0\right)-p_2^{\infty }-c_1}\end{array}}\\ {\displaystyle \{\begin{array}{l}{\displaystyle c_1=\frac{\left({\text{λ}}_2-{\alpha }_{22}\right)\left(p_2\left(0\right)-p_2^{\infty }\right)-{\alpha }_{21}\left(p_1\left(0\right)-p_1^{\infty }\right)}{{\text{λ}}_2-{\text{λ}}_1}}\\ {\displaystyle c_2=p_2\left(0\right)-p_2^{\infty }-c_1}\end{array}}\end{array}$

The power spectrum can be calculated using the same reasoning as with the model n⁴ previously because it is a sequential model with three states instead of five states. In particular, the single-channel autocovariance is based on the conditional probability p₂(t∣0) and the steady state probability $p_2^{\infty }$

$C(t)=p_2^{\infty }·p_2(t\mid 0)-{(p_2^{\infty })}^2$

Considering the unique state (p₀, p₁, p₂) = (0, 0, 1) when the channel is open

$\{\begin{array}{l}{c}_1=\frac{({\text{λ}}_2-{\alpha }_{22})(1-p_2^{\infty })+{\alpha }_{21}{p}_1^{\infty }}{{\text{λ}}_2-{\text{λ}}_1}\\ c_2=1-p_2^{\infty }-c_1\end{array}$

As a result, the conditional probability p₂(t|0) is given by

$p_2(t\mid 0)=p_2^{\infty }+c_1{e}^{{\text{λ}}_{1}t}+c_2{e}^{{\text{λ}}_{2}t}$

The autocovariance is given by

$\begin{array}{l}C(t)=p_2^{\infty }\left[p_2^{\infty }+c_1{e}^{{\text{λ}}_{1}t}+c_2{e}^{{\text{λ}}_{2}t}\right]-{(p_2^{\infty })}^2\\ =p_2^{\infty }{c}_1{e}^{{\text{λ}}_{1}t}+p_2^{\infty }{c}_2{e}^{{\text{λ}}_{2}t}\end{array}$

This implies autocovariance C_IK(t) of the total fluctuating current

$\begin{array}{l}{C}_{I_{\text{K}}}(t)=N_{\text{K}}{i}_{\text{K}}^{2}C(t)\\ =N_{\text{K}}{i}_{\text{K}}^2{p}_2^{\infty }(c_1{e}^{{\text{λ}}_{1}t}+c_2{e}^{{\text{λ}}_{2}t})\end{array}$

Following the same conventions as previously, the power spectrum is

$\begin{array}{l}{S}_{I_{\text{K}}}(\omega )=2S_{I_{\text{K}}}^{\pm }(\omega )\\ =4\Re {\displaystyle {\int }_0^{+\infty }}{C}_{I_{\text{K}}}(t)e^{-i\omega t}dt\\ =4N_{\text{K}}{i}_{\text{K}}^2{p}_2^{\infty }\Re {\displaystyle {\int }_0^{+\infty }}(c_1{e}^{{\text{λ}}_{1}t}+c_2{e}^{{\text{λ}}_{2}t})e^{-i\omega t}dt\end{array}$

There are two integrals to be calculated for each number q = 1, 2

$S_q(\omega )={\displaystyle {\int }_0^{+\infty }}{c}_q{e}^{{\text{λ}}_{q}t}{e}^{-i\omega t}dt$

Using the non-unitary Fourier transform $\mathcal{F}\left[e^{-a\left|t\right|}\right]=\frac{2a}{a^2+{\omega }^2}$ with angular frequencies and divided by 2 because of the positive half of the time [0; + ∞]

$S_q(\omega )=c_q\frac{-{\text{λ}}_q}{{\text{λ}}_q^2+{\omega }^2}$

Then, the power spectrum is deduced

$S_{I_{\text{K}}}\left(\omega \right)=4N_{\text{K}}{i}_{\text{K}}^2{p}_2^{\infty }\left(c_1\frac{-{\text{λ}}_1}{{\text{λ}}_1^2+{\omega }^2}+c_2\frac{-{\text{λ}}_2}{{\text{λ}}_2^2+{\omega }^2}\right)$

By using the time constants ${\tau }_1=-\frac{1}{{\lambda }_1}$ and ${\tau }_2=-\frac{1}{{\lambda }_2}$

$\begin{array}{l}{S}_{I_{\text{K}}}\left(\omega \right)=4N_{\text{K}}{i}_{\text{K}}^2{p}_2^{\infty }\left(c_1\frac{1/{\tau }_1}{1/{\tau }_1^2+{\omega }^2}+c_2\frac{1/{\tau }_2}{1/{\tau }_2^2+{\omega }^2}\right)\\ =4N_{\text{K}}{i}_{\text{K}}^2{p}_2^{\infty }\left(c_1\frac{{\tau }_1}{1+{\omega }^2{\tau }_1^2}+c_2\frac{{\tau }_2}{1+{\omega }^2{\tau }_2^2}\right)\end{array}$

Figure 3A illustrates a voltage clamp Markov simulation for a 5 mV depolarization of the spontaneous potassium current fluctuations superimposed on the theoretically predicted power spectrum S_IK, the two Lorentzians, the two corner frequencies and the squared admittance |Ŷ|². The simulation was done with the Markov model and QSA multi-sinusoidal stimuli of amplitude 0.25 mV. At this potential, the predicted power spectrum S_IK (magenta curve) is well described by a single Lorentzian function for frequencies greater than the corner frequency. The square of the admittance is obtained with an input of constant amplitude, which makes it possible to compare it to other curves independent of the input. Figure 3B shows that simulations done without stimulus gave identical power spectra.

Figure 3

Figure 3. Markov simulations of the p₂ model. The power spectra approximated by 128 iterations are represented by the blue scatterplot. The predicted power spectrum S_IK (magenta curve) accurately fits the blue scatterplot. The two Lorentzians compose the predicted power spectrum with S₁ (brown curve) and S₂ (orange curve). The corresponding corner frequencies ω₁ and ω₂ are represented by vertical lines. (A) Simulations for a 5 mV depolarization with QSA multi-sinusoidal stimulus of amplitude 0.25 mV and frequencies (2, 3, 10, 21, 35, 50, 76, 104, 134, 143, 223, 239, 285, 388, 405, 515, 564, 636, 815, 892, 982) Hertz. The squared admittance |Ŷ|² (green curve) matches the linear Markov responses after adjustment of the vertical offset. (B) Simulations for a 5 mV depolarization without stimulus. The fluctuation power spectra are the same as those in (A) but there is no admittance. (C) Simulations for a 55 mV depolarization without stimulus. The rate constants have been manually selected to approximate the frequency domain responses of the n⁴ model, which can be compared to Figure 2C.

3.5 Neuronal functions compared to neuronal fluctuations

Individual neurons are characterized both by deterministic neuronal functions, such as those revealed by linear and quadratic analysis (QSA), and by random fluctuations, such as those produced by Markov models. Both behaviors play an essential role in subthreshold neuronal processing to generate action potentials and, consequently, for neuronal networks in general.

Both behaviors are fundamentally related to the nonlinear kinetic processes underlying ion channels, such as the n⁴ model, its n² simplification, or the generalization of n² to the p₂ model without exponentiation, as described above. In contrast, the fluctuation-dissipation theorem implies that ion channel fluctuations should be described by a linearization of the underlying nonlinear channel kinetics. To address these issues, the sequential kinetic model p₂ described above is an approximate alternative to the Hodgkin–Huxley model that may be useful because it is not based on exponentiation. In all cases, Markov models are considered an adequate approach to simulate the fluctuations of a single-channel. Other neuronal noise fluctuation models are of interest, such as stochastic differential equations (SDEs), but are not considered in this paper. Since sequential kinetic models include Hodgkin–Huxley models with exponentiation as a special case, and have been shown to describe kinetic behavior as well as or better than the Hodgkin–Huxley model, it is appropriate to compare the simple p₂ model with the n⁴ model. This will be done in the following voltage clamp simulations for the model with nonlinear exponentiation n⁴ and the model without exponentiation p₂ to provide a comparison of Markov fluctuations, their linear and quadratic (QSA) stimulated behavior. A fundamental question is whether models with exponentiation of the Hodgkin–Huxley type or models without exponentiation are the best descriptions of a single neuron behavior. These fluctuation simulations determine whether or not the spontaneous current noise at a fixed voltage clamped membrane potential should be based on a nonlinear exponentiation such that n⁴.

The control case for nonlinear exponentiation is the Hodgkin–Huxley n⁴ model. The frequency domain analysis is performed on simulated data measurements at two different depolarized membrane potentials, namely 5 mV and 55 mV, such as one could record from a voltage clamped biological neuron. These two voltage clamp potentials were chosen to be either near the resting value or full activation of the voltage-dependent ionic conductances. Simulations were performed using MATLAB R2022b (MathWorks) with the nonlinear classical Hodgkin–Huxley n⁴ model described in above sections. Two modes have been considered, deterministic and stochastic, based on the programs of Goldwyn et al. (2011) and Goldwyn and Shea-Brown (2011). The deterministic mode was done using Euler's method for ODE, while the stochastic mode was based on Gillespie's algorithm. The results of the numerical simulations are superimposed on an analytical form for the impedance (Equation 2) and the Markov fluctuations (Equation 6). Thus, the simulations are analogous to measured data from a biological neuron, which can be analyzed in the frequency domain using the methods described above. As demonstrated by Magnani et al. (2013), this type of frequency domain analysis is an especially efficient and useful way to determine the accuracy of a particular model, since frequency domain analysis of data from biological neurons is a sensitive experimental measure that can be rigorously compared to model predictions. When data are produced by a model, as in this paper, comparisons between numerical simulations and analytical expressions of impedance and Markov fluctuations are used to check the accuracy of calculations. In addition, the QSA method provides a nonlinear analysis independent of a particular type of model or experiment, and thus can be applied to compare simulated data for different models.

Figure 4 illustrates the frequency analysis of the n⁴ model for a 5 mV depolarization. The upper left plot (Figure 4A) shows a typical low frequency linear admittance anti-resonance (smooth line) generated by a QSA multi-sinusoidal stimulus of amplitude 0.25 mV. The upper right plot (Figure 4B) shows the QSA matrix as a 3D representation, an intersection of lines on the plane for any two linear frequencies ω_i, ω_j represents an interactive quadratic frequency ω_i + ω_j for i, j ∈ {−N, …, −1, +1, …, +N} and the color code indicates the amplitude of the quadratic response. Remarkably, the quadratic response shows no anti-resonance. The lower left plot (Figure 4C) shows the power spectra of the Markov simulation superimposed on the analytical estimate of the Equation 6. It is clear that the analytical estimation provides an excellent control of the Markov simulation, both of which reveal the characteristics of a low-pass filter. The lower right plot (Figure 4D) shows the quadratic power spectra averaged over several ODE simulations for different random sets of QSA frequencies. Indeed, since the QSA frequency sets have few frequencies, it is necessary to average several measurements to increase the accuracy of the power spectra. Individual amplitudes of the various quadratic responses are represented at their particular frequencies, namely S_P(|ω_i| + |ω_j|) at frequency sums (red points), S_M(||ω_i| − |ω_j||) at frequency differences (blue points), S_D(2|ω_k|) at frequency doubling (orange points) and S_R(|ω_j|) for the mean squared quadratic output by matrix columns (black points). All of these responses exhibit low-pass filter characteristics, but they flatten out by reaching a constant value at high frequencies, with the exception of S_D(2|ω_k|) at frequency doubling for which the amplitudes decrease at high frequencies. Unlike the QSA matrix which is based on ratios of outputs to inputs, quadratic power spectra are evaluated directly from output measurements like Markov power spectra.

Figure 4

Figure 4. Frequency analysis of the n⁴ model for a 5 mV depolarization. (A) Linear admittance (smooth line) generated by a QSA multi-sinusoidal stimulus of amplitude 0.25 mV and frequencies (2, 3, 10, 21, 35, 50, 76, 104, 134, 143, 223, 239, 285, 388, 405, 515, 564, 636, 815, 892, 982) Hertz. The linear responses at stimulating frequencies are marked by small circles. (B) QSA matrix in 3D representation obtained with the same stimulus as the linear function, but the plot is cut at 500 Hz for a better readability. The color code indicates the amplitude of the quadratic response. Each value (ω_i, ω_j, |Q_{i, j}|) in the 3D plot represents the magnitude of the quadratic response to a frequency interaction. (C) Power spectra of the Markov simulation (128 iterations) superimposed on the analytical estimate. The frequencies are continuous up to twice the maximum stimulus frequency to include the highest QSA frequency (frequency doubling). The analytical estimation provides an excellent control of the Markov simulation. (D) Quadratic power spectra averaged over several ODE simulations (128 iterations) for different sets of QSA frequencies randomized up to 1,000 Hz. Quadratic responses are represented at their particular frequencies, namely S_P(|ω_i| + |ω_j|) at frequency sums (red points), S_M(||ω_i| − |ω_j||) at frequency differences (blue points), S_D(2|ω_k|) at frequency doubling (orange points) and S_R(|ω_j|) for the mean squared quadratic output by matrix columns (black points). Unlike the QSA matrix which is based on ratios of outputs to inputs, quadratic power spectra are evaluated directly from output measurements like Markov power spectra.

The amplitudes of the quadratic power spectra are smaller than those of the Markov power spectra. Indeed, quadratic responses are an order of magnitude smaller than linear responses, while the fluctuation-dissipation theorem relates spontaneous fluctuations to linear responses. In particular, Markov power spectra are based on autocorrelation, which corresponds to the second order cumulant, whereas higher order spectra would involve at least the third order cumulant as explained by Mendel (1991). At this end of this paper, the discussion provides some simulations to compare the amplitudes of quadratic responses and spontaneous fluctuations.

It is well known that nonlinear systems can generate frequency mixing processes, such as those producing interactive frequencies |ω_i| + |ω_j|, ||ω_i| − |ω_j|| and 2|ω_k|. Frequency mixing behavior is fundamental in some scientific fields, such as nonlinear optics, as described for example by Boyd (2008). Frequency mixing introduces complexity into the response, which may contain frequencies that are not present in the stimulus. Thus, neurons mix stimulus oscillations into quadratic responses that may appear less uniform than linear responses. In particular, the power spectra S_P(|ω_i| + |ω_j|) and S_M(||ω_i| − |ω_j||) show scatter plots with a lot of dispersion that look like stochastic fluctuations. It should be noted that although different random sets of QSA frequencies were used for averaging, each individual QSA matrix is obtained from a deterministic ODE. Randomizing the frequencies only extends the range of frequencies for analysis, implying that the dispersion is due to the quadratic neuronal function rather than a lack of frequencies. In contrast, frequency doubling generates power spectra S_D(2|ω_k|) that approximately follow a smooth line. This is because frequency doubling is calculated along the frequency diagonal (ω_k, ω_k) on the 3D representation, which decreases smoothly from low frequencies (color-coded red for height) to high frequencies (color-coded blue for height). Similarly, the power spectra S_R(|ω_j|) for the mean squared quadratic output by matrix columns appear to be quite smooth for the reason that the 3D representation is smooth and summed column by column. Remarkably, the power spectra S_D(2|ω_k|) tends to zero at high frequencies due to the attenuation of quadratic responses along the diagonal (ω_k, ω_k), while the power spectra S_R(|ω_j|) flatten at high frequencies because of the residual quadratic responses around the horizontal (x, 0) and vertical axes (0, y). Therefore, the irregularity of the power spectra S_P(|ω_i| + |ω_j|) and S_M(||ω_i| − |ω_j||) are fundamentally due to non-trivial quadratic frequency interactions between ω_i and ω_j when they do not follow the smooth shape of the 3D representation (by diagonal or by columns).

Figure 5A shows the superposition of five curves of the n⁴ model for a 5 mV depolarization. This allows a comparison of the linear, nonlinear and fluctuation amplitudes behavior. The curves were scaled to unity at low frequencies. The curve S_IK represents the analytical estimate of Markov fluctuations. The curve S_L represents the linear power spectrum at fundamental frequencies. The curve S_D represents the quadratic power spectrum at frequency doubling. The curve S_R represents the mean squared quadratic output by matrix columns. The multi-sinusoidal power spectra S_L, S_D, S_R were averaged over different random sets of non-overlapping frequencies as in Figure 4. The curve $|{Ŷ}_m{|}^2$ represents the squared admittance modified from the Equation 1 without the membrane capacitance nor the frequency independent terms, namely

${\hat{Y}}_m=g_{\text{K}}\left[4n_0^3\left(V_0-V_{\text{K}}\right)\widehat{D_n}\right]$

Figure 5

Figure 5. Power spectral analysis of the n⁴ model. Comparison between the analytical estimate of Markov fluctuations S_IK (magenta), the linear power spectrum S_L (green), the quadratic power spectra S_D (orange) and S_R (black), the modified squared admittance $|{Ŷ}_m{|}^2$ (blue). (A) Simulations for a 5 mV depolarization. (B) Simulations for a 55 mV depolarization.

The contribution of the capacitance to the actual admittance is in part responsible for an anti-resonance that leads to high frequency responses that are dominant, thus preventing a direct comparison of the linearized admittance with nonlinear responses or fluctuation behavior. The square of the admittance is obtained with an input of constant amplitude, which makes it possible to compare it to other curves independent of the input.

The fall of the modified squared admittance $|{Ŷ}_m{|}^2$ is more marked than that of the Markov fluctuations S_IK, which is also the case for the linear power spectrum S_L over a limited range of frequencies before the reversal. The fall of the quadratic power spectra S_D and S_R is more marked than that of the Markov fluctuations S_IK but less than that of the modified squared admittance, as well as compared to the linear power spectrum S_L before the reversal. As explained previously, S_D and S_R follow the smooth 3D representation of the QSA matrices, in particular S_D decreases indefinitely while S_R flattens at high frequencies.

These results show that the evoked nonlinear responses have a lower frequency behavior than spontaneous fluctuations, suggesting that Markov fluctuations have a different origin than the nonlinearity evoked by stimulus despite nonlinearity of rate constants or exponentiation n⁴ in the stochastic equations.

Figure 6 shows the frequency analysis of the n⁴ model for a 55 mV depolarization, which is very different from the previous 5 mV depolarization, clearly illustrating the voltage dependence of the potassium channel in both the linear and nonlinear analyses. The upper left plot shows linear admittance (smooth line) generated by a QSA multi-sinusoidal stimulus of amplitude 0.25 mV, which reveals the onset of a high frequency anti-resonance minimum. The upper right plot shows the QSA matrix as a 3D representation, which is increased in a slightly shifted bandwidth from the previous depolarization of 5 mV. The lower left plot shows the power spectra of the Markov simulation superimposed on the analytical estimate, which is similar to a low-pass filter as for a depolarization of 5 mV. The lower right plot shows the quadratic power spectra averaged over several ODE simulations for different random sets of QSA frequencies. The individual components S_P, S_M, S_D confirms the slightly shifted bandwidth observed on the QSA matrix, revealing more precisely the marked resonance. In particular, the nonlinear resonant frequencies are clearly lower than that of the linear anti-resonance. However, the S_R component based on each matrix column has low-pass filter characteristics. Thus, low-pass behavior is observed at 5 and 55 mV for the Markov simulation of and the column-mean-square S_R. Interestingly, the nonlinear behavior for highly activated conductances shows resonance in the quadratic responses at frequencies different from the anti-resonance observed for the linear response.

Figure 6

Figure 6. Frequency analysis of the n⁴ model for a 55 mV depolarization. The plots were generated using the same presentation as in Figure 4. (A) Linear function. (B) Quadratic function. (C) Markov power spectra (M = 128). (D) Quadratic power spectra (M = 128).

Figure 5B (55 mV) shows the superposition of five curves of the n⁴ model for a 55 mV depolarization, using the same presentation as in Figure 5A (5 mV). The modified squared admittance $|{Ŷ}_m{|}^2$ nearly superimposes on the Markov fluctuations S_IK, which is an effect of the voltage-dependent potassium conductance induced by the 55 mV depolarization.

Also, by comparing Figure 2C (55 mV) and Figure 2B (5 mV), each of the four relaxation time constants clearly appears as a function of the membrane potential and at large depolarizations the slowest time constant is dominant, i.e. the component S₁ fits the Markov fluctuations S_IK accurately.

Thus, at 55 mV depolarization, the fluctuation behavior S_IK is dominated by S₁ and has a kinetic behavior similar to that of a modified squared admittance $|{Ŷ}_m{|}^2$, which does not occur at the less depolarized potential 5 mV shown above.

The linear power spectrum S_L is similar to the modified squared admittance $|{Ŷ}_m{|}^2$ at low frequencies, but it does not fall at high frequencies due to the membrane capacitance and frequency independent terms. The column-mean-square component S_R has a higher frequency content than the linear behavior, although less dramatic than the marked scaled resonance of the frequency doubling component S_D.

In summary, the power spectra of the QSA responses for the n⁴ model are clearly different than those of Markov fluctuations and the power spectra of the linear responses do not generally follow that of the spontaneous fluctuations. Thus, Markov fluctuations are not predicted by small signal evoked linear or quadratic responses. Interestingly, Markov fluctuations at depolarized membrane potentials are dominated by the linear time constant, this is not generally the case for lesser depolarizations, as illustrated between a 5 and 55 mV depolarizations.

3.6 Observations on the linear fluctuation-dissipation theorem

As discussed by Stevens (1972), the power n⁴ in Hodgkin–Huxley equations leads to two interpretations of the origin of potassium conductance noise. The first case involves Markov simulations and power spectra with four Lorentzian terms, as illustrated in Figure 2B. The second case involves linearization of the Hodgkin–Huxley equations with a single relaxation time, as illustrated by the term $g_{\text{K}}\left[4n_0^3(V_0-V_{\text{K}})\widehat{D_n}+n_0^4\right]$ in Equation 1.

The linear fluctuation-dissipation theorem states that the spontaneous fluctuations should have the same kinetic behavior as a response to a small signal, namely a linear response. The fact that Markov simulations have four time constants whereas the linearization has a single relaxation time does not agree with the linear fluctuation-dissipation theorem. Thus, the linear fluctuation-dissipation theorem does not hold for the n⁴ model. This suggests exploring nonlinear extensions of the linear fluctuation-dissipation theorem in the physics literature. However, it is also interesting to consider the case when Markov simulations have fewer time constants, as with the p₂ model described in this article.

3.7 Comparison between n⁴ and p₂ models

3.7.1 Depolarization of 55 mV

As discussed above, sequential kinetic models can produce the exponentiation of the n⁴ Hodgkin–Huxley model, so that the probability of the channel opening coincides with n⁴. Similarly, for the n² model, the probability of the channel being open coincides with n². However, the p₂ model generalizes the rate constants of the n² model, so that the probability of the channel opening is not necessarily an exponentiation. Thus, sequential kinetic models are more general than the Hodgkin–Huxley model, as that they do not require exponentiation. In these models, the nonlinear voltage dependence of neuronal ionic conductances lies essentially in the rate constants between sequential states and not in an exponentiation functional dependence of the conductance gating variable.

If such non-exponentiation-based sequential kinetic models also fit the data, then their nonlinear behavior may more realistically describe the underlying molecular mechanisms of ion channel activity. Indeed, in accordance with the principle of parsimony, these models free themselves from the constraint that the probability of channel opening must be an exponentiation.

Figures 3C, 7 show an example of a p₂ model, in which the rate constants have been manually selected to approximate the frequency domain responses of the n⁴ model, for a depolarization of 55 mV. As expected, both models show similar resonance behavior for the linear and quadratic responses. However, there are quantitative differences. First, the S_D quadratic power spectra of the p₂ model (Figure 7) show not one, but two resonances reflected by a peak and bump, unlike the n⁴ model (Figure 6). Second, the Markov fluctuations of the p₂ model (Figure 3C) are also different with the indication of more than one time constant illustrated by the inflection just before the final slope of S₂, unlike the n⁴ model (Figure 2C) for which S₁ fits S_IK accurately.

Figure 7

Figure 7. Frequency analysis of the p₂ model for a 55 mV depolarization. The rate constants have been manually selected to approximate the frequency domain responses of the n⁴ model, which can be compared to Figure 6. (A) Linear function. (B) Quadratic function. (C) Markov power spectra (M = 128). (D) Quadratic power spectra (M = 128).

Figure 8B shows the superposition of five curves of the p₂ model for a 55 mV depolarization, using the same presentation as for the n⁴ model in Figure 5B. As with the n⁴ model, the modified squared admittance $|{Ŷ}_m{|}^2$ nearly superimposes on the Markov fluctuations S_IK. However, the power spectra are quite different between the two models. The inflection of the Markov fluctuations S_IK is matched by the modified squared admittance $|{Ŷ}_m{|}^2$, which is consistent with the fact that the analytical expressions for Markov fluctuations and admittance both have two relaxation processes in the p₂ model. In contrast, the n⁴ model has only one relaxation process due to the linearization of the exponentiation. Similarly, the modified squared admittance $|{Ŷ}_m{|}^2$ agrees well, at low frequencies, with the linear power spectrum S_L and column-mean-square component S_R. Interestingly, the frequency doubling component S_D exhibits a double resonance, suggesting the existence of more than one relaxation process.

Figure 8

Figure 8. Power spectral analysis of the p₂ model. Comparison between the analytical estimate of Markov fluctuations S_IK (magenta), the linear power spectrum S_L (green), the quadratic power spectra S_D (orange) and S_R (black), the modified squared admittance $|{Ŷ}_m{|}^2$ (blue). (A) Simulations for a 5 mV depolarization. (B) Simulations for a 55 mV depolarization.

3.7.2 Depolarization of 5 mV

Figure 9 shows the same simulations as in Figure 7 but for a depolarization of 5 mV. The results are quite similar to the n⁴ model in Figure 4. The upper left plot shows a typical low frequency linear admittance anti-resonance. The upper right plot shows the QSA matrix as a 3D representation, which shows no anti-resonance. The lower left plot shows the power spectra of the Markov simulation superimposed on the analytical estimate, both of which reveal the characteristics of a low-pass filter. The lower right plot shows the quadratic power spectra averaged over several ODE simulations, all of these responses exhibit low-pass filter characteristics, but they flatten out by reaching a constant value at high frequencies, with the exception of S_D(2|ω_k|) at frequency doubling for which the amplitudes decrease at high frequencies.

Figure 9

Figure 9. Frequency analysis of the p₂ model for a 5 mV depolarization. The rate constants have been manually selected to approximate the frequency domain responses of the n⁴ model, which can be compared to Figure 4. (A) Linear function. (B) Quadratic function. (C) Markov power spectra (M = 128). (D) Quadratic power spectra (M = 128).

Figure 8A shows that, as with the n⁴ model at 5 mV in Figure 5A, the linear behavior represented by the modified squared admittance $|{Ŷ}_m{|}^2$ does not accurately describe quadratic power spectra or Markov fluctuations. In particular, as with the n⁴ model at 5 mV, the modified squared admittance and quadratic power spectra of the p₂ model have lower frequency components than the Markov fluctuations.

3.8 Validity

Depolarizations at 5 and 55 mV were selected to investigate potassium channel dynamics under contrasting conditions for both the n⁴ model and the p₂ model. At 5 mV, $n_{\infty }^4\simeq 0.0247$ and $p_2^{\infty }\simeq 0.0297$ indicate negligible activation, reflecting the near-closed state of the channels. At 55 mV, $n_{\infty }^4\simeq 0.596$ and $p_2^{\infty }\simeq 0.565$ correspond to partial activation, illustrating the dynamics of channel opening under significant depolarization. These voltage levels were chosen to explore and compare the behavior of the two models across a physiologically relevant range under identical stimulation conditions.

4 Discussion

4.1 Spontaneous fluctuations and nonlinear responses

Neuronal fluctuations influence both transient and steady state responses, which fundamentally set the thresholds for neuronal impulse propagation. The fluctuating behavior of neurons in the central nervous system can be divided into three categories: (1) subthreshold synaptic bombardment from multiple networks, (2) additional nonlinear components generated for a particular stimulus, and (3) intrinsic spontaneous noise independent of external stimulus. Computational neuroscience is an attempt to understand the network behavior of the nervous system. However, it is highly dependent on neuronal models determined from numerous systems using very restricted data sets. The tools for extracting experimental data are impressive, ranging from current and voltage clamp, single channel and noise analysis, multi-electrode arrays recording activity on a population of cells, imaging techniques to estimate activity in intact systems, transfer function analysis and a variety of nonlinear approaches such as Wiener and Volterra kernels. The QSA method used in this article refines the quantitative characteristics of neuronal models in the frequency domain.

The above simulations shown for the potassium conductance suggest that the non-smooth fluctuations of the QSA power spectra S_P and S_M generate an alternative kind of noise due to the complexity of nonlinear interactions, which is not identical to the spontaneous fluctuations simulated by a stochastic Markov process. Although each individual QSA characterization uses a limited number of non-overlapping frequencies, this is overcome relatively well by averaging the power spectra for different stimulus sets. The p₂ sequential model of the ion channel suggests that the fluctuation-dissipation theorem may be, in part, valid since both the linear and fluctuation spectra have similar relaxation times (corner frequencies), which is not the case for the n⁴ model.

QSA analysis of a Markov model requires the averaging of individual trace measurements using the same multi-sinusoidal stimulus, which will converge to the deterministic response if a sufficient number of averages are performed. QSA analysis of individual trace measurements leads to highly variable QSA power spectra, unlike the same procedures performed on deterministic models, which are invariant. As real biological neuronal cells are intrinsically fluctuating, QSA experiments on patch clamped neurons have been performed on averaged trace measurements (Magnani and Moore, 2011), leading to a deterministic response that averages out the spontaneous fluctuation behavior of ion channels. The above-mentioned QSA power spectra were performed with random QSA stimulus frequencies applied to ODEs in order to compare them with the Markov noise power spectra. Alternatively, a single set of QSA stimulus frequencies can be applied to a set of Markov simulations, in which case the power spectra of the QSA matrix coefficients for the same stimulus can be averaged to determine how prominent the Markov model's responses are at the nonlinear interactive frequencies. This can be called an average QSA Markov noise power spectrum.

Figure 10 shows the power spectra of the QSA matrix coefficients of the p₂ model for a 55 mV depolarization, comparing two different stimulus amplitudes 4 mV (top) and 1 mV (bottom). Stimulus amplitudes are relatively larger than those used previously in deterministic analyses, as quadratic responses must overcome spontaneous fluctuations, otherwise the quality of noisy signals is too degraded. The two plots on the left column represent the power spectra of the usual deterministic ODE equations, which serve as a reference. As expected, they are similar at 4 and 1 mV. The two plots in the middle column represent the power spectra of the QSA Markov noise (QSA analysis performed on each individual Markov simulation). Surprisingly, the result is similar to that of ODE at 4 mV, but dramatically different at 1 mV. In particular, the enhanced diagonal follows harmonic frequencies. The right column shows the same power spectra as the middle column, but with the diagonal (harmonic frequencies) removed (set to zero). The result remains dramatically different from the ODE at 1 mV, which means that both the diagonal and cross-terms responses represent another quadratic function. Namely, spontaneous Markov fluctuations modify the neuronal quadratic function as the stimulus amplitude decreases. It is interesting to note that harmonic frequencies of the QSA Markov noise are also enhanced at 4 mV, but the matrix is essentially similar to that of ODE. The enhanced diagonal may be due to the fact that the modified squared admittance $|{Ŷ}_m{|}^2$ predicts the Markov fluctuations S_IK at a 55 mV depolarization, as shown in Figure 8B. Each squared frequency component generates frequency doubling. Thus, Markov fluctuations contribute to the harmonic frequencies in the QSA analysis of each individual Markov simulation, as well as at, to a lesser extent, to the interactive frequencies since the Markov noise cloud is not a smooth curve.

Figure 10

Figure 10. Power spectra of the QSA matrix coefficients of the p₂ model for a 55 mV depolarization. (Upper row) Stimulus amplitude 4 mV. (Lower row) Stimulus amplitude 1 mV. (Left column) Deterministic ODE equations. (Middle column) QSA analysis performed on each individual Markov simulation for 16 iterations. (Right column) Identical to the middle column but with the diagonal (harmonic frequencies) removed (set to zero).

Figure 11 supports this hypothesis as illustrated with the power spectra of the QSA matrix coefficients of the p₂ model for a 5 mV depolarization, comparing two different stimulus amplitudes 4 mV (top) and 1 mV (bottom). At this low depolarization level, the modified squared admittance $|{Ŷ}_m{|}^2$ no longer predicts the Markov fluctuations S_IK, as shown in Figure 8A. In this case, the diagonal of the QSA Markov noise is smaller and the power spectra of the QSA matrix coefficients are similar to those of the ODE at 4 mV and 1 mV. Nevertheless, the diagonal remains somewhat enhanced in the QSA Markov noise, likely because Markov fluctuations and modified squared admittance $|{Ŷ}_m{|}^2$ have a similar appearance.

Figure 11

Figure 11. Power spectra of the QSA matrix coefficients of the p₂ model for a 5 mV depolarization. (Upper row) Stimulus amplitude 4 mV. (Lower row) Stimulus amplitude 1 mV. (Left column) Deterministic ODE equations. (Middle column) QSA analysis performed on each individual Markov simulation for 16 iterations. (Right column) Identical to the middle column but with the diagonal (harmonic frequencies) removed (set to zero).

Figures 12–15 reinforce this interpretation by varying the surface area of the membrane of the p₂ model for a 55 mV depolarization and stimulus amplitude 1 mV. The larger the surface area of the membrane, the lower the noise effect. The smallest areas $A_{\text{K}}=50\text{μ}{\text{m}}^2$, $A_{\text{K}}=500\text{μ}{\text{m}}^2$, $A_{\text{K}}=5,000\text{μ}{\text{m}}^2$ show a behavior similar to Figure 10, that is to say spontaneous Markov fluctuations modify the neuronal quadratic function. However, the largest area $A_{\text{K}}=50,000\text{μ}{\text{m}}^2$ shows a behavior that partially recovers the ODE and better without the diagonal. This suggests that increasing the stimulus amplitude or decreasing the noise amplitude have similar effects on the quadratic function expressed. However, the two approaches are not equivalent since the amplitude of the stimulus tends to be modulated by the inputs of a neuron while the amplitude of the noise depends on the anatomy of a neuron.

Figure 12

Figure 12. Power spectra of the QSA matrix coefficients of the p₂ model for a membrane surface area $A_{\text{K}}=50\text{μ}{\text{m}}^2$, a 55 mV depolarization and a stimulus amplitude 1 mV. (Top left) Deterministic ODE equations. (Top right) A single trace (blue curve) from Markov simulations in time domain (s, pA) with linear (green curve) and quadratic (red curve) analyses. (Bottom left) QSA analysis performed on each individual Markov simulation for 16 iterations. (Bottom right) Identical to the bottom left but with the diagonal (harmonic frequencies) removed (set to zero).

Figure 13

Figure 13. Power spectra of the QSA matrix coefficients of the p₂ model for a membrane surface area $A_{\text{K}}=500\text{μ}{\text{m}}^2$, a 55 mV depolarization and a stimulus amplitude 1 mV. The plots were generated using the same presentation as in Figure 12.

Figure 14

Figure 14. Power spectra of the QSA matrix coefficients of the p₂ model for a membrane surface area $A_{\text{K}}=5,000\text{μ}{\text{m}}^2$, a 55 mV depolarization and a stimulus amplitude 1 mV. The plots were generated using the same presentation as in Figure 12.

Figure 15

Figure 15. Power spectra of the QSA matrix coefficients of the p₂ model for a membrane surface area $A_{\text{K}}=50,000\text{μ}{\text{m}}^2$, a 55 mV depolarization and a stimulus amplitude 1 mV. The plots were generated using the same presentation as in Figure 12.

This suggests that responses to neuronal stimuli would be different for quiet versus noisy neurons. Quiet neurons would exhibit complex nonlinear frequency responses, whereas the nonlinear responses of noisy neurons would tend to have additional frequency components, particularly at harmonic frequencies of the stimulus. Thus, individual neurons within a neuronal network would process input stimuli according to background synaptic activity, which is likely to be the main determinant of neuronal noise.

4.2 Potential applications

The spectral analysis of spontaneous fluctuations and nonlinear responses in ion channels, particularly potassium channels, provides an effective approach for understanding fundamental mechanisms of neuronal dynamics. These methods, which combine computational modeling and data analysis, are closely aligned with the principles of neuroinformatics by enabling the systematic interpretation of biological signals. This methodology suggests practical applications in areas relevant to real neuroinformatics scenarios.

The development of computational models of neurons provides a natural platform for applying spectral analysis, offering a way to connect ion channel dynamics with broader neuronal behavior. For example, Martin and Pedersen (2024) investigates the critical roles of HCN and M-type potassium channels in modeling and analyzing cAMP-induced mixed-mode oscillations in cortical neurons. Their use of a Hodgkin–Huxley-type model presents an opportunity to apply spectral analysis methods to characterize the behavior of different types of potassium channels in the context of cortical neurons.

The study by Pettersen et al. (2014) investigates power spectral densities (PSDs) of the soma potential, soma membrane current, and the single neuron contribution to the electroencephalogram (EEG). Their analytical approach, based on linear neuronal cable theory, explores how power laws may arise from the intrinsic biophysical properties of single neurons, independent of complex network interactions. Integrating spectral analysis methods within this theoretical basis could help clarify the interplay between ion channel dynamics and macroscopic neuronal activity patterns. Such integration holds potential for applications requiring accurate modeling of neuronal behavior, such as in EEG analysis or brain-computer interfaces.

Traumatic brain injuries and repetitive head impacts can result in altered neuronal excitability, with potassium channels playing an important role. Chapman et al. (2023) developed an in silico model of mouse CA1 pyramidal cells to identify ion channels contributing to hypoexcitability in this context. The model used morphology from the neuromorpho archive (Hines, 2009) and simulations performed with the NEURON platform (Ascoli et al., 2007). Applying spectral analysis methods to such realistic models could help identify spectral signatures of potassium channel dysfunction. Such analyses align with the goals of precision medicine and targeted therapeutic strategies.

Operator learning provides effective tools for modeling the complex dynamics of neurons. Centofanti et al. (2024) demonstrated the potential of operator learning, including Deep Operator Networks (DeepONets), Fourier Neural Operators (FNOs), and Wavelet Neural Operators (WNOs), to model the Hodgkin–Huxley system. Using square pulse stimuli, these methods map input currents to transmembrane potentials, capturing the dynamics of the Hodgkin–Huxley model. The square pulse stimuli could be extended to QSA multi-sinusoidal inputs to more precisely probe subthreshold nonlinear responses, although this may require some methodological adjustments. Extending their approach to a stochastic Markov model appears to be another challenge. Building on these ideas, spectral analysis methods could be combined with operator learning approaches to support applications in neuroscience and AI.

The principles of quantum mechanics provide innovative tools for understanding and modeling neuronal dynamics. The quantized Hodgkin–Huxley model proposed by Gonzalez-Raya et al. (2019) associates potassium conductance with a memristor driven by sinusoidal currents. Extending this method with QSA multi-sinusoidal driving could broaden its scope by probing nonlinear interactions across multiple frequencies. Additionally, Bradley et al. (2020) describe the use of quantum states to represent classical probability distributions over sets of sequences. Their use of a density operator to encode probabilistic information offers a potential approach to compare stochastic spontaneous fluctuations to nonlinear responses encoded in the QSA matrix, which is Hermitian. Together, these concepts offer promising directions for exploring quantum-inspired models of neuronal dynamics.

4.3 Comparing the p₂ model with existing approaches

The p₂ model is constructed as a variant of the n⁴ five-state Markov chain with two independent rate constants, formulated as a three-state Markov chain and generalized to four independent rate constants. In this way, the p₂ model avoids explicit exponentiation, yet it can represent the minimal degree of nonlinearity n² as well as models without exponentiation. Moreover, the rate constants of the p₂ model can be tuned to preserve the key dynamical features of the n⁴ model, as evaluated through spectral analysis methods.

The relevance of decreasing the degree of exponentiation to n² for potassium permeability is exemplified by the work of Frankenhaeuser and Huxley (1964) in their study of myelinated nerve fibers. This adaptation, motivated by experimental observations and computational simplicity, preserved the essential dynamical behaviors. Similarly, the p₂ model generalizes this approach by providing a flexible three-state Markov representation with independent rate constants.

The master equation of the p₂ model given in Equation 8 is consistent with established methods in the literature. For example, Güler (2015) describes a similar master equation in the context of the minimal diffusion formulation of Markov chain ensembles and its application to ion channel clusters. This connection highlights the relevance of such master equations in modeling ion channel dynamics and their stochastic behaviors.

Another example illustrating a similarity with the p₂ model formulation is provided by Schmandt and Galán (2012), who developed the stochastic-shielding approximation of Markov chains to efficiently simulate random ion-channel gating. They describe a three-state Markov chain and use variables N_k(t) and N_ij(t), representing the number of elements in state k and the number transitioning from state i to j, respectively, similar to the variables described in the above section on stochastic Markov simulations.

An insightful study is provided by Güler (2013), which investigates the impact of introducing noise into rate constants in Hodgkin–Huxley models. This approach modifies the standard rate constants α and β with stochastic terms, while the p₂ model applies fixed multiplicative adjustments to these rates. This work could inspire a potential extension of the p₂ model to incorporate stochastic variability in the rate constants.

4.4 Experimental considerations

The simulations and theoretical analyses presented in this article involve two key aspects for experimental validation: the deterministic QSA method and stochastic Markov models, both of which provide complementary insights into neuronal dynamics.

The QSA method was specifically designed to support both model-based and experimental measurements. In particular, the multi-sinusoidal stimuli are carefully generated to avoid frequency overlap at first and second orders, enabling a direct evaluation of the QSA matrix from input-output data. In this article, the QSA method has been extended by averaging multiple measurements to improve the accuracy of the power spectra. However, each individual measurement remains the response to a stimulus without overlap, ensuring that this extension does not alter the original experimental protocol. Consequently, the spectral analyses of linear and quadratic responses described in the above sections could potentially be applied to experimental measurements. Importantly, the QSA method has already been validated using experimental data from whole-cell patch-clamp recordings of prepositus hypoglossi nucleus neurons (Magnani and Moore, 2011; Magnani et al., 2013) and medial entorhinal cortex neurons (Magnani et al., 2014). Furthermore, these validated data were compared with conductance-based models of the corresponding neurons, emphasizing the relevance of the method for bridging experimental and theoretical approaches.

The validation of stochastic Markov models, such as the n⁴ and p₂ models, for describing ion channel kinetics presents a significant experimental challenge, as it inherently involves solving an inverse channel–fitting problem (Cannon and D'Alessandro, 2006). From a biophysical perspective, Markov models are considered more flexible than the Hodgkin–Huxley model because they incorporate discrete states that better capture the complexity of voltage-gated ion channels. For example, Andreozzi et al. (2019) compared Hodgkin–Huxley and Markov models for a sodium channel (Na_V1.5) with experimental data, indicating limitations of the Hodgkin–Huxley model and providing simplified Markov models as effective tools to approximate the complexity of ion channel kinetics.

The n⁴ exponentiation in the Hodgkin–Huxley model is commonly interpreted as describing a potassium ion channel composed of four independent gates. If each gate is open with probability n, the channel as a whole is open with probability n⁴. However, gating independence is merely an assumption. This nuance is addressed by Mannuzzu and Isacoff (2000), who show that the late gating steps of Shaker potassium channels involve cooperative interactions between subunits. Consistently, the p₂ model does not rely on the concept of independent gates. Instead, it relies on independent rate constants. This approach allows the p₂ model to represent behaviors such as n², as well as variants without explicit exponentiation. Moreover, this article has demonstrated that the p₂ and n⁴ models can produce similar neuronal responses, suggesting that the underlying kinetic schemes are more essential than gating independence for capturing channel behavior.

The p₂ model (or an extension) could potentially be further adapted to approximate experimental data as demonstrated here theoretically to approximate the n⁴ model. Power spectra could be compared to investigate discrepancies between fluctuation and QSA measurements (linear, quadratic). Earlier work by Fishman et al. (1983) compared fluctuation power spectra and linear admittance measurements of sodium current kinetics in squid axons, demonstrating significant discrepancies between the two methods and highlighting the nonlinear nature of microscopic sodium channel kinetics. This approach could shed light on the dynamics of ion channels and contribute to the validation of theoretical models.

Molecular dynamics simulations offer a promising approach to improve Markov models. Catacuzzeno et al. (2024) proposed an innovative method to build Markov models of ion channel permeation from molecular dynamics simulations. Combining these simulations with spectral analysis methods and the p₂ model remains a complex challenge but holds potential to achieve a more accurate representation of ion channel behavior under physiological conditions.

4.5 Implications for neuromorphic systems

Neuromorphic systems represent an innovative approach in neuroinformatics, based on biologically inspired principles to mimic dynamic processes of neuronal systems. By leveraging non-von Neumann architectures, they enable efficient modeling of complex neuronal functions. The potential of these systems has been demonstrated in various contexts: Sun et al. (2020) explored hippocampal dynamics using scalable digital neuromorphic models; Yang et al. (2021a) developed BiCoSS, a platform integrating multigranular spiking networks for advanced cognitive tasks; Yang et al. (2021b) applied neuromorphic principles to cerebellar motor learning. These studies collectively highlight the ability of neuromorphic systems to bridge biological relevance with computational efficiency, contributing to a deeper understanding of neuronal mechanisms.

The inherent complexity of neuronal systems, particularly the nonlinear dynamics of voltage-gated ion channels, poses significant challenges for neuromorphic systems, which must balance biological realism with computational efficiency. Neuronal processing is not limited to high-level brain functions but also encompasses computations at the level of individual neurons and voltage-gated ion channels, as characterized in this article through spectral methods such as linear impedance, quadratic interactions (QSA matrix), and microscopic noise (power spectra), which means that simplifying while preserving key aspects of neuronal computations in neuromorphic systems is a delicate task. For example, Sun et al. (2020) used the Izhikevich model (Izhikevich, 2003) in their neuromorphic implementation of hippocampal networks, achieving computational efficiency while simplifying the underlying biophysical complexity of ion channels. As pointed out in Yang et al. (2018), the Izhikevich model, while computationally efficient, does not fully capture the full range of ionic conductance dynamics, which may limit its use in studies requiring high biological accuracy. However, platforms such as BiCoSS (Yang et al., 2021a) demonstrate how neuromorphic systems can enhance biological representation and offer flexibility to address such challenges.

Biological neurons have highly specific properties, which present another significant challenge for neuromorphic systems to replicate their dynamics. As explained by Idoux et al. (2008), the prepositus hypoglossi nucleus (PHN) contains distinct neuronal populations, such as type B and type D neurons, which differ in properties like action potential shape and oscillatory behavior, highlighting the diversity within a single brain region. These neurons also exhibit varying electrotonic lengths: type B neurons are more compact, allowing distal dendritic inputs to significantly influence somatic activity, whereas type D neurons have extended electrotonic lengths, emphasizing the computational role of dendrites. Furthermore, the persistent sodium conductance (g_NaP), present in all PHN neurons, modulates activity differently depending on its somatic or dendritic localization, further illustrating the finely tuned specificity of neuronal mechanisms. Emulating multicompartment neurons with dendritic computations is essential for advancing the biological realism of neuromorphic systems, as highlighted by Yang et al. (2019). The platform BiCoSS (Yang et al., 2021a) supports multi-level, multigranule computing, including both point neuron and compartmental neuron. Furthermore, Beaubois et al. (2024) explored multicompartment Hodgkin–Huxley models in neuromorphic applications to better capture the complex morphological and functional properties of neurons.

Addressing the duality between membrane potential activity (continuous, analog) and spiking activity (discrete, event-driven) remains a significant challenge for neuromorphic systems, as these processes are intrinsically linked yet distinct in neuronal computation. Several methods have been developed to explore this duality. For example, Ris et al. (2001) showed that neurons in the medial vestibular nucleus exhibit resonance in the modulation of spike discharge, where spiking activity is closely tied to underlying membrane potential dynamics. Moreover, Recio-Spinoso et al. (2005) demonstrated that second-order Wiener kernels applied to spiking activity in auditory nerve fibers capture nonlinear interactions and temporal dynamics that first-order analyses cannot, providing deeper insights into the relationship between input sound stimuli and neuronal output. In another approach, Magnani et al. (2014) used quadratic sinusoidal analysis (QSA) in neurons of the medial entorhinal cortex to study how subthreshold nonlinearities influence spiking responses, highlighting the intricate link between membrane potential dynamics and suprathreshold activity. Additionally, Inoue et al. (2021) proposed a data-driven method to estimate latent variables and parameters of the Izhikevich neuron model using only spike-train data, employing the replica exchange particle-Gibbs with ancestor sampling (REPGAS) method. These methods are useful tools for understanding the interplay between subthreshold and suprathreshold dynamics and could provide insights for enhancing the biological realism of neuromorphic systems, as explored in Sun et al. (2020) and Yang et al. (2021a,b).

Maintaining eye movement stability in response to head or visual field movements relies on the precise coordination of neural mechanisms that integrate sensory inputs and sustain motor commands. Yang et al. (2021b) examined the optokinetic reflex (OKR) in the context of a neuromorphic cerebellar model, exploring how their cerebellar-inspired architecture adapts to retinal slip signals through mechanisms of supervised motor learning and synaptic plasticity. These mechanisms were applied in simulations of OKR adaptability to replicate dynamic responses like gain modulation and temporal synchronization in eye movement control. Building on biological insights, Magnani et al. (2013) investigated prepositus hypoglossi nucleus neurons, essential for maintaining eye position by integrating head velocity signals, highlighting the role of nonlinear dynamics driven by persistent sodium conductance (g_NaP) and active dendritic structures for sustaining the effects needed for gaze stabilization. Furthermore, Koulakov et al. (2002) introduced a robust model of the neural integrator based on bistable units, demonstrating how stability in eye position control can be achieved without fine-tuning of recurrent synaptic weights, which may be particularly relevant for neuromorphic systems. These studies show the importance of voltage-gated mechanisms and individual neuron dynamics in achieving eye movement stability, which have potential relevance for the design of neuromorphic systems that aim to replicate complex sensory-motor coordination.

Grid cells in the medial entorhinal cortex play a central role in spatial representation and navigation, providing a coordinate-like system to encode the position of an organism within its environment. These circuits are of particular relevance for neuromorphic systems due to their role in modeling spatial and cognitive functions. Krishna et al. (2021) represents a significant advancement in the hardware implementation of biologically inspired spatial navigation systems, presenting a biomimetic FPGA-based model that combines grid cells and place cells. However, the model abstracts away the biophysical details of individual neurons, like voltage-gated ion conductances, which are crucial for capturing the temporal dynamics of neurons, such as phase precession and resonance properties. From a biological perspective, Magnani et al. (2014) investigated the nonlinear properties of layer II stellate neurons in the rat medial entorhinal cortex (MEC) using quadratic sinusoidal analysis (QSA). The study revealed how dendritic filtering shapes somatic nonlinearities, with responses strongly influenced by the frequency content of stimuli. Persistent sodium conductance (g_NaP) was identified as a key source of nonlinearity, and near-threshold linear and nonlinear responses reliably predicted suprathreshold behavior, including action potential modulation. These findings emphasize the importance of biologically realistic modeling to connect single neuron voltage-gated mechanisms with scalable implementations of spatial navigation circuits. Thus, multigranular platforms like BiCoSS (Yang et al., 2021a) show particular promise for integrating such biophysical details into neuromorphic designs.

5 Conclusion

In conclusion, neurons clearly have linear and nonlinear responses to input stimuli. Linear responses are not just mirror images of the stimulus, they are clearly capable of enhancing certain frequencies due to linear resonance behavior. The linear response of resonance for the Hodgkin–Huxley potassium conductance is only present if the steady state value of n is voltage dependent, namely $\frac{d{n}_{\infty }}{d{V}_0}>0$, i.e. the resonance is a linear response that depends on the nonlinear property of the conductance. In general, nonlinear responses are considerably more complex, as shown by the presence of new interactive frequencies in the neuronal response that are not present in the input signal. The fluctuations present in neurons are due to voltage-dependent random mechanisms for which probabilities are controlled both linearly and nonlinearly and, as suggested above, ongoing synaptic activity can alter the nature of nonlinear responses (in the sense that QSA matrix can reflect fluctuations for small stimulus amplitudes).

Data availability statement

The raw simulation data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

CM: Investigation, Writing – original draft. LM: Investigation, Writing – original draft.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Acknowledgments

The authors thank the editor and the reviewers for their comments and suggestions on this work.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Andreozzi, E., Carannante, I., D'Addio, G., Cesarelli, M., and Balbi, P. (2019). Phenomenological models of NaV1.5. A side by side, procedural, hands-on comparison between Hodgkin–Huxley and kinetic formalisms. Sci. Rep. 9:17493. doi: 10.1038/s41598-019-53662-9

PubMed Abstract | Crossref Full Text | Google Scholar

Ascoli, G. A., Donohue, D. E., and Halavi, M. (2007). Neuromorpho.org: a central resource for neuronal morphologies. J. Neurosci. 27, 9247–9251. doi: 10.1523/JNEUROSCI.2055-07.2007

PubMed Abstract | Crossref Full Text | Google Scholar

Baravalle, R., Rosso, O. A., and Montani, F. (2017). A path integral approach to the Hodgkin–Huxley model. Phys A: Stat. Mech. Appl. 486, 986–999. doi: 10.1016/j.physa.2017.06.016

PubMed Abstract | Crossref Full Text | Google Scholar

Bear, M. F., Connors, B. W., and Paradiso, M. A. (2016). Neuroscience, 4th Edn. Philadelphia, PA: Wolters Kluwer.

Google Scholar

Beaubois, R., Cheslet, J., Ikeuchi, Y., Branchereau, P., and Levi, T. (2024). Real-time multicompartment Hodgkin–Huxley neuron emulation on SOC FPGA. Front. Neurosci. 18:1457774. doi: 10.3389/fnins.2024.1457774

PubMed Abstract | Crossref Full Text | Google Scholar

Boyd, R. W. (2008). Nonlinear Optics. Cambridge, MA: Academic Press.

Google Scholar

Bradley, T.-D., Stoudenmire, E. M., and Terilla, J. (2020). Modeling sequences with quantum states: a look under the hood. Mach. Learn. Sci. Technol. 1:035008. doi: 10.1088/2632-2153/ab8731

Crossref Full Text | Google Scholar

Cannon, R. C., and D'Alessandro, G. (2006). The ion channel inverse problem: Neuroinformatics meets biophysics. PLoS Comput. Biol. 2:e91. doi: 10.1371/journal.pcbi.0020091

PubMed Abstract | Crossref Full Text | Google Scholar

Catacuzzeno, L., Leonardi, M. V., Franciolini, F., Domene, C., Michelucci, A., Furini, S., et al. (2024). Building predictive Markov models of ion channel permeation from molecular dynamics simulations. Biophys. J. 123, 3832–3843. doi: 10.1016/j.bpj.2024.09.030

PubMed Abstract | Crossref Full Text | Google Scholar

Centofanti, E., Ghiotto, M., and Pavarino, L. F. (2024). Learning the Hodgkin–Huxley model with operator learning techniques. Comput. Methods Appl. Mech. Eng. 432:117381. doi: 10.1016/j.cma.2024.117381

Crossref Full Text | Google Scholar

Chapman, D. P., Vicini, S., Burns, M. P., and Evans, R. (2023). Single neuron modeling identifies potassium channel modulation as potential target for repetitive head impacts. Neuroinformatics 21, 501–516. doi: 10.1007/s12021-023-09633-7

PubMed Abstract | Crossref Full Text | Google Scholar

Cole, K. S. (1941). Rectification and Inductance in the Squid Giant Axon. J. Gen. Physiol. 25, 29–51. doi: 10.1085/jgp.25.1.29

PubMed Abstract | Crossref Full Text | Google Scholar

Conti, F. (1984). Noise analysis and single-channel recordings. Top. Membr. Transp. 22, 371–405. doi: 10.1016/S0070-2161(08)60478-5

Crossref Full Text | Google Scholar

Conti, F., and Wanke, E. (1975). Channel noise in nerve membranes and lipid bilayers. Q. Rev. Biophys. 8, 451–506. doi: 10.1017/S0033583500001967

PubMed Abstract | Crossref Full Text | Google Scholar

Dayan, P., and Abbott, L. F. (2001). Theoretical Neuroscience. Cambridge, MA: The MIT Press.

Google Scholar

Destexhe, A., and Rudolph-Lilith, M. (2012). Neuronal Noise. Cham: Springer US. doi: 10.1007/978-0-387-79020-6

Crossref Full Text | Google Scholar

Elinder, F., Frankenhaeuser, B., and Århem, P. (2001). Non-stationary fluctuation analysis of the Na current in myelinated nerve fibers of Xenopus laevis: experiments and stochastic simulations. BioSystems. 62, 13–28. doi: 10.1016/S0303-2647(01)00134-4

PubMed Abstract | Crossref Full Text | Google Scholar

Fishman, H., Leuchtag, H., and Moore, L. (1983). Fluctuation and linear analysis of Na-current kinetics in squid axon. Biophys. J. 43, 293–307. doi: 10.1016/S0006-3495(83)84353-7

PubMed Abstract | Crossref Full Text | Google Scholar

Fishman, H. M., Poussart, D. J. M., Moore, L. E., and Siebenga, E. (1977). K+ Conduction description from the low frequency impedance and admittance of squid axon. J. Membr. Biol. 32, 255–290. doi: 10.1007/BF01905222

PubMed Abstract | Crossref Full Text | Google Scholar

Fourier, J. B. J. (1822). Théorie Analytique de La Chaleur. Paris: Firmin Didot, père et fils.

Google Scholar

Fox, R. (1997). Stochastic versions of the Hodgkin–Huxley equations. Biophys. J. 72, 2068–2074. doi: 10.1016/S0006-3495(97)78850-7

PubMed Abstract | Crossref Full Text | Google Scholar

Fox, R. F., and Lu, Y.-n. (1994). Emergent collective behavior in large numbers of globally coupled independently stochastic ion channels. Phys. Rev. E 49, 3421–3431. doi: 10.1103/PhysRevE.49.3421

PubMed Abstract | Crossref Full Text | Google Scholar

Frankenhaeuser, B., and Huxley, A. F. (1964). The action potential in the myelinated nerve fibre of Xenopus laevis as computed on the basis of voltage clamp data. J. Physiol., 171, 302–315. doi: 10.1113/jphysiol.1964.sp007378

PubMed Abstract | Crossref Full Text | Google Scholar

Gillespie, D. T. (1977). Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81, 2340–2361. doi: 10.1021/j100540a008

Crossref Full Text | Google Scholar

Goldwyn, J. H., Imennov, N. S., Famulare, M., and Shea-Brown, E. (2011). Stochastic differential equation models for ion channel noise in Hodgkin–Huxley neurons. Phys. Rev. E 83:041908. doi: 10.1103/PhysRevE.83.041908

PubMed Abstract | Crossref Full Text | Google Scholar

Goldwyn, J. H., and Shea-Brown, E. (2011). The what and where of adding channel noise to the Hodgkin–Huxley equations. PLoS Comput. Biol. 7:e1002247. doi: 10.1371/journal.pcbi.1002247

PubMed Abstract | Crossref Full Text | Google Scholar

Gonzalez-Raya, T., Cheng, X.-H., Egusquiza, I. L., Chen, X., Sanz, M., Solano, E., et al. (2019). Quantized single-ion-channel Hodgkin–Huxley model for quantum neurons. Physi. Rev. Appl. 12:014037. doi: 10.1103/PhysRevApplied.12.014037

Crossref Full Text | Google Scholar

Güler, M. (2013). An investigation of the stochastic Hodgkin–Huxley models under noisy rate functions. Neural Comput. 25, 2355–2372. doi: 10.1162/NECO_a_00487

PubMed Abstract | Crossref Full Text | Google Scholar

Güler, M. (2015). Minimal diffusion formulation of Markov chain ensembles and its application to ion channel clusters. Phys. Rev. E 91:062116. doi: 10.1103/PhysRevE.91.062116

PubMed Abstract | Crossref Full Text | Google Scholar

Hines, M. (2009). Neuron and python. Front. Neuroinform. 3:2009. doi: 10.3389/neuro.11.001.2009

PubMed Abstract | Crossref Full Text | Google Scholar

Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. doi: 10.1113/jphysiol.1952.sp004764

PubMed Abstract | Crossref Full Text | Google Scholar

Idoux, E., Eugène, D., Chambaz, A., Magnani, C., White, J. A., and Moore, L. E. (2008). Control of neuronal persistent activity by voltage-dependent dendritic properties. J. Neurophysiol. 100, 1278–1286. doi: 10.1152/jn.90559.2008

PubMed Abstract | Crossref Full Text | Google Scholar

Inoue, H., Hukushima, K., and Omori, T. (2021). Estimation of neuronal dynamics of Izhikevich neuron models from spike-train data with particle Markov chain Monte Carlo method. J. Phys. Soc. Jpn. 90:104801. doi: 10.7566/JPSJ.90.104801

PubMed Abstract | Crossref Full Text | Google Scholar

Izhikevich, E. (2003). Simple model of spiking neurons. IEEE Trans. Neural Netw. 14, 1569–1572. doi: 10.1109/TNN.2003.820440

PubMed Abstract | Crossref Full Text | Google Scholar

Koulakov, A. A., Raghavachari, S., Kepecs, A., and Lisman, J. E. (2002). Model for a robust neural integrator. Nat. Neurosci. 5, 775–782. doi: 10.1038/nn893

PubMed Abstract | Crossref Full Text | Google Scholar

Krishna, A., Mittal, D., Virupaksha, S. G., Nair, A. R., Narayanan, R., Thakur, C. S., et al. (2021). Biomimetic FPGA-based spatial navigation model with grid cells and place cells. Neural. Netw. 139, 45–63. doi: 10.1016/j.neunet.2021.01.028

PubMed Abstract | Crossref Full Text | Google Scholar

Linaro, D., Storace, M., and Giugliano, M. (2011). Accurate and fast simulation of channel noise in conductance-based model neurons by diffusion approximation. PLoS Comput. Biol. 7:e1001102. doi: 10.1371/journal.pcbi.1001102

PubMed Abstract | Crossref Full Text | Google Scholar

Magnani, C., Economo, M. N., White, J. A., and Moore, L. E. (2014). Nonlinear properties of medial entorhinal cortex neurons reveal frequency selectivity during multi-sinusoidal stimulation. Front. Cell. Neurosci. 8:239. doi: 10.3389/fncel.2014.00239

PubMed Abstract | Crossref Full Text | Google Scholar

Magnani, C., Eugène, D., Idoux, E., and Moore, L. E. (2013). Vestibular integrator neurons have quadratic functions due to voltage dependent conductances. J. Comput. Neurosci. 35, 243–259. doi: 10.1007/s10827-013-0451-y

PubMed Abstract | Crossref Full Text | Google Scholar

Magnani, C., and Moore, L. E. (2011). Quadratic sinusoidal analysis of voltage clamped neurons. J. Comput. Neurosci. 31, 595–607. doi: 10.1007/s10827-011-0325-0

PubMed Abstract | Crossref Full Text | Google Scholar

Mannuzzu, L. M., and Isacoff, E. Y. (2000). Independence and cooperativity in rearrangements of a potassium channel voltage sensor revealed by single subunit fluorescence. J. Gen. Physiol. 115, 257–268. doi: 10.1085/jgp.115.3.257

PubMed Abstract | Crossref Full Text | Google Scholar

Marmarelis, P. Z., and Naka, K.-I. (1972). White-noise analysis of a neuron chain: an application of the Wiener theory. Science 175, 1276–1278. doi: 10.1126/science.175.4027.1276

PubMed Abstract | Crossref Full Text | Google Scholar

Marmont, G. (1949). Studies on the axon membrane. I. A new method. J. Cell. Comp. Physiol. 34, 351–382. doi: 10.1002/jcp.1030340303

PubMed Abstract | Crossref Full Text | Google Scholar

Martin, M., and Pedersen, M. G. (2024). Modelling and analysis of camp-induced mixed-mode oscillations in cortical neurons: critical roles of HCN and m-type potassium channels. PLoS Comput. Biol. 20:e1011559. doi: 10.1371/journal.pcbi.1011559

PubMed Abstract | Crossref Full Text | Google Scholar

Mauro, A. (1970). Subthreshold behavior and phenomenological impedance of the squid giant axon. J. Gen. Physiol. 55, 497–523. doi: 10.1085/jgp.55.4.497

PubMed Abstract | Crossref Full Text | Google Scholar

Mendel, J. M. (1991). Tutorial on higher-order statistics spectra in signal processing and system theory: theoretical results and some applications. Proc. IEEE 79, 278–305. doi: 10.1109/5.75086

Crossref Full Text | Google Scholar

Murphey, C. R., Moore, L. E., and Buchanan, J. T. (1995). Quantitative analysis of electrotonic structure and membrane properties of NMDA-activated lamprey spinal neurons. Neural Comput. 7, 486–506. doi: 10.1162/neco.1995.7.3.486

PubMed Abstract | Crossref Full Text | Google Scholar

O'Donnell, C., and van Rossum, M. C. W. (2014). Systematic analysis of the contributions of stochastic voltage gated channels to neuronal noise. Front. Comput. Neurosci. 8:105. doi: 10.3389/fncom.2014.00105

PubMed Abstract | Crossref Full Text | Google Scholar

Pettersen, K. H., Lindén, H., Tetzlaff, T., and Einevoll, G. T. (2014). Power laws from linear neuronal cable theory: Power spectral densities of the soma potential, soma membrane current and single-neuron contribution to the EEG. PLoS Comput. Biol. 10:e1003928. doi: 10.1371/journal.pcbi.1003928

PubMed Abstract | Crossref Full Text | Google Scholar

Poussart, D., and Ganguly, U. S. (1977). Rapid measurement of system kinetics–an instrument for real-time transfer function analysis. Proc. IEEE 65, 741–747. doi: 10.1109/PROC.1977.10556

Crossref Full Text | Google Scholar

Poussart, D. J. M. (1969). Nerve membrane current noise: direct measurements under voltage clamp. Proc. Nat. Acad. Sci. 64, 95–99. doi: 10.1073/pnas.64.1.95

PubMed Abstract | Crossref Full Text | Google Scholar

Ramlow, L., and Lindner, B. (2024). Noise intensity of a Markov chain. Phys. Rev. E 110:014139. doi: 10.1103/PhysRevE.110.014139

PubMed Abstract | Crossref Full Text | Google Scholar

Recio-Spinoso, A., Temchin, A. N., van Dijk, P., Fan, Y.-H., and Ruggero, M. A. (2005). Wiener-kernel analysis of responses to noise of chinchilla auditory-nerve fibers. J. Neurophysiol. 93, 3615–3634. doi: 10.1152/jn.00882.2004

PubMed Abstract | Crossref Full Text | Google Scholar

Ris, L., Hachemaoui, M., Vibert, N., Godaux, E., Vidal, P. P., Moore, L. E., et al. (2001). Resonance of spike discharge modulation in neurons of the guinea pig medial vestibular nucleus. J. Neurophysiol. 86, 703–716. doi: 10.1152/jn.2001.86.2.703

PubMed Abstract | Crossref Full Text | Google Scholar

Schmandt, N. T., and Galán, R. F. (2012). Stochastic-shielding approximation of Markov chains and its application to efficiently simulate random ion-channel gating. Phys. Rev. Lett. 109, 118101. doi: 10.1103/PhysRevLett.109.118101

PubMed Abstract | Crossref Full Text | Google Scholar

Siekmann, I., Fackrell, M., Crampin, E. J., and Taylor, P. (2016). Modelling modal gating of ion channels with hierarchical Markov models. Proc. R. Soc. A Math. Phys. Eng. Sci. 472:20160122. doi: 10.1098/rspa.2016.0122

PubMed Abstract | Crossref Full Text | Google Scholar

Sigworth, F. J. (1980). The variance of sodium current fluctuations at the node of Ranvier. J. Physiol. 307, 97–129. doi: 10.1113/jphysiol.1980.sp013426

PubMed Abstract | Crossref Full Text | Google Scholar

Stevens, C. F. (1972). Inferences about membrane properties from electrical noise measurements. Biophys. J. 12, 1028–1047. doi: 10.1016/S0006-3495(72)86141-1

PubMed Abstract | Crossref Full Text | Google Scholar

Sun, W., Wang, J., Zhang, N., and Yang, S. (2020). Scalable implementation of hippocampal network on digital neuromorphic system towards brain-inspired intelligence. Appl. Sci. 10:2857. doi: 10.3390/app10082857

Crossref Full Text | Google Scholar

Vandenberg, C. A., and Bezanilla, F. (1991). A sodium channel gating model based on single channel, macroscopic ionic, and gating currents in the squid giant axon. Biophys. J. 60, 1511–1533. doi: 10.1016/S0006-3495(91)82186-5

PubMed Abstract | Crossref Full Text | Google Scholar

Verveen, A. A., and Derksen, H. E. (1968). Fluctuation phenomena in nerve membrane. Proc. IEEE 56, 906–916. doi: 10.1109/PROC.1968.6443

Crossref Full Text | Google Scholar

Victor, J., and Shapley, R. (1980). A method of nonlinear analysis in the frequency domain. Biophys. J., 29, 459–483. doi: 10.1016/S0006-3495(80)85146-0

PubMed Abstract | Crossref Full Text | Google Scholar

Yang, S., Deng, B., Wang, J., Li, H., Lu, M., Che, Y., et al. (2019). Scalable digital neuromorphic architecture for large-scale biophysically meaningful neural network with multi-compartment neurons. IEEE Trans. Neural Netw. Learn. Syst. 31, 148–162. doi: 10.1109/TNNLS.2019.2899936

PubMed Abstract | Crossref Full Text | Google Scholar

Yang, S., Wang, J., Deng, B., Liu, C., Li, H., Fietkiewicz, C., et al. (2018). Real-time neuromorphic system for large-scale conductance-based spiking neural networks. IEEE Trans. Cybern. 49, 2490–2503. doi: 10.1109/TCYB.2018.2823730

PubMed Abstract | Crossref Full Text | Google Scholar

Yang, S., Wang, J., Hao, X., Li, H., Wei, X., Deng, B., et al. (2021a). Bicoss: toward large-scale cognition brain with multigranular neuromorphic architecture. IEEE Trans. Neural Netw. Learn. Syst. 33, 2801–2815. doi: 10.1109/TNNLS.2020.3045492

PubMed Abstract | Crossref Full Text | Google Scholar

Yang, S., Wang, J., Zhang, N., Deng, B., Pang, Y., Azghadi, M. R., et al. (2021b). Cerebellumorphic: large-scale neuromorphic model and architecture for supervised motor learning. IEEE Trans. Neural Netw. Learn. Syst. 33, 4398–4412. doi: 10.1109/TNNLS.2021.3057070

PubMed Abstract | Crossref Full Text | Google Scholar