<figure><img src="image_1.png" /></figure>

Check for updates

OPEN ACCESS

EDITED BY
Pawel Oswiecimka,
Polish Academy of Sciences, Poland

REVIEWED BY
Shuangming Yang,
Tianjin University, China
Anis Yuniati,
National Taiwan Normal University, Taiwan

*CORRESPONDENCE
Christophe Magnani
christophe.magnani@u-paris.fr

RECEIVED 29 July 2024
ACCEPTED 18 December 2024
PUBLISHED 15 January 2025

CITATION
Magnani C and Moore LE (2025) Power spectral analysis of voltage-gated channels in neurons. *Front. Neuroinform.* 18:1472499.
doi: 10.3389/fninf.2024.1472499

COPYRIGHT
© 2025 Magnani and Moore. This is an open-access article distributed under the terms of the [Creative Commons Attribution License (CC BY)](https://creativecommons.org/licenses/by/4.0/). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Power spectral analysis of voltage-gated channels in neurons

Christophe Magnani* and Lee E. Moore

Centre Borelli, Université Paris Cité, UMR 9010, CNRS, Paris, France

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_2$, 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_2$ (without exponentiation) and $n^4$ (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.

KEYWORDS

Hodgkin-Huxley, Markov, voltage-gated ion channels, neuronal noise, admittance, quadratic sinusoidal analysis

## 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 $\hat{u}(\omega)$ of the frequency variable $\omega$ where $\hat{u}$ 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_K = g_K^{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^4$) 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^4$ model is reduced to the minimal degree of nonlinearity $n^2$. Then, a novel sequential chemical state model named $p_2$ is introduced as a generalization of the $n^2$ model without exponentiation. Linear analysis, fluctuation simulations and theoretical power spectra are applied to the $p_2$ model by adapting the previous methods. Linear and quadratic functions are compared to neuronal fluctuations by adapting the previous methods. The $p_2$ and $n^4$ 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_2$ and $n^4$ 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_L - I_K - I_{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

$$\left\{ \begin{aligned} I_L &= g_L(V - V_L) \\ I_K &= g_K n^4(V - V_K) \\ I_{Na} &= g_{Na} m^3 h(V - V_{Na}) \end{aligned} \right.$$

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 $\alpha_i(V)$ and $\beta_i(V)$

$$\left\{ \begin{aligned} \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{aligned} \right.$$

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

$$\left\{ \begin{aligned} 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{aligned} \right.$$

where

$$\tau_n = (\alpha_n + \beta_n)^{-1}, \quad \tau_m = (\alpha_m + \beta_m)^{-1}, \quad \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_0 = (V_0, n_0, m_0, h_0)$ a steady state and $\frac{dX}{dt} = F(X)$ the system of differential equations. Clearly, $n_0 = n_{\infty}(V_0)$, $m_0 = m_{\infty}(V_0)$ and $h_0 = h_{\infty}(V_0)$ where $V_0$ 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) \sim \delta X^T DF(X_0)$$

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

$$\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\hat{\delta n} + (\alpha_n + \beta_n)\hat{\delta n} = \left[ \frac{d\alpha_n}{dV} - n_0 \frac{d(\alpha_n + \beta_n)}{dV} \right] \hat{\delta V}$$

This can be solved in the frequency domain as

$$\hat{\delta n} = \frac{\frac{d\alpha_n}{dV} - n_0 \frac{d(\alpha_n + \beta_n)}{dV}}{i\omega + \alpha_n + \beta_n} \hat{\delta V}$$

and similarly for $\hat{\delta m}$ and $\hat{\delta h}$. In order to simplify notations, this fraction will be denoted as $\hat{D}_n$, $\hat{D}_m$ and $\hat{D}_h$ respectively.

Similarly, linearization of the ionic currents is given by

$$\left\{ \begin{aligned} \delta I_L &= g_L \delta V \\ \delta I_K &= 4g_K n_0^3 (V_0 - V_K) \delta n + g_K n_0^4 \delta V \\ \delta I_{Na} &= 3g_{Na} m_0^2 h_0 (V_0 - V_{Na}) \delta m + g_{Na} m_0^3 (V_0 - V_{Na}) \delta h \\ &\quad + g_{Na} m_0^3 h_0 \delta V \end{aligned} \right.$$

which can be solved in the frequency domain as

$$\left\{ \begin{aligned} \delta \hat{I}_L &= g_L \delta \hat{V} \\ \delta \hat{I}_K &= g_K [4n_0^3(V_0 - V_K)\hat{D}_n + n_0^4]\delta \hat{V} \\ \delta \hat{I}_{Na} &= g_{Na}[3m_0^2h_0(V_0 - V_{Na})\hat{D}_m + m_0^3(V_0 - V_{Na})\hat{D}_h \\ &\quad + m_0^3h_0]\delta \hat{V} \end{aligned} \right.$$
Similarly, linearization of the current across the lipid bilayer is given by

$$C_m \frac{d\delta V}{dt} = \delta I - \delta I_L - \delta I_K - \delta I_{Na}$$

or equivalently in the frequency domain

$$C_m i \omega \delta \hat{V} = \hat{\delta I} - \hat{\delta I_L} - \hat{\delta I_K} - \hat{\delta I_{Na}}$$

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

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

$$\hat{Y} = C_m i \omega + g_L + g_K [4n_0^3(V_0 - V_K)\hat{D}_n + n_0^4] + g_{Na} [3m_0^2h_0(V_0 - V_{Na})\hat{D}_m + m_0^3(V_0 - V_{Na})\hat{D}_h + m_0^3h_0] \quad (1)$$

Equivalently, the impedance is defined as

$$\hat{Z} = \frac{1}{\hat{Y}} \qquad (2)$$

Note that the frequency $f$ in Hz is interchangeable with angular frequency $\omega = 2\pi 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_0$. 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 $|\hat{Z}(f, V_0)|$ 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.

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{dn_\infty}{dv} > 0$ for the potassium conductance [i.e. $n_\infty(V_0)$ does increase with depolarizing levels $V_0$]. 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 Quadratic sinusoidal analysis

Biological neurons and their models are fundamentally nonlinear systems that sometimes significantly contradict the linear superposition principle. System identification methods

<figure><img src="image_2.png" /><figcaption>FIGURE 1
Impedance functions for the Hodgkin-Huxley model. Different membrane potential displacements $V_0$ were applied from the resting potential of zero. Curves from top to bottom at lowest frequencies: -5, 0, 5, 10, 2, 15, 2, 25 in mV. Abscissa: frequency $f$ in kHz (logarithmic scale). Ordinate: magnitude $|\hat{Z}(f, V_0)|$ in $\Omega$ (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$).</figcaption></figure>
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 = 1$, $f_2 = 2$, $f_3 = 3$, $f_4 = 4$ (in Hz) generate unwanted overlaps because $f_1 + f_2 = f_3$ or $f_1 + f_4 = f_2 + f_3$. 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^2 = -1$ and $\omega$ is the angular frequency. More generally, a multi-sinusoidal signal of $N$ frequencies can be expressed as a sum of complex exponentials

$$x(t) = \sum_{k \in \Gamma} x_k e^{i\omega_k t}$$

where $\Gamma = \{-N, \dots, -1, 1, \dots, +N\}$ enumerates integers between $-N$ and $+N$ (zero excluded), $x_k$ are complex Fourier coefficients such that $x_{-k} = \bar{x}_k$ and $\omega_k$ are angular frequencies such that $\omega_{-k} = -\omega_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 = |x_k|e^{i\theta_k}$ where $|x_k|$ is the amplitude (nonnegative real number) and $\theta_k$ is the phase (between 0 and $2\pi$). 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_0$ is the DC component, $y_1(t)$ is the linear component and $y_2(t)$ is the quadratic component. More specifically

$$y(t) = y_0 + \sum_{k \in \Gamma} l_k x_k e^{i\omega_k t} + \sum_{i,j \in \Gamma} b_{i,j} x_i e^{i\omega_i t} x_j e^{i\omega_j t} \quad (3)$$

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 \in \Gamma$ so that all DC components are encoded in $y_0$.

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

$$\mathbf{x} = \sum_{k \in \Gamma} x_k \mathbf{e}_k$$

The corresponding time dependent vectors are defined by

$$\mathbf{x}_t = \sum_{k \in \Gamma} x_k e^{i\omega_k t} \mathbf{e}_k$$

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

$$y(t) = y_0 + \mathbf{L}\mathbf{x}_t + \mathbf{x}_t^T\mathbf{B}\mathbf{x}_t \quad (4)$$

By noticing that

$$\sum_{i,j \in \Gamma} b_{i,j} x_i e^{i\omega_i t} x_j e^{i\omega_j t} = \sum_{i,j \in \Gamma} b_{-i,j} \overline{x_i e^{i\omega_i t} x_j e^{i\omega_j t}}$$

we are led to define the QSA matrix $\mathbf{Q}$ by

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

so that

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

Remarkably, the QSA matrix $\mathbf{Q}$ is Hermitian (Magnani and Moore, 2011), which means that

$$\bar{\mathbf{Q}}^T = \mathbf{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 $\omega_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 $\omega_i + \omega_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 $|\hat{u}(\omega)|^2$ where $\hat{u}$ is the Fourier transform of u and $\omega$ 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 $|\omega_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 $|\omega_i| \neq |\omega_j|$:

$$S_P(|\omega_i| + |\omega_j|) = |\hat{y}_2(|\omega_i| + |\omega_j|)|^2$$

* Quadratic power spectrum at frequency differences for $|\omega_i| \neq |\omega_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} \sum_{i \in \Gamma} |Q_{i,j} \bar{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 $\Omega^{(m)}$ for $m = 1, \dots, M$. Each set of first order frequencies $\Omega^{(m)}$ determines a set of second order frequencies $\Xi^{(m)}$. The global sets of first and second order frequencies are defined by merging the individual sets, respectively

$$\left\{ \begin{array}{l} \Omega = \bigcup_{m=1}^{M} \Omega^{(m)} \\ \\ \Xi = \bigcup_{m=1}^{M} \Xi^{(m)} \end{array} \right.$$

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

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

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 $\omega \in \Omega^{(m)}$ and zero elsewhere. This allows to calculate the averaged linear power spectrum for $\omega \in \Omega$

$$S_L(\omega) = \frac{1}{\mathcal{N}_{\Omega}(\omega)} \sum_{m=1}^{M} |\hat{y}_1^{(m)}(\omega)|^2$$

where $\omega \in \Omega$ 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 $\xi \in \Xi^{(m)}$ and zero elsewhere. This allows to calculate the averaged quadratic power spectrum for $\xi \in \Xi$

$$S_2(\xi) = \frac{1}{\mathcal{N}_{\Xi}(\xi)} \sum_{m=1}^{M} |\hat{y}_2^{(m)}(\xi)|^2$$

where $\xi \in \Xi$ 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 $\omega \in \Omega$

$$S_R(\omega) = \frac{1}{\mathcal{N}_{\Omega}(\omega)} \sum_{m=1}^{M} S_R^{(m)}(\omega)$$

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

It should be noted that although different sets $\Omega^{(m)}$ and $\Xi^{(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 \underset{\beta}{\stackrel{\alpha}{\rightleftharpoons}} 1$$

where $\alpha_n(V)$ and $\beta_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^4$.

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^4$; 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_0$ 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 $\alpha_n(V_0)$ and $\beta_n(V_0)$ 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

$$0 \underset{\beta_n}{\stackrel{4\alpha_n}{\rightleftharpoons}} 1 \underset{2\beta_n}{\stackrel{3\alpha_n}{\rightleftharpoons}} 2 \underset{3\beta_n}{\stackrel{2\alpha_n}{\rightleftharpoons}} 3 \underset{4\beta_n}{\stackrel{\alpha_n}{\rightleftharpoons}} 4 \quad (5)$$

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_0 + p_1 + p_2 + p_3 + p_4$. The channel is open when all four subunits are open, which corresponds to state 4 with probability $p_4$. When the system is in the state $k$, there are $k$ open subunits and $(4-k)$ closed subunits. For $k < 4$, the transition $k \to k+1$ opens one of the $(4-k)$ closed subunits, which corresponds to the rate $(4-k)\alpha_n$. Similarly, for $k > 0$, the transition $k \to k-1$ closes one of the $k$ open subunits, which corresponds to the rate $k\beta_n$.

More generally, arbitrary rate constants could be chosen other than integer multiples of $\alpha_n$ and $\beta_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 \in [0; T]$ discretized by $\Delta t$. Let $N_k(t)$ be the number of channels in state $k$ at time $t$. In particular, $N = N_0 + N_1 + N_2 + N_3 + N_4$. In the limit, for a large population,

the proportion of channels in state $k$ tends to the corresponding probability

$$\lim_{N \to \infty} \frac{N_k(t)}{N} = p_k(t)$$

Let $\Delta N_k(t)$ be the variation of the number of channels in state $k$ during $\Delta t$. Let $\Delta N_{i \to j}(t)$ be the number of channels switching from state $i$ to state $j$ during $\Delta t$.

To ensure validity, it will always be assumed that $p_k$, $N_k$, $\Delta N_k$, $\Delta N_{i \to j}$ are zero for indices outside the range $\{0, 1, 2, 3, 4\}$.

The proportion of channels switching from state $k$ to state $k \pm 1$ during $\Delta t$ is determined by the transition rates

$$\begin{cases} \Delta N_{k \to k+1} = (4-k)\alpha_n N_k \Delta t \\ \Delta N_{k \to k-1} = k \beta_n N_k \Delta t \end{cases}$$

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

$$\Delta N_k = \Delta N_{k-1 \to k} + \Delta N_{k+1 \to k} - \Delta N_{k \to k-1} - \Delta N_{k \to k+1}$$

Replacing expressions

$$\Delta N_k = (5-k)\alpha_n N_{k-1} \Delta t + (k+1)\beta_n N_{k+1} \Delta t - k\beta_n N_k \Delta t - (4-k)\alpha_n N_k \Delta t$$

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

$$\frac{dp_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

$$\left\{ \begin{aligned} \frac{dp_0}{dt} &= -4\alpha_n p_0 + \beta_n p_1 \\ \frac{dp_1}{dt} &= 4\alpha_n p_0 - (3\alpha_n + \beta_n)p_1 + 2\beta_n p_2 \\ \frac{dp_2}{dt} &= 3\alpha_n p_1 - 2(\alpha_n + \beta_n)p_2 + 3\beta_n p_3 \\ \frac{dp_3}{dt} &= 2\alpha_n p_2 - (\alpha_n + 3\beta_n)p_3 + 4\beta_n p_4 \\ \frac{dp_4}{dt} &= \alpha_n p_3 - 4\beta_n p_4 \end{aligned} \right.$$

The conductance of the population of channels is calculated as $g_k f$ 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_4$ that a channel is open, which corresponds to the deterministic conductance $g_k p_4$.

Although the master equation does not use exponentiation, the equivalence with the $n^4$ 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 = \binom{4}{k} 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. Linoar 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^4$), 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_\infty$ 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_1, t_2)$ characterizes the similarity of $X_t$ between times $t_1$ and $t_2$. Denoting by $\mathbf{E}$ the expected value of a random variable

$$r(t_1, t_2) = \mathbf{E}[X_{t_1} X_{t_2}]$$

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

$$\begin{aligned}
r(t_1, t_2) &= \sum_{i,j \in \{0,1\}} i \cdot j \cdot p(X_{t_1} = i, X_{t_2} = j) \\
&= p(X_{t_1} = 1, X_{t_2} = 1) \\
&= p(X_{t_2} = 1 | X_{t_1} = 1)p(X_{t_1} = 1)
\end{aligned}$$

or more concisely

$$r(t_1, t_2) = p(t_1)p(t_2 | t_1)$$

The autocorrelation can also be reformulated between times $t_0$ and $t_0 + s$ to make the time lag $s$ explicit

$$r(t_0, t_0 + s) = p(t_0)p(t_0 + s | t_0)$$

Here, $p(t_0+s | t_0)$ is the conditional probability that the channel is open at time $t_0+s$ provided that it is open at time $t_0$. Importantly, if the channel is open at time $t_0$ then $p(t_0) = 1$. In this special case, the Markov state $(p_0, p_1, p_2, p_3, p_4) = (0, 0, 0, 0, 1)$ is unique due to the constraint $p_0 + p_1 + p_2 + p_3 + p_4 = 1$. The same reasoning is valid for other sequential models with fewer or more Markov states. Thus, the time evolution of $p(t_0 + s | t_0)$ from the unique state $(0, 0, 0, 0, 1)$ between $t_0$ and $t_0+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 | 0)$$

Averaging over the time origin $t_0$ gives the autocorrelation as a function of the time lag

$$\begin{aligned}
r(s) &= \langle r(t_0, t_0 + s) \rangle_{t_0} \\
&= \langle p(t_0)p(s | 0) \rangle_{t_0} \\
&= \langle p(t_0) \rangle_{t_0} p(s | 0) \\
&= p_\infty p(s | 0)
\end{aligned}$$
By definition, the autocovariance $C(t_1, t_2)$ characterizes the similarity of $X_t - E[X_t]$ between times $t_1$ and $t_2$. Namely

$$ 
\begin{aligned}
C(t_1, t_2) &= \mathbf{E}[(X_{t_1} - \mathbf{E}[X_{t_1}])(X_{t_2} - \mathbf{E}[X_{t_2}])] \\
&= \mathbf{E}[(X_{t_1} - p_{\infty})(X_{t_2} - p_{\infty})] \\
&= \mathbf{E}[X_{t_1}X_{t_2}] - p_{\infty}\mathbf{E}[X_{t_1}] - p_{\infty}\mathbf{E}[X_{t_2}] + p_{\infty}^2 \\
&= r(t_1, t_2) - p_{\infty}^2
\end{aligned}
$$

The autocovariance can also be reformulated between times $t_0$ and $t_0 + 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_0$ gives the autocovariance as a function of the time lag

$$ 
\begin{aligned}
C(s) &= \langle C(t_0, t_0 + s) \rangle_{t_0} \\
&= r(s) - p_{\infty}^2 \\
&= p_{\infty}p(s | 0) - p_{\infty}^2
\end{aligned}
$$

When the channel is open, the single-channel current is given by Ohm's law, namely $i_K = \gamma_K(V - V_K)$ where $\gamma_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_K X_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_K(t) = \sum_{n=1}^{N_K} i_K X_t^n $$

By definition, the autocovariance $C_{I_K}(t_1, t_2)$ of the total fluctuating current between times $t_1$ and $t_2$ is given by

$$ 
\begin{aligned}
C_{I_K}(t_1, t_2) &= \mathbf{E}[(I_{K}(t_1) - \mathbf{E}[I_{K}(t_1)])(I_{K}(t_2) - \mathbf{E}[I_{K}(t_2)])] \\
&= \mathbf{E}\left[\left(\sum_{n=1}^{N_K} i_K X_{t_1}^n - \mathbf{E}\left[\sum_{n=1}^{N_K} i_K X_{t_1}^n\right]\right)\left(\sum_{m=1}^{N_K} i_K X_{t_2}^m - \mathbf{E}\left[\sum_{m=1}^{N_K} i_K X_{t_2}^m\right]\right)\right] \\
&= i_K^2 \mathbf{E}\left[\left(\sum_{n=1}^{N_K} X_{t_1}^n - \sum_{n=1}^{N_K} \mathbf{E}[X_{t_1}^n]\right)\left(\sum_{m=1}^{N_K} X_{t_2}^m - \sum_{m=1}^{N_K} \mathbf{E}[X_{t_2}^m]\right)\right] \\
&= i_K^2 \mathbf{E}\left[\left(\sum_{n=1}^{N_K} (X_{t_1}^n - \mathbf{E}[X_{t_1}^n])\right)\left(\sum_{m=1}^{N_K} (X_{t_2}^m - \mathbf{E}[X_{t_2}^m])\right)\right] \\
&= i_K^2 \sum_{n=1}^{N_K} \sum_{m=1}^{N_K} \mathbf{E}[(X_{t_1}^n - p_{\infty})(X_{t_2}^m - p_{\infty})]
\end{aligned}
$$

It can be noticed that

$$ 
\begin{aligned}
\mathbf{E}[(X_{t_1}^n - p_{\infty})(X_{t_2}^m - p_{\infty})] &= \mathbf{E}[X_{t_1}^n X_{t_2}^m] - p_{\infty}\mathbf{E}[X_{t_1}^n] - p_{\infty}\mathbf{E}[X_{t_2}^m] + p_{\infty}^2 \\
&= \mathbf{E}[X_{t_1}^n X_{t_2}^m] - p_{\infty}^2
\end{aligned}
$$

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

$$ \mathbf{E}[(X_{t_1}^n - p_{\infty})(X_{t_2}^m - p_{\infty})] = \mathbf{E}[X_{t_1}^n]\mathbf{E}[X_{t_2}^m] - 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

$$ \mathbf{E}[(X_{t_1}^n - p_{\infty})(X_{t_2}^m - p_{\infty})] = \mathbf{E}[X_{t_1}X_{t_2}] - 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{aligned}
C_{I_K}(t_1, t_2) &= i_K^2 \sum_{n=1}^{N_K} C(t_1, t_2) \\
&= N_K i_K^2 C(t_1, t_2)
\end{aligned}
$$

The autocovariance can also be reformulated between times $t_0$ and $t_0 + s$ to make the time lag $s$ explicit

$$ C_{I_K}(t_0, t_0 + s) = N_K i_K^2 C(t_0, t_0 + s) $$

Averaging over the time origin $t_0$ gives the autocovariance as a function of the time lag

$$ 
\begin{aligned}
C_{I_K}(s) &= \langle N_K i_K^2 C(t_0, t_0 + s) \rangle_{t_0} \\
&= N_K i_K^2 \langle C(t_0, t_0 + s) \rangle_{t_0} \\
&= N_K i_K^2 C(s)
\end{aligned}
$$

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

$$ S_{I_K}^{\pm}(\omega) = \int_{-\infty}^{\infty} C_{I_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_0$ is independent of a time lag

$$ 
\begin{aligned}
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{aligned}
$$

Thus, the power spectrum can be decomposed into two parts

$$ 
\begin{aligned}
S_{I_K}^{\pm}(\omega) &= \int_{-\infty}^{0} C_{I_K}(s) e^{-i\omega s} ds + \int_{0}^{+\infty} C_{I_K}(s) e^{-i\omega s} ds \\
&= \int_{0}^{+\infty} C_{I_K}(-s) e^{i\omega s} ds + \int_{0}^{+\infty} C_{I_K}(s) e^{-i\omega s} ds \\
&= \int_{0}^{+\infty} C_{I_K}(s) [e^{i\omega s} + e^{-i\omega s}] ds \\
&= 2 \int_{0}^{+\infty} C_{I_K}(s) \cos(\omega s) ds
\end{aligned}
$$
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_{IK}(\omega) = 2S_{IK}^{\pm}(\omega) = 4\Re \int_{0}^{+\infty} C_{IK}(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 | 0) - p_{\infty}^2$ based on the two components $p(s | 0)$ and $p_{\infty}$. 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_0$

$$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_0 = 1$

$$n(t | 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{aligned} p(t | 0) &= n^4(t | 0) \\ &= [n_{\infty} + (1 - n_{\infty})e^{-t/\tau_n}]^4 \\ &= \sum_{q=0}^{4} \binom{4}{q} n_{\infty}^{4-q} (1 - n_{\infty})^q e^{-qt/\tau_n} \end{aligned}$$

The autocovariance does follow

$$\begin{aligned} C_{IK}(t) &= N_K I_K^2 C(t) \\ &= N_K I_K^2 [p_{\infty}p(t | 0) - p_{\infty}^2] \\ &= N_K I_K^2 [n_{\infty}^4 p(t | 0) - (n_{\infty}^4)^2] \\ &= N_K I_K^2 n_{\infty}^4 [p(t | 0) - n_{\infty}^4] \\ &= N_K I_K^2 n_{\infty}^4 \sum_{q=1}^{4} \binom{4}{q} n_{\infty}^{4-q} (1 - n_{\infty})^q e^{-qt/\tau_n} \end{aligned}$$

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{aligned} S_q(\omega) &= \int_0^{+\infty} \binom{4}{q} n_{\infty}^{4-q} (1 - n_{\infty})^q e^{-qt/\tau_n} e^{-i\omega t} dt \\ &= \binom{4}{q} n_{\infty}^{4-q} (1 - n_{\infty})^q \int_0^{+\infty} e^{-qt/\tau_n} e^{-i\omega t} dt \end{aligned}$$

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

$$\begin{aligned} S_q(\omega) &= \binom{4}{q} n_{\infty}^{4-q} (1 - n_{\infty})^q \frac{\tau_n}{q^2/\tau_n^2 + \omega^2} \\ &= \binom{4}{q} n_{\infty}^{4-q} (1 - n_{\infty})^q \frac{q\tau_n}{q^2 + \omega^2\tau_n^2} \end{aligned}$$

Applying $S_q$ to each $q$

$$\left\{ \begin{aligned} 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{aligned} \right.$$

the power spectrum is obtained

$$S_{IK}(\omega) = 4N_K I_K^2 n_{\infty}^4 [S_1(\omega) + S_2(\omega) + S_3(\omega) + S_4(\omega)] \quad (6)$$

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\omega_q$ correspond to the poles.

The parameters of Table 1 were reused with the potassium channel density $\rho_K = 18/\mu\text{m}^2$, the membrane area $A_K = 500\mu\text{m}^2$, the number of potassium channels $N_K = \rho_K A_K$ and the potassium single-channel conductance $\gamma_K = g_K/N_K$.

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_1$, $S_2$, $S_3$, $S_4$ 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 $|\hat{Y}|^2$ (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><img src="image_3.png" /></figure>

## 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.
<table><caption>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₀.</caption><tr><td>C<sub>m</sub></td><td>1 µF/cm<sup>2</sup></td></tr><tr><td>g<sub>L</sub></td><td>0.3 mS/cm<sup>2</sup></td></tr><tr><td>V<sub>L</sub></td><td>10.6 mV</td></tr><tr><td>g<sub>K</sub></td><td>36 mS/cm<sup>2</sup></td></tr><tr><td>V<sub>K</sub></td><td>-12 mV</td></tr><tr><td>g<sub>Na</sub></td><td>120 mS/cm<sup>2</sup></td></tr><tr><td>V<sub>Na</sub></td><td>120 mV</td></tr><tr><td>α<sub>n</sub>(V)</td><td>0.01 · (10 − V) / (e<sup>(10.001 − V)/10</sup> − 1)</td></tr><tr><td>β<sub>n</sub>(V)</td><td>0.125 · e<sup>-V/80</sup></td></tr><tr><td>α<sub>m</sub>(V)</td><td>0.1 · (25 − V) / (e<sup>(25 − V)/10</sup> − 1)</td></tr><tr><td>β<sub>m</sub>(V)</td><td>4 · e<sup>-V/18</sup></td></tr><tr><td>α<sub>h</sub>(V)</td><td>0.07 · e<sup>-V/20</sup></td></tr><tr><td>β<sub>h</sub>(V)</td><td>1 / (e<sup>(30 − V)/10</sup> + 1)</td></tr><tr><td>Stimulus amplitude</td><td>0.0125 mV for Hodgkin–Huxley model, default 0.25 mV for n<sup>4</sup> and p<sub>2</sub> models</td></tr><tr><td>Stimulus frequencies</td><td>0.2, 0.7, 2, 3, 10, 21, 35, 50, 76, 104, 134, 143, 223, 239, 285, 388, 405, 515, 564, 636, 815, 892, 982 (units in Hz)</td></tr><tr><td>Stimulus phases</td><td>Pseudo-random values [0; π] rad</td></tr></table>

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 \underset{\beta_n}{\overset{2\alpha_n}{\rightleftharpoons}} 1 \quad \underset{2\beta_n}{\overset{\alpha_n}{\rightleftharpoons}} 2$$

Thus, the n² model has only two independent rate constants, αₙ and βₙ. 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:

$$p_0 \underset{k_2}{\overset{k_1}{\rightleftharpoons}} p_1 \quad \underset{k_4}{\overset{k_3}{\rightleftharpoons}} p_2 \quad (7)$$

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

More generally, it may be convenient to rewrite k₁, k₂, k₃, k₄ in terms of αₙ and βₙ such that k₁ = Aαₙ, k₃ = αₙ, k₄ = Bβₙ and k₂ = βₙ, 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

$$\left\{ \begin{aligned} \frac{dp_0}{dt} &= -k_1 p_0 + k_2 p_1 \\ \frac{dp_1}{dt} &= k_1 p_0 - (k_2 + k_3) p_1 + k_4 p_2 \\ \frac{dp_2}{dt} &= k_3 p_1 - k_4 p_2 \end{aligned} \right. \quad (8)$$

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

$$\left\{ \begin{aligned} \frac{dp_1}{dt} &= -(k_1 + k_2 + k_3) p_1 + (-k_1 + k_4) p_2 + k_1 \\ \frac{dp_2}{dt} &= k_3 p_1 - k_4 p_2 \end{aligned} \right. \quad (9)$$

This can be rewritten in matrix form

$$\frac{d}{dt} \begin{pmatrix} p_1 \\ p_2 \end{pmatrix} = \begin{pmatrix} \alpha_{11} & \alpha_{12} \\ \alpha_{21} & \alpha_{22} \end{pmatrix} \begin{pmatrix} p_1 \\ p_2 \end{pmatrix} + \begin{pmatrix} -\alpha_{22} - \alpha_{12} \\ 0 \end{pmatrix} \quad (10)$$

where

$$\left\{ \begin{aligned} \alpha_{11} &= -(k_1 + k_2 + k_3) \\ \alpha_{12} &= -k_1 + k_4 \\ \alpha_{21} &= k_3 \\ \alpha_{22} &= -k_4 \end{aligned} \right.$$
At steady state

$$ 
\begin{aligned}
\begin{pmatrix} p_1^\infty \\ p_2^\infty \end{pmatrix} &= \begin{pmatrix} \alpha_{11} & \alpha_{12} \\ \alpha_{21} & \alpha_{22} \end{pmatrix}^{-1} \begin{pmatrix} \alpha_{22} + \alpha_{12} \\ 0 \end{pmatrix} \\
&= \frac{1}{\alpha_{11}\alpha_{22} - \alpha_{12}\alpha_{21}} \begin{pmatrix} \alpha_{22} & -\alpha_{12} \\ -\alpha_{21} & \alpha_{11} \end{pmatrix} \begin{pmatrix} \alpha_{22} + \alpha_{12} \\ 0 \end{pmatrix} \quad (11)
\end{aligned}
$$

$$ p_1^\infty = \frac{\alpha_{22}(\alpha_{22} + \alpha_{12})}{\alpha_{11}\alpha_{22} - \alpha_{12}\alpha_{21}}, \quad p_2^\infty = \frac{-\alpha_{21}(\alpha_{22} + \alpha_{12})}{\alpha_{11}\alpha_{22} - \alpha_{12}\alpha_{21}} \quad (12) $$

## 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

$$ 
\left\{ 
\begin{aligned}
\frac{dp_1}{dt} &= \alpha_{11}p_1 + \alpha_{12}p_2 - (\alpha_{22} + \alpha_{12}), \\
\frac{dp_2}{dt} &= \alpha_{21}p_1 + \alpha_{22}p_2
\end{aligned}
\right.
\qquad (13)
$$

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

$$ 
\left\{ 
\begin{aligned}
\delta \frac{dp_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 \\
&\quad + \alpha_{12} \delta p_2 - \frac{d(\alpha_{22} + \alpha_{12})}{dV} \delta V \\
\delta \frac{dp_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{aligned}
\right.
$$

or equivalently

$$ 
\left\{ 
\begin{aligned}
\delta \frac{dp_1}{dt} - \alpha_{11} \delta p_1 - \alpha_{12} \delta p_2 &= \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] \delta V \\
\delta \frac{dp_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{aligned}
\right.
$$

Applying the Fourier transform with angular frequency ω

$$ 
\left\{ 
\begin{aligned}
(i\omega - \alpha_{11})\hat{\delta p}_1 - \alpha_{12}\hat{\delta p}_2 &= \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] \hat{\delta V} \\
(i\omega - \alpha_{22})\hat{\delta p}_2 - \alpha_{21}\hat{\delta p}_1 &= \left[ \frac{d\alpha_{21}}{dV} p_1^\infty + \frac{d\alpha_{22}}{dV} p_2^\infty \right] \hat{\delta V}
\end{aligned}
\right.
$$

Then $\hat{\delta p}_2$ can be deduced from

$$ 
\left\{ 
\begin{aligned}
\alpha_{21}(i\omega - \alpha_{11})\hat{\delta p}_1 - \alpha_{21}\alpha_{12}\hat{\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]\hat{\delta V} \\
(i\omega - \alpha_{11})(i\omega - \alpha_{22})\hat{\delta p}_2 - (i\omega - \alpha_{11})\alpha_{21}\hat{\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]\hat{\delta V}
\end{aligned}
\right.
$$

By summation

$$ 
\begin{aligned}
&[(i\omega - \alpha_{11})(i\omega - \alpha_{22}) - \alpha_{21}\alpha_{12}]\hat{\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]\hat{\delta V} \\
&\quad + (i\omega - \alpha_{11})\left[\frac{d\alpha_{21}}{dV}p_1^\infty + \frac{d\alpha_{22}}{dV}p_2^\infty\right]\hat{\delta V}
\end{aligned}
$$

This implies the rate of variation

$$ \frac{\hat{\delta p}_2}{\hat{\delta V}} = \frac{\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) + (i\omega - \alpha_{11}) \left( \frac{d\alpha_{21}}{dV} p_1^\infty + \frac{d\alpha_{22}}{dV} p_2^\infty \right)}{(i\omega - \alpha_{11})(i\omega - \alpha_{22}) - \alpha_{21}\alpha_{12}} $$

The current across the lipid bilayer is given by

$$ I = C_m \frac{dV}{dt} + g_L(V - V_L) + g_K p_2(V - V_K) $$

It is linearized by

$$ \delta I = C_m \frac{d(\delta V)}{dt} + g_L \delta V + g_K [\delta p_2(V_0 - V_K) + p_2^\infty \delta V] $$

Applying the Fourier transform with angular frequency ω

$$ \hat{\delta I} = C_m i \omega \hat{\delta V} + g_L \hat{\delta V} + g_K [\hat{\delta p}_2(V_0 - V_K) + p_2^\infty \hat{\delta V}] $$

This provides the admittance of the p₂ model

$$ \hat{Y} = \frac{\hat{\delta I}}{\hat{\delta V}} = C_m i \omega + g_L + g_K \left( \frac{\hat{\delta p}_2}{\hat{\delta V}} (V_0 - V_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( p_1 - p_1^\infty \right) = \left( \begin{matrix} \alpha_{11} & \alpha_{12} \\ \alpha_{21} & \alpha_{22} \end{matrix} \right) \left( \begin{matrix} p_1 - p_1^\infty \\ p_2 - p_2^\infty \end{matrix} \right) $$

More concisely

$$ \frac{dp}{dt} = Mp $$

$$ \text{where } \mathbf{p} = \begin{pmatrix} p_1 - p_1^\infty \\ p_2 - p_2^\infty \end{pmatrix} \text{ and } \mathbf{M} = \begin{pmatrix} \alpha_{11} & \alpha_{12} \\ \alpha_{21} & \alpha_{22} \end{pmatrix}. $$

The general solution is a linear combination

$$ \mathbf{p} = c_1 e^{\lambda_1 t} \mathbf{v}_1 + c_2 e^{\lambda_2 t} \mathbf{v}_2 \quad (14) $$

where $\lambda_1, \lambda_2$ are eigenvalues of $M, v_1, v_2$ are eigenvectors of M and $c_1, c_2$ are constants. Indeed
$$
\begin{aligned}
\mathbf{M}\mathbf{p} &= c_1 e^{\lambda_1 t} \mathbf{M}\mathbf{v}_1 + c_2 e^{\lambda_2 t} \mathbf{M}\mathbf{v}_2 \\
&= c_1 e^{\lambda_1 t} \lambda_1 \mathbf{v}_1 + c_2 e^{\lambda_2 t} \lambda_2 \mathbf{v}_2 \\
&= \frac{dp}{dt}
\end{aligned}
$$

The trace $T = \alpha_{11} + \alpha_{22}$ and determinant $D = \alpha_{11}\alpha_{22} - \alpha_{12}\alpha_{21}$ of $\mathbf{M}$ determine the eigenvalues

$$
\lambda_{1,2} = \frac{T \pm \sqrt{T^2 - 4D}}{2}
$$

$$
\lambda_{1,2} = \frac{\alpha_{11} + \alpha_{22} \pm \sqrt{(\alpha_{11} + \alpha_{22})^2 - 4(\alpha_{11}\alpha_{22} - \alpha_{12}\alpha_{21})}}{2}
$$

Since $p$ represents probabilities, $\lambda_1$ and $\lambda_2$ must be negative to
have exponential decreases rather than exponential increases. This
determines two time constants for the $p_2$ model

$$
\tau_1 = -\frac{1}{\lambda_1}, \quad \tau_2 = -\frac{1}{\lambda_2}
$$

The eigenvectors **v**<sub>1</sub> and **v**<sub>2</sub> are obtained by solving the equation
**Mv** = **λv** and substituting eigenvalues **λ** = **λ**<sub>1</sub> and **λ** = **λ**<sub>2</sub>

$$
\mathbf{v}_1 = \begin{pmatrix} \frac{\lambda_1 - \alpha_{22}}{\alpha_{21}} \\ 1 \end{pmatrix}, \quad \mathbf{v}_2 = \begin{pmatrix} \frac{\lambda_2 - \alpha_{22}}{\alpha_{21}} \\ 1 \end{pmatrix}
$$

Inserting **v**<sub>1</sub> and **v**<sub>2</sub> into Equation 14

$$
\left( \begin{array}{c} p_1 - p_1^\infty \\ p_2 - p_2^\infty \end{array} \right) = c_1 e^{\lambda_1 t} \left( \begin{array}{c} \frac{\lambda_1 - \alpha_{22}}{\alpha_{21}} \\ 1 \end{array} \right) + c_2 e^{\lambda_2 t} \left( \begin{array}{c} \frac{\lambda_2 - \alpha_{22}}{\alpha_{21}} \\ 1 \end{array} \right)
$$

The time $t = 0$ determines the coefficients $c_1$ and $c_2$ as follows

$$
\left\{
\begin{aligned}
p_1(0) - p_1^\infty &= c_1 \frac{\lambda_1 - \alpha_{22}}{\alpha_{21}} + c_2 \frac{\lambda_2 - \alpha_{22}}{\alpha_{21}} \\
p_2(0) - p_2^\infty &= c_1 + c_2
\end{aligned}
\right.
\\[1em]
\left\{
\begin{aligned}
p_1(0) - p_1^\infty &= c_1 \frac{\lambda_1 - \alpha_{22}}{\alpha_{21}} + (p_2(0) - p_2^\infty - c_1) \frac{\lambda_2 - \alpha_{22}}{\alpha_{21}} \\
c_2 &= p_2(0) - p_2^\infty - c_1
\end{aligned}
\right.
\\[1em]
\left\{
\begin{aligned}
c_1 \left( \frac{\lambda_2 - \lambda_1}{\alpha_{21}} \right) &= p_2(0) \frac{\lambda_2 - \alpha_{22}}{\alpha_{21}} - p_2^\infty \frac{\lambda_2 - \alpha_{22}}{\alpha_{21}} - p_1(0) + p_1^\infty \\
c_2 &= p_2(0) - p_2^\infty - c_1
\end{aligned}
\right.
\\[1em]
\left\{
\begin{aligned}
c_1 &= \frac{(\lambda_2 - \alpha_{22})(p_2(0) - p_2^\infty) - \alpha_{21}(p_1(0) - p_1^\infty)}{\lambda_2 - \lambda_1} \\
c_2 &= p_2(0) - p_2^\infty - c_1
\end{aligned}
\right.
$$

The power spectrum can be calculated using the same
reasoning as with the model $n^4$ 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_2(t \mid 0)$ and the steady state probability $p_2^\infty$

$$
C(t) = p_2^\infty \cdot p_2(t \mid 0) - (p_2^\infty)^2
$$

Considering the unique state ($p_0, p_1, p_2$) = (0, 0, 1) when the
channel is open

$$
\begin{cases}
c_1 = \frac{(\lambda_2 - \alpha_{22})(1 - p_2^\infty) + \alpha_{21}p_1^\infty}{\lambda_2 - \lambda_1} \\
c_2 = 1 - p_2^\infty - c_1
\end{cases}
$$

As a result, the conditional probability $p_2(t|0)$ is given by

$$
p_2(t | 0) = p_2^{\infty} + c_1 e^{\lambda_1 t} + c_2 e^{\lambda_2 t}
$$

The autocovariance is given by

$$
\begin{aligned}
C(t) &= p_2^\infty [p_2^\infty + c_1 e^{\lambda_1 t} + c_2 e^{\lambda_2 t}] - (p_2^\infty)^2 \\
&= p_2^\infty c_1 e^{\lambda_1 t} + p_2^\infty c_2 e^{\lambda_2 t}
\end{aligned}
$$

This implies autocovariance $C_{I_K}(t)$ of the total fluctuating
current

$$
\begin{aligned}
C_{I_K}(t) &= N_K I_K^2 C(t) \\
&= N_K I_K^2 p_2^\infty (c_1 e^{\lambda_1 t} + c_2 e^{\lambda_2 t})
\end{aligned}
$$

Following the same conventions as previously, the power
spectrum is

$$
\begin{aligned}
S_{I_K}(\omega) &= 2S_{I_K}^{\pm}(\omega) \\
&= 4\Re \int_0^{+\infty} C_{I_K}(t) e^{-i\omega t} dt \\
&= 4N_K I_K^2 p_2^\infty \Re \int_0^{+\infty} (c_1 e^{\lambda_1 t} + c_2 e^{\lambda_2 t}) e^{-i\omega t} dt
\end{aligned}
$$

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

$$
S_q(\omega) = \int_0^{+\infty} c_q e^{\lambda_q t} e^{-i\omega t} dt
$$

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

$$
S_q(\omega) = c_q \frac{-\lambda_q}{\lambda_q^2 + \omega^2}
$$

Then, the power spectrum is deduced

$$
S_{I_K}(\omega) = 4N_K I_K^2 p_2^\infty \left( c_1 \frac{-\lambda_1}{\lambda_1^2 + \omega^2} + c_2 \frac{-\lambda_2}{\lambda_2^2 + \omega^2} \right)
$$

By using the time constants $\tau_1 = -\frac{1}{\lambda_1}$ and $\tau_2 = -\frac{1}{\lambda_2}$

$$
\begin{aligned}
S_{I_K}(\omega) &= 4N_K I_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_K I_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{aligned}
$$

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 $|\hat{Y}|^2$. 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.

## 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^4$ model, its $n^2$ simplification, or the generalization of $n^2$ to the $p_2$ 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_2$ 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_2$ model with the $n^4$ model. This will be done in the following voltage clamp simulations for the model with nonlinear exponentiation $n^4$ and the model without exponentiation $p_2$ 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^4$.

The control case for nonlinear exponentiation is the Hodgkin-Huxley $n^4$ 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^4$ 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^4$ 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 $\omega_i, \omega_j$ represents an interactive quadratic frequency $\omega_i + \omega_j$ for $i,j \in \{-N, \dots, -1, +1, \dots, +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(|\omega_i| + |\omega_j|)$ at frequency sums (red points), $S_M(|\omega_i| - |\omega_j|)$ at frequency differences (blue points), $S_D(2|\omega_k|)$ at frequency doubling (orange points) and $S_R(|\omega_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|\omega_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.

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
<figure><img src="image_4.png" /><figcaption>FIGURE 3</figcaption></figure>

Markov simulations of the $p_2$ 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_1$ (brown curve) and $S_2$ (orange curve). The corresponding corner frequencies $\omega_1$ and $\omega_2$ 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 $|\hat{Y}|^2$ (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^4$ model, which can be compared to Figure 2C.

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 $|\omega_i| + |\omega_j|$, $||\omega_i| - |\omega_j||$ and $2|\omega_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(|\omega_i| + |\omega_j|)$ and $S_M(||\omega_i| - |\omega_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|\omega_k|)$ that approximately follow a smooth line. This is because frequency doubling is calculated along the frequency diagonal $(\omega_k, \omega_k)$ on the
<figure><img src="image_5.png" /><figcaption>FIGURE 4</figcaption></figure>

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 (ωᵢ, ωⱼ, |Qᵢⱼ|) 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<sub>P</sub>(|ωᵢ| + |ωⱼ|) at frequency sums (red points), S<sub>M</sub>(||ωᵢ| - |ωⱼ||) at frequency differences (blue points), S<sub>D</sub>(2|ωₖ|) at frequency doubling (orange points) and S<sub>R</sub>(|ωⱼ|) 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.

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<sub>R</sub>(|ω<sub>j</sub>|) 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<sub>D</sub>(2|ω<sub>k</sub>|) tends to zero at high frequencies due to the attenuation of quadratic responses along the diagonal (ω<sub>k</sub>, ω<sub>k</sub>), while the power spectra S<sub>R</sub>(|ω<sub>j</sub>|) 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<sub>P</sub>(|ω<sub>i</sub>| + |ω<sub>j</sub>|) and S<sub>M</sub>(||ω<sub>i</sub>| - |ω<sub>j</sub>|) are fundamentally due to non-trivial quadratic frequency interactions between ω<sub>i</sub> and ω<sub>j</sub> 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<sub>IK</sub> represents the analytical estimate of Markov fluctuations. The curve S<sub>L</sub> represents the linear power spectrum at fundamental frequencies. The curve S<sub>D</sub> represents the quadratic power spectrum at frequency doubling.

The curve S<sub>R</sub> represents the mean squared quadratic output by matrix columns. The multi-sinusoidal power spectra S<sub>L</sub>, S<sub>D</sub>, S<sub>R</sub> were averaged over different random sets of non-overlapping frequencies as in Figure 4. The curve |ŷ<sub>m</sub>|² represents the squared admittance modified from the Equation 1 without the membrane capacitance nor the frequency independent terms, namely

$$ \hat{Y}_m = g_K [4n_0^3 (V_0 - V_K) \hat{D}_n] $$

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 |ŷ<sub>m</sub>|² is more marked than that of the Markov fluctuations S<sub>IK</sub>, which is also the case for the linear power spectrum S<sub>L</sub> over a limited range of frequencies before the reversal. The fall of the quadratic power spectra S<sub>D</sub> and S<sub>R</sub> is more marked than that of the Markov fluctuations S<sub>IK</sub> but less than that of the modified squared
<figure><img src="image_6.png" /></figure>

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^4$ in the stochastic equations.

Figure 6 shows the frequency analysis of the $n^4$ 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 5B (55 mV) shows the superposition of five curves of the $n^4$ model for a 55 mV depolarization, using the same presentation as in Figure 5A (5 mV). The modified squared admittance $|\hat{Y}_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_1$ fits the Markov fluctuations $S_{IK}$ accurately.

Thus, at 55 mV depolarization, the fluctuation behavior $S_{IK}$ is dominated by $S_1$ and has a kinetic behavior similar to that of a modified squared admittance $|\hat{Y}_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 $|\hat{Y}_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^4$ 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.
<figure><img src="image_7.png" /><figcaption>FIGURE 6
Frequency analysis of the $n^4$ 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$).</figcaption></figure>

## 3.6 Observations on the linear fluctuation-dissipation theorem

As discussed by Stevens (1972), the power $n^4$ 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_K[4n_0^3(V_0 - V_K)\bar{D}_n + n_0^4]$ 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^4$ 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_2$ model described in this article.

## 3.7 Comparison between $n^4$ and $p_2$ models

### 3.7.1 Depolarization of 55 mV

As discussed above, sequential kinetic models can produce the exponentiation of the $n^4$ Hodgkin-Huxley model, so that the

probability of the channel opening coincides with $n^4$. Similarly, for the $n^2$ model, the probability of the channel being open coincides with $n^2$. However, the $p_2$ model generalizes the rate constants of the $n^2$ 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_2$ model, in which the rate constants have been manually selected to approximate the frequency domain responses of the $n^4$ 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_2$ model (Figure 7) show not one, but two resonances reflected by a peak and bump, unlike the $n^4$ model (Figure 6). Second, the Markov fluctuations of the $p_2$ model (Figure 3C) are also different with the indication of more than one time constant illustrated by
<figure><img src="image_8.png" /><figcaption>FIGURE 7</figcaption></figure>

Frequency analysis of the $p_2$ model for a 55 mV depolarization. The rate constants have been manually selected to approximate the frequency domain responses of the $n^4$ 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$).

the inflection just before the final slope of $S_2$, unlike the $n^4$ model (Figure 2C) for which $S_1$ fits $S_{IK}$ accurately.

Figure 8B shows the superposition of five curves of the $p_2$ model for a 55 mV depolarization, using the same presentation as for the $n^4$ model in Figure 5B. As with the $n^4$ model, the modified squared admittance $|\hat{Y}_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 $|\hat{Y}_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_2$ model. In contrast, the $n^4$ model has only one relaxation process due to the linearization of the exponentiation. Similarly, the modified squared admittance $|\hat{Y}_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.

## 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^4$ 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|\omega_k|)$ at frequency doubling for which the amplitudes decrease at high frequencies.

Figure 8A shows that, as with the $n^4$ model at 5 mV in Figure 5A, the linear behavior represented by the modified squared admittance $|\hat{Y}_m|^2$ does not accurately describe quadratic power spectra or Markov fluctuations. In particular, as with the $n^4$ model at 5 mV, the modified squared admittance and quadratic power spectra of the $p_2$ 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^4$ model and the $p_2$ 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.
<figure><img src="image_9.png" /></figure>

# 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_2$ 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^4$ 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_2$ 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 $|\hat{Y}_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 11 supports this hypothesis as illustrated with the power spectra of the QSA matrix coefficients of the $p_2$ 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 $|\hat{Y}_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 $|\hat{Y}_m|^2$ have a similar appearance.

Figures 12–15 reinforce this interpretation by varying the surface area of the membrane of the $p_2$ 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_K = 50\mu m^2$, $A_K = 500\mu m^2$, $A_K = 5,000\mu 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_K = 50,000\mu 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.

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
<figure><img src="image_10.png" /><figcaption>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).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>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).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIGURE 11</figcaption></figure>

Power spectra of the QSA matrix coefficients of the $p_2$ 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).

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_2$ model with existing approaches

The $p_2$ model is constructed as a variant of the $n^4$ 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_2$ model avoids explicit exponentiation, yet it can represent the minimal degree of nonlinearity $n^2$ as well

as models without exponentiation. Moreover, the rate constants of the $p_2$ model can be tuned to preserve the key dynamical features of the $n^4$ model, as evaluated through spectral analysis methods.

The relevance of decreasing the degree of exponentiation to $n^2$ 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_2$ model generalizes this approach by providing a flexible three-state Markov representation with independent rate constants.

The master equation of the $p_2$ 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_2$ 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.
<figure><img src="image_13.png" /><figcaption>FIGURE 12</figcaption></figure>

Power spectra of the QSA matrix coefficients of the $p_2$ model for a membrane surface area $A_K = 50\mu 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><img src="image_14.png" /><figcaption>FIGURE 13</figcaption></figure>

Power spectra of the QSA matrix coefficients of the $p_2$ model for a membrane surface area $A_K = 500\mu 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><img src="image_15.png" /></figure>

FIGURE 14

Power spectra of the QSA matrix coefficients of the *p*₂ model for a membrane surface area *A*<sub>K</sub> = 5,000µm², a 55 mV depolarization and a stimulus amplitude 1 mV. The plots were generated using the same presentation as in Figure 12.

<figure><img src="image_16.png" /></figure>

FIGURE 15

Power spectra of the QSA matrix coefficients of the *p*₂ model for a membrane surface area *A*<sub>K</sub> = 50,000µm², a 55 mV depolarization and a stimulus amplitude 1 mV. The plots were generated using the same presentation as in Figure 12.
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^4$ and $p_2$ 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 (NaV1.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^4$ 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^4$. 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_2$ model does not rely on the concept of independent gates. Instead,

it relies on independent rate constants. This approach allows the $p_2$ model to represent behaviors such as $n^2$, as well as variants without explicit exponentiation. Moreover, this article has demonstrated that the $p_2$ and $n^4$ models can produce similar neuronal responses, suggesting that the underlying kinetic schemes are more essential than gating independence for capturing channel behavior.

The $p_2$ model (or an extension) could potentially be further adapted to approximate experimental data as demonstrated here theoretically to approximate the $n^4$ 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_2$ 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{dn}{dV_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.

## References

Andreozzi, E., Carannante, I., D'Addio, G., Cesarelli, M., and Balpi, B. (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

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

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

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

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

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

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

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

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

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

## 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.

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

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

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

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

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

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

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

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

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/BF019 05222

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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.304592

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