Frequency Switches at Transition Temperature in Voltage-Gated Ion Channel Dynamics of Neural Oscillators^*

Phase transitions are important physical phenomena. A typical example is the transition between the ferromagnetic and paramagnetic phases of magnetic materials at the Curie point. In this study, we consider phase transitions in neural information processing in the brain, whereby neural activity transits to another neural activity in the network. Understanding of neural phase transition mechanisms is important for clarifying cognitive functions, such as perception and learning, in humans and higher animals, since cognitive functions largely depend on transitions of neural activity in the brain.^[1–3]

To reveal mechanisms of cognitive functions, we must investigate phase transitions found in single neurons at the membrane potential level. Mathematical analysis of the neural network has shown that the transition from slow wave sleep to rapid eye movement sleep is induced by the switching of stabilities of equilibria established in the saddle node on an invariant circle (SNIC) bifurcation.^[4] Similarly, a physiological experiment has found that transient neural activities, such as membrane potential fluctuations in single neurons, induce phase transitions in the neural networks with several synchronous states.^[5]

A mathematical aspect of phase transitions in single neuron models appears to be less clearly understood, although the physiological aspect has been well characterized. Chiu et al.^[6] found transition temperatures (i.e., breaks in the Arrhenius plot) for sodium conductance and the time constant of sodium inactivation in rabbits. Schwarz^[7] demonstrated that various kinetic parameters of sodium permeability in frog nerve and muscle fibers undergo temperature-induced phase transitions at low temperature. These results are consistent with the hypothesis that temperature-dependent phase transitions found in the membrane lipids.

At the phase transition temperature, single-ion channels in the artificial membranes of unmodified lipid bilayers undergo sudden and dramatic increases in ionic permeability.^[8] This finding implies that ion-specific channel phase transitions, rather than lipid phase transitions, are responsible for the sudden and dramatic increase in ionic permeability. Thus phase transitions in the voltage-gated ion channels of single neurons require further investigation and clarification.

In this Letter, we reveal that both nonlinear and thermodynamic mechanisms affect the voltage-gated ion channel dynamics in the Morris–Lecar (ML) model.^[9] These mechanisms are responsible for two types of switches of frequency at different thermodynamic temperatures (Fig. 1(a)). Specifically, the ML neuron model undergoes different transitions of the bifurcation from the big saddle-homoclinic orbit (BSHCO) type to the SNIC and from the SNIC type to small SHCO (SSHCO) type. In particular, we analyze the temperature-dependent frequency immediately after a saddle-node (SN) equilibrium has disappeared with an increase in the input current (Fig. 1(b)). Let us note that in the BSHCO and SSHCO bifurcations, the SN equilibrium exists inside and outside the limit cycle (LC) attractor, respectively, before disappearing.^[10]

Zoom In Zoom Out Reset image size

We now describe the ML model, defined as

$$\begin{eqnarray}\begin{array}{lll}\frac{dx}{dt} & = & F(x)=({F}_{V}(x),{F}_{W}(x))\\ & = & \left(\begin{array}{c}\frac{1}{{C}_{{\rm{m}}}}({I}_{{\rm{ion}}}(x)+{I}_{0})\\ \frac{\mu }{{\tau }_{{\rm{w}}}(V)}({w}_{\infty }(V)-W)\end{array}\right),\end{array}\end{eqnarray} \tag{ 1 }$$

where I_ion(V, W) = −g_L (V − V_L) − g_KW(V − V_K) − g_Cam_∞(V) (V − V_Ca), m_∞(V) = (1 + tanh [(V − V₁)/V₂])/2, w_∞(V) = (1 + tanh [(V − V₃)/V₄])/2, τ_w(V) = 1/(ϕcosh[(V − V₃)/(2V₄)]), x ∈ (V, W), V (mV) and W represent the membrane potential and the gating variable of the K⁺ channel, respectively, C_m is the membrane capacity, I₀ is the constant input current, V_Ca, V_K and V_L denote the equilibrium potentials of Ca²⁺, K⁺ and the leak currents, respectively, g_Ca, g_K and g_L are the maximum conductances of the corresponding ionic currents, and m_∞(V) denotes the calcium channel permeability P_Ca ∝ m_∞(V).^[11] In transition state theory,^[12] the V-term in the hyperbolic tangent function of m_∞(V) represents the potential energy difference ΔH_m(V) when the membrane potential is V. More precisely, m_∞(V) = [exp((V − V₁)/V₂)/(exp((V − V₁)/V₂) + exp(−(V − V₁)/V₂))] = [k₊₂(V)/(k₊₂(V) + k₋₂(V))], where k₊₂ and k₋₂ are the reaction coefficients of the closed-to-open and open-to-closed transition states, respectively, V₁ is the potential energy difference between the open (or closed) state and the transition state at V = 0 in the Ca²⁺ channel kinetics, V₂ is an assumed constant, although in a thermodynamic interpretation at temperature T, it represents the Boltzmann constant (as explained in detail in the following). A similar thermodynamic interpretation is applied to w_∞(V) and τ_w(V). For τ_w(V), and ϕ implicitly involves the thermodynamic variables of temperature and entropy, although we here set ϕ as a constant. The potential energy difference in the K⁺ channel is denoted by ΔH_w(V), where V₃ is the potential energy difference between the open (or closed) state and the transition state at V = 0. These kinetics depend on the temperature through the relationship $\mu ={Q}_{10}^{(T-{T}_{{\rm{c}}})/10}$ in Eq. (1) where Q₁₀ = 3,^[13] and T is the temperature influencing the voltage-gated channel dynamics. Let T_c = 0^○C be the environmental temperature of the physiological experiment. Thus we assume that the transition states of the channel dynamics vary under constant pressure and volume conditions at T = T_c, according to $[{\tau }_{{\rm{w}}}(T+10)/{\tau }_{{\rm{w}}}(T)]={Q}_{10}^{(T-{T}_{{\rm{c}}})/10}$.

In this work, the parameters are set as C_m = 20 μF/cm², V_K = −80 mV, V_L = −60 mV, V_Ca = 120 mV, g_L = 2 μS/cm², g_K = 8 μS/cm², g_Ca = 4 μS/cm², V₁ = −1.2, V₂ = 18, V₃ = 12 and V₄ = 17.4.^[9] Such a parameter setting establishes a nonlinear dynamical mechanism in which the SNIC (at T = 10^○C), BSHCO (at T = 3^○C) and SSHCO (at T = 26.4^○C) bifurcations occur by an upward shift of the cubic V-nullcline as I₀ (μA/cm²) increases (Fig. 1(b)). We observe that as I₀ passes a critical value I_sn ~ 39.6935 μA/cm², a transition occurs from quiescence to repetitive firing. We also observe two different types of switches of the frequency at different temperatures. Around T = 5^○C, the frequency drastically decreases to its low or 0. At T = 26.4^○C, it undergoes a sudden and drastic increase (Fig. 1(a)).

To analyze how each bifurcation transition affects its relevant frequency switch, we compute the first derivatives with respect to T, which relate to the phase response curve (PRC). The PRC is widely used in neuroscience to characterize firing dynamics in the neuron. Introducing a small perturbation G(x) = (0, γμF_W(x)ΔT) ( ≪ 1) into Eq. (1), we obtain

$$\begin{eqnarray}&&\frac{dx}{dt}=F(x)+\epsilon G(x),\end{eqnarray} \tag{ 2 }$$

where γ = [log(Q₁₀)/10]. We now define a small perturbation in the phase direction on the periodic orbit x_p(t) as q(t_s), where t_s = t is a slow time, x_p(t) is a stable T_p-periodic solution for the unperturbed dynamics, [dx_p/dt] = F(x_p) with x_p (t + T_p) = x_p(t) for all t. Let us consider the vector Z(t) = (Z_V(t), Z_W(t)) called the adjoint solution to the following linearization around the T_p-periodic solution x_p(t): [dZ/dt] = −[D_x F(xp)]^trZ, where D_xF represents the Jacobian matrix of F, and the superscript tr indicates the transpose of a matrix, and Z, frequently regarded as the PRC, is assumed to be normalized by infinitesimal pulsed perturbations.^[14] However, in this study, we denote the PRC by Z_V for convenience. We thus find a differential equation that governs the time evolution of q. Averaging the differential equation over the period T_p, we obtain the following evolution equation for q

$$\begin{eqnarray}\begin{array}{lll}\frac{dq}{d{t}_{{\rm{s}}}} & = & \frac{\epsilon }{{T}_{{\rm{p}}}}\displaystyle {\int }_{0}^{{T}_{{\rm{p}}}}Z({t}^{^{\prime} }+q({t}_{{\rm{s}}}))\\ & & \cdot G({x}_{{\rm{p}}}({t}^{^{\prime} }+q({t}_{{\rm{s}}})))d{t}^{^{\prime} }.\end{array}\end{eqnarray} \tag{ 3 }$$

Introducing a phase variable θ = [(t + q(t_s))/T_p] with θ ∈ [0, 1) and defining the new firing rate as [dθ/dt] = f(T₀ + ΔT),^[14] Eq. (3) becomes

$$\begin{eqnarray}\begin{array}{lll}\frac{d\theta }{dt} & = & \frac{1}{{T}_{{\rm{p}}}}+\frac{\epsilon }{{T}_{{\rm{p}}}}\displaystyle {\int }_{0}^{1}\tilde{Z}({\theta }^{^{\prime} })\cdot G({\tilde{x}}_{{\rm{p}}}({\theta }^{^{\prime} }))d{\theta }^{^{\prime} }\\ & = & f+\gamma \mu \langle {Z}_{W}{F}_{W}\rangle \Delta T,\end{array}\end{eqnarray} \tag{ 4 }$$

where $\tilde{Z}(\theta )=Z({T}_{{\rm{p}}}\theta )$, ${\tilde{x}}_{{\rm{p}}}(\theta )={x}_{{\rm{p}}}({T}_{{\rm{p}}}\theta )$, f = f(T₀) = [1/T_p], and $\langle {Z}_{W}{F}_{W}\rangle =[1/{T}_{{\rm{p}}}]\displaystyle {\int }_{0}^{1}{Z}_{W}{F}_{W}d{\theta }^{^{\prime} }$. The new frequency is straightforwardly expanded in ΔT as follows:

$$\begin{eqnarray}&&\frac{d\theta }{dt}=f({T}_{0}+\Delta T)\simeq f({T}_{0})+\frac{df({T}_{0})}{dT}\Delta T.\end{eqnarray} \tag{ 5 }$$

From Eqs. (4) and (5), we obtain the relationships

$$\begin{eqnarray}\begin{array}{lll}\frac{df}{dT} & = & \gamma \langle {Z}_{W}\frac{dW}{dt}\rangle \\ & = & \gamma f(1-\frac{1}{{C}_{{\rm{m}}}}({\langle {Z}_{V}{I}_{{\rm{ion}}}\rangle }_{\theta }+{\langle {Z}_{V}\rangle }_{\theta }{I}_{0})),\end{array}\end{eqnarray} \tag{ 6 }$$

where ${\langle {Z}_{V}{I}_{{\rm{ion}}}\rangle }_{\theta }=\displaystyle {\int }_{0}^{1}{Z}_{V}{I}_{{\rm{ion}}}d{\theta }^{^{\prime} }$ and ${\langle {Z}_{V}\rangle }_{\theta }=\displaystyle {\int }_{0}^{1}{Z}_{V}d{\theta }^{^{\prime} }$. In Eq. (6), we used Z(t) · F(x_p(t)) = 1 for any t according to Ref. [14]. Therefore, we have found that the derivative of the firing rate with respect to temperature is represented as a function relating to the PRC.

First, let us calculate γ⟨Z_W[dW/dt]⟩ (corresponding to [df/dT]), where ||[df/dT] − γ⟨Z_W[dW/dt]⟩|| = O(10⁻¹). This calculation elucidates the effects of the PRC on the slope of the firing rate-temperature curve. In Fig. 2(a), the SN equilibrium remains at I₀ = 39.693. The slope of the frequency with respect to T (γ⟨Z_W [dW/dt]⟩) diverges to the negative infinity as T approaches the transition temperature of the BSHCO to SNIC, after which it discontinuously resets to 0. Meanwhile, in the SNIC to SSHCO bifurcation transition, the slope of the frequency with respect to T discontinuously changes from 0 to some positive value. Following the disappearance of the SN equilibrium, γ⟨Z_W[dW/dt]⟩ (> 0 around T = 0^○C) develops a negative valley at the BSHCO to SNIC bifurcation transition and returns to positive in the temperature range of the SSHCO bifurcation. These results indicate a conspicuous difference between dynamical properties at the bifurcation transitions from BSHCO to SNIC and from SNIC to SSHCO.

Zoom In Zoom Out Reset image size

We now show how this conspicuous difference arises by evaluating Eq. (6) (Fig. 2(b)). In the high temperature domain (T > 26.4^○C) of the SSHCO bifurcation, ⟨Z_W[dW/dt]⟩_θ > 0, because [⟨Z_V⟩_θI₀/C_m] > [⟨Z_V I_ion⟩_θ/C_m]. By contrast, in the temperature domain of the SNIC and BSHCO bifurcations, [⟨Z_V⟩_θI₀/C_m] < [⟨Z_V I_ion⟩_θ/C_m] and ⟨Z_W[dW/dt]⟩_θ can be negative. In fact, the PRC in the SSHCO and SNIC bifurcations develop different shapes, as reported in Ref. [15]. Furthermore, the PRC in SSHCO is larger by O(10¹) than in the SNIC bifurcation. Such differences between the PRC shapes and sizes in the SNIC and SSHCO bifurcations cause γ⟨Z_W[dW/dt]⟩ to switch from positive values to negative values and vice versa.

For a thermodynamic interpretation of the bifurcation transitions, let us compute τ_w(V) and m_∞(V). For τ_w(V), we observe two breaks in the Arrhenius plot of ${\langle {\tau }_{{\rm{w}}}\rangle }_{\theta }$ ($=\displaystyle {\int }_{0}^{1}{\tau }_{{\rm{w}}}(V({\theta }^{^{\prime} }))d{\theta }^{^{\prime} }$), occurring at different transition temperatures (Fig. 3(a)). However, in the ⟨τ_w⟩ plot, no break appeared (data not shown). This indicates that the size of τ_w, not the average τ_w over one period, is responsible for obtaining two different discontinuities in the Arrhenius plot.

Zoom In Zoom Out Reset image size

We also plotted the potential energies ΔH_w versus temperature. Along the horizontal axis of Fig. 3(b), the ⟨ΔH_w⟩−T curves are almost symmetric to the f−T curves. At I₀ = 39.693, the potential energy reaches the thermodynamic equilibrium through the first bifurcation transition, followed by a drastic decline as T increases further. At each bifurcation transition, we have identified different discontinuities in the first derivative of the potential energy with respect to temperature. At I₀ = 39.693, [d⟨ΔH_w⟩/dT] switches from +∞ to 0 through the bifurcation transition from BSHCO to SNIC. Furthermore, it rapidly ascends from −∞ following the second bifurcation transition, i.e., from SNIC to SSHCO.

To gain further insights, we now numerically calculate ⟨m_∞⟩. At I₀ = 39.693, as T increases, ⟨m_∞⟩ is reset to 0 at the BSHCO to SNIC bifurcation transition, followed by a rapid positive increase following the SNIC to SSHCO bifurcation transition. The shape of this curve is similar to that of the f–T relation (Fig. 3(d)), and reveals breaks in permeability at the phase transition temperature of the corresponding bifurcation transition. The first derivative of ⟨m_∞⟩ with respect to T reveals different discontinuities at different bifurcation transitions, similar to those of Fig. 2(a). To verify that these discontinuities represent first-order phase transitions, we calculated the average of the potential energy term in m_∞, ⟨ΔH_m⟩ (results not shown). Similar to Figs. 3(b) and 3(c), [d⟨ΔH_m⟩/dT] developed discontinuities at different bifurcation transitions. Therefore, we infer that these two frequency switches, at which different bifurcation transitions occur, are the first-order phase transition temperatures.

We have largely neglected the important roles of bifurcation transitions in characterizing neuronal activity, because various types of bifurcations have been understood by analyzing the stabilities of the equilibria in the system. For instance, in a modified version of the Hodgkin–Hulxey (HH) model, a SHCO bifurcation of the SN equilibrium was found to be embedded in a chaotic attractor.^[16] Nevertheless, it remains unclear why the SHCO bifurcation causes the intervals between spikes to expand. Although the original HH model is employed for schematically analyzing different types of bifurcations, its intrinsic parameters^[17,18] do not have to be subjected to a rigorous thermodynamic interpretation. In this work, we have suggested that the two different bifurcation transitions of the BSHCO to SNIC and SNIC to SSHCO are related to first-order phase transition phenomena at the level of the voltage-gated channel dynamics.

As mentioned above, the equilibrium state of the calcium channel in Eq. (1) largely depends on the ion permeability of the membrane. Figure 3(d) displays sudden and drastic decrease and increase of the permeability at T = 5^○C and 26.4^○C, respectively. The drastic increase in the ion permeability is consistent with a previously reported result, which relates the appearance of a single-ion channel to the dramatic increase in its relevant ion permeability at the phase transition temperature of unmodified lipid bilayer membranes.^[8] To account for the additional phase transition in the membrane lipids, a nonlinear term expressing the lipid domain interaction can be inserted into the potential energy term of the transition state. This term is expected to reasonably explain the hysteresis caused by lipid phase transitions. However, in the present study, we have found that nonlinear dynamical characteristics, particularly bifurcation transitions, significantly influence the potential energies of the voltage-gated channels in the ML model, regardless of lipid phase transitions.

In conclusion, we have confirmed that frequency switches found in the ML system, where the bifurcation transitions occur from BSHCO to SNIC and from SNIC to SSHCO as the temperature increases, represent different first-order phase transition phenomena. The nonlinear and thermodynamic characteristics of the two phase transitions are markedly different. In the BSHCO to SNIC bifurcation transition, the potential energy reaches thermodynamic equilibrium (ΔH_w = 0 or ΔH_m = 0) as the temperature increases. The first derivatives of the potential energy diverge towards +∞ immediately before this transition. Regarding the nonlinear dynamics, the slopes of the firing rate with respect to temperature can predominantly be composed of ⟨Z_V I_ion⟩_θ rather than ⟨Z_V⟩_θ. On the other hand, in the SNIC to SSHCO bifurcation transition, both the potential energies and their first derivatives drop towards −∞ immediately following the transition. Interestingly, ⟨Z_V⟩_θ, rather than ⟨Z_V I_ion⟩_θ, determines [df/dT]. Thus we have identified two first-order phase transition temperatures, at which transitions are induced by potential energy changes in biophysical channel systems, independent of lipid phase transitions.