Bistability at the onset of neuronal oscillations

Abstract

The Hodgkin–Huxley (HH) model and squid axon (bathed in reduced Ca²⁺) fire repetitively for steady current injection. Moreover, for a current-range just suprathreshold, repetitive firing coexists with a stable steady state. Neuronal excitability, as such, shows bistability and hysteresis providing the opportunity for the system to perform as switchable between firing and non-firing states with transient input and providing the backbone as a dynamical mechanism for bursting oscillations. Some conditions for bistability can be derived by intricate analysis (bifurcation theory) and characterized by simulation, but conditions for emergence and robustness of such bistability do not typically follow from intuition. Here, we demonstrate with a semi-quantitative two-variable, V–w, reduction of the HH model features that promote/reduce bistability. Visualization of flow and trajectories in the V–w phase plane provides an intuitive grasp for bistability. The geometry of action potential recovery involves a late phase during which the dynamic negative feedback of ${I}_{\rm Na}$ inactivation and ${I}_{\rm K}$ activation over/undershoot, respectively, their resting values, thereby leading to hyperexcitabilty and an intrinsically generated opportunity to by-pass the spiral-like stable rest state and initiate the next spike upstroke. We illustrate control of bistability and dependence of the degree of hysteresis on recovery timescales and gating properties. Our dynamical dissection reveals the strongly attracting depolarized phase of the spike, enabling approximations like the resetting feature of adapting integrate-and-fire models. We extend our insights and show that the Morris–Lecar model can also exhibit robust bistability.

Appendices

Appendix A: Numerical methods

1.1 Time courses

The time courses and limit cycles for HH, V-w and ML model were solved by the 4^th-order Runge–Kutta method with step size $dt=0.01$ ms in MATLAB, except that the unstable limit cycles in Fig. 2B, D and 3B, C were obtained by grabbing initial data on an unstable limit cycle from a one-parameter bifurcation diagram in XPP AUTO and running backward in time. Simulations were run for $T=$ 200 ms with step size $dt=$ 0.01 ms. Trajectories were plotted starting from 100 ms after the initial transients decay. The nullclines and fixed points were solved by MATLAB build-in nonlinear solver ‘lsqnonlin’ with the default stopping criteria according to Eq. (1) for the V-w model. The version of MATLAB is R2021b.

1.2 Phase planes

In most cases, we added the heat map or flow vectors onto the phase plane to show the flow more clearly. When plotting the flow vectors in the phase plane, we divided $V$ by a factor of 100 due to the scaling difference in membrane potential $V$ and the ${I}_{\rm K}$ gating variable, $w$. The factor was chosen based on the amplitude of $V$ (around 100 mV). The white arrows for the flow represent the normalization of the vector$(dV/(100\cdot dt) , dw/dt)$. So they only show the directions of the flow but carry no information about the magnitudes. The heatmap is a continuous visualization of the directions of the flows. We calculated the angle of the scaled flow (the white arrow) with the closer horizontal direction (either $(\mathrm{1,0})$ or$(-\mathrm{1,0})$) on a finer grid and displayed it by color, where yellow represents 90 degrees (vertical) and blue represents 0 degree (horizontal). By the nature of nullclines, it is expected that the heatmap is yellow near the $V$-nullcline, while blue near the w-nullcline.

1.3 Bifurcation diagrams

One-parameter, two-parameter bifurcation diagrams and frequency diagrams were usually computed using AUTO as incorporated into XPP. The one-parameter bifurcation diagrams and the frequency plots calculated by XPP AUTO were also reproduced from MATLAB using Ting-Hao Hsu’s MATLAB Interface for XPP produced diagrams (available on the webpage for XPP: http://www.math.pitt.edu/~bard/xpp/xpp.html).

In order to calculate DoB, we also exported the data of two-parameter bifurcation diagrams from XPP AUTO and interpolated it in MATLAB so that the sample points in the second parameter were consistent for LP and HB.

Due to the difficulty in finding appropriate numerical parameters to track some bifurcation points in AUTO, some of the ${I}_{LP}$ curves in two-parameter bifurcation diagrams (such as Fig. 5Bb) were calculated in MATLAB by a different method and combined with AUTO outputs. We drove the V-w model with a slow linearly ramping down ${I}_{\rm app}$ and identified ${I}_{LP}$ by the point where oscillation terminated. The simulation time duration was set long enough (10,000 ms) to attain a reasonable resolution.

Appendix B: Measure for degree of bistability (DoB)

The definition of the measure degree of bistability (DoB) varies from figure to figure depending on the parameter of interest. The general rule is that for a parameter the value of DoB is given by the length of bistable interval with respect to the parameter normalized by its mean. For each one-parameter bifurcation diagram, the measure was computed and reported in the figure caption. For each two-parameter bifurcation diagram, the dependence of the measure on a second parameter was visualized and usually included as a plot.

To be more specific, in the ${I}_{\rm app}$ bifurcation diagram for the HH model and V-w models (Fig. 2E), the system is bistable in the interval [${I}_{LP}$, ${I}_{HB1}$]. Then the degree of bistability for ${I}_{\rm app}$ is defined as the length of this interval over its mean value, which is $Do{B}_{Iapp}=\frac{{I}_{HB1}-{I}_{LP}}{\left({I}_{HB1}+{I}_{LP}\right)/2}$. Furthermore, for a two-parameter bifurcation diagram like Fig. 2F, where a second parameter temperature is considered, we calculated the value of DoB for each fixed temperature and the corresponding ${I}_{LP}$ and ${I}_{HB1}$ values. Similarly, Fig. 5Ac is a $\overline{{g}_{\rm K}}$ bifurcation diagram, which shows the system is bistable in the interval [${\overline{{g}_{\rm K}}}_{HB2}$, ${\overline{{g}_{\rm K}}}_{LP}$]. Following the same rule, $Do{B}_{\overline{{g}_{\rm K}}}=\frac{{\overline{{g}_{\rm K}}}_{LP}-{\overline{{g}_{\rm K}}}_{HB2}}{\left({\overline{{g}_{\rm K}}}_{LP}+{\overline{{g}_{\rm K}}}_{HB2}\right)/2}$. In Fig. 5Bb, ${I}_{\rm app}$ is the second parameter of interest, so DoB for $\overline{{g}_{\rm K}}$ was plotted as a function of ${I}_{\rm app}$ in Fig. 5Bc.

We also considered other measures, such as the range of bistability over the range of oscillation. However, since the termination of oscillation is not necessarily via a Hopf bifurcation, the measure does not reflect the relative range of bistability well.

Appendix C: Hodgkin–Huxley model and the validity of the reduction to V-w model

3.1 Hodgkin–Huxley model

We considered the four-variable Hodgkin–Huxley model (Hodgkin and Huxley 1952):

$${\begin{array}{l}{C}_{M}\frac{dV}{dt}={I}_{\rm app}-{I}_{\rm Na}-{I}_{\rm K}-{I}_{L}={I}_{\rm app}-{\overline{g}}_{\rm Na}{m}^{3}h \left(V-{E}_{\rm Na}\right)\\ -{\overline{g} }_{\rm K}{n}^{4}\left(V-{E}_{\rm K}\right)-{\overline{g}}_{L}\left(V-{E}_{L}\right)\\ \frac{dm}{dt}=\phi \cdot \frac{{m}_{\infty }\left(V\right)-m}{{\tau }_{m}\left(V\right)}\\\frac{dh}{dt}=\phi \cdot \frac{{h}_{\infty }\left(V\right)-h}{{\tau }_{h}\left(V\right)}\\ \frac{dn}{dt}=\phi \cdot \frac{{n}_{\infty }\left(V\right)-n}{{\tau }_{n}\left(V\right)}\end{array}}$$

(2)

where ${I}_{\rm app}$ is the applied current; ${I}_{\rm Na}=-{\overline{g} }_{\rm Na}{m}^{3}h\left(V-{E}_{\rm Na}\right)$, ${I}_{\rm K}=-{\overline{g} }_{\rm K}{n}^{4}\left(V-{E}_{\rm K}\right)$ and ${I}_{L}=-{g}_{L}\left(V-{E}_{L}\right)$ are the ionic currents. The maximal conductance for Na⁺, K⁺ and leak currents is

${\overline{g} }_{\rm Na}=120 mS/{\rm cm}^{2}$, ${\overline{g} }_{\rm K}=36 mS/{\rm cm}^{2}$ and ${g}_{L}=0.3 mS/{\rm cm}^{2}$, respectively. The reversal potentials are ${V}_{\rm Na}=50mV$, ${V}_{\rm Na}=-77mV$, ${V}_{L}=-54.4mV$. The membrane capacitance is ${C}_{m}=1 \mu F/{\rm cm}^{2}$. The correction factor for temperature, $\mathrm{Temp}$, is:

$$\begin{array}{c}\phi ={{Q}_{10}}^{\frac{\mathrm{Temp}-{\mathrm{T}}_{\mathrm{ref}}}{10}}\end{array}$$

for the squid giant axon, the reference temperature ${\mathrm{T}}_{\mathrm{ref}}$ is $6.3^\circ{\rm C}$, and ${Q}_{10}=3$. In this paper, unless otherwise specified, we set Temp = 18.5 $^\circ{\rm C}$, then $\phi =3.82$. Figure 

Fig. 8

Gating functions of $m$, $h$ and $n$ for the HH model at reference temperature, 6.3 °C

8 shows the steady-state functions for the gating variables, ${m}_{\infty }\left(V\right)$, ${n}_{\infty }\left(V\right)$, ${h}_{\infty }\left(V\right)$ (left), along with time constant functions ${\tau }_{m}\left(V\right)$, ${\tau }_{n}\left(V\right)$, ${\tau }_{h}\left(V\right)$ (right). Notice that the time constant ${\tau }_{m}\left(V\right)$ is much smaller than the other two; $m$ is a fast process.

For $X=m$, $h$ or $n$,

$$\left\{\begin{array}{c}{X}_{\infty }\left(V\right)=\frac{{\alpha }_{X}\left(V\right)}{{\alpha }_{X}\left(V\right)+{\beta }_{X}\left(V\right)}\\ {\tau }_{\infty }\left(V\right)=\frac{1}{{\alpha }_{X}\left(V\right)+{\beta }_{X}\left(V\right)}\end{array}\right.$$

for each specific $X$,

$$\begin{array}{c}{\alpha }_{m}\left(V\right)=0.1\left(V+40\right)/\left(1-\mathrm{exp}\left(-\left(V+40\right)/10\right)\right)\end{array}$$

$$\begin{array}{c}{\beta }_{m}\left(V\right)=4\mathrm{exp}\left(-\left(V+65\right)/18\right)\end{array}$$

$$\begin{array}{c}{\alpha }_{h}\left(V\right)=0.07\mathrm{exp}\left(-\left(V+65\right)/20\right)\end{array}$$

$$\begin{array}{c}{\beta }_{h}\left(V\right)=1/\left(1+\mathrm{exp}\left(-\left(V+35\right)/10\right)\right)\end{array}$$

$$\begin{array}{c}{\alpha }_{n}\left(V\right)=0.01\left(V+55\right)/\left(1-\mathrm{exp}\left(-\left(V+55\right)/10\right)\right)\end{array}$$

$$\begin{array}{c}{\beta }_{n}\left(V\right)=0.125\mathrm{exp}\left(-\left(\mathrm{V}+65\right)/80\right)\end{array}$$

3.2 n–h plane projection

In reducing the HH model to the V-w model, we assumed $h=1-sn$. The validity of this linear relationship is further tested numerically for different steady driving currents (Fig. 9). The trajectories of the reduced model (the V-w model), lie on the line $1-sw,$ and are close to the projections of trajectories of the HH model in the $n-h$ plane, especially when the injected current is small. Meanwhile, the oscillation amplitudes of the HH model are well-predicted by the reduced model. Although the approximation accuracy reduces as the current becomes strong, the impact on our study is slight, since the bistability of interest happens near the onset of oscillation where ${I}_{\rm app}$ is small. Another important condition for the approximation to be valid is that the time constants of the gating variables $h$ and $n$ are comparable.

Fig. 9

HH model trajectories (solid curves with arrows) projected onto the $n - h$ plane for ${I}_{\rm app}$=15 (black), 30 (red), 70 (blue) ${\mu}{A}/{\mathrm{cm}}^{2}$ lie close to the V-w model trajectories that, of course, lie on the anti-diagonal line $\left(1-sw\right) \, \mathrm{vs}\, w$ (gray). The peak and trough of the V-w model’s oscillations (open squares and asterisks, respectively) compare well with those of the HH model

Analytically, the approximation argument can be strengthened under the following assumptions: (1) ${h}_{\infty }\left(V\right)+s{n}_{\infty }\left(V\right)=1$ and (2) ${\tau }_{h}\left(V\right)={\tau }_{n}\left(V\right)$. Then in the HH model, by summing up the $h$-equation and the $n$-equation (multiplied by $s$), we have

$$\frac{dh}{dt}+s\cdot \frac{dn}{dt}=\phi \cdot \left[\frac{{h}_{\infty }\left(V\right)-h}{{\tau }_{h}\left(V\right)}+s\cdot \frac{{n}_{\infty }\left(V\right)-n}{{\tau }_{n}\left(V\right)}\right]$$

$$\Rightarrow \frac{d(h+sn)}{dt}=\phi \cdot \frac{({h}_{\infty }(V)+s{n}_{\infty }(V))-(h+sn)}{{\tau }_{h}(V)}$$

Denote $x=h+sn$. It follows that

$$\begin{aligned}\frac{dx}{dt}=&\phi \cdot \frac{{(h}_{\infty }(V)+s{n}_{\infty }(V))-x}{{\tau }_{h}(V)}\\=&\phi \cdot \frac{1-x}{{\tau }_{h}(V)}=\frac{1-x}{{\tau }_{h}(V)/\phi }\end{aligned}$$

Suppose the initial condition is set such that ${x}_{0}={h}_{0}+s{n}_{0}=1$, then the solution will always lie on the straight line. Otherwise, $h+sn$ converges to 1 exponentially with the effective time constant ${\tau }_{h}(V)/\phi$. That is, the trajectory still converges to the line. It can be further concluded from the analysis that making the slope drive-dependent can yield more precise approximation. Since the accuracy is not crucial for our discussion on bistability, the simple case where $s$ is a constant for all input conditions is used henceforth. Similar analysis was done in Krinskiĭ and Kokoz (1973), but with a different approximation of the linear relationship between $h$ and $n$.

3.3 Comparing critical I_app values for HH and V-w models

Comparing with the HH model, larger oscillation amplitude and higher frequency are seen in the V-w model over most of the ${I}_{\rm app}$ range (Fig. 2E). Moreover, critical values such as the onset of bistable oscillation (${I}_{LP}$), the emergence of an oscillation (${I}_{HB1}$) and the loss of an oscillation (${I}_{HB2}$) are shifted upward in the V-w model. Values of ${I}_{\rm app}$ and features at the Hopf bifurcation for the onset of oscillations of the HH and V-w models are compared in Table.

Table 1 Comparison of critical values of ${I}_{\rm app}$ and the membrane potentials and gating variables of the full HH model and reduced V-w model at the first Hopf bifurcation (${{I}_{\rm app}=I}_{HB1}$). Also listed are values of ${I}_{\rm app}$ at the second Hopf bifurcation and at the limit point of oscillatory solutions (see Results). Temperature: 18.5° C

1.

Notice that the steady-state membrane potentials of the two models are almost the same and thus the values of ${m}_{HB1}$ are close. The differences of the activation of ${I}_{\rm K}$, ${n}_{HB1}$ and ${w}_{HB1}$, can account for the higher ${I}_{\rm app}$ value to achieve same ${V}_{HB1}$ level.

Appendix D: Timescale separation during depolarized and recovery phase in the V-w model

During a spike’s depolarized phase, the trajectory is near to the V-nullcline’s right branch where the outward ${I}_{\rm K}$ and inward ${I}_{\rm Na}$ are both large, on the order of 500 $\mu A/{\rm cm}^{2}$ (Fig. 10 right middle). Somewhat away from the nullcline these opposing currents create a strong leftward or rightward push. Compared with the ionic currents, ${I}_{\rm app}$ is small and can be neglected. The opposing currents are nearly in balance, while their sum is outward ~ 100 $\mu A/{\rm cm}^{2}$(Fig. 10 right bottom), accounting for the subsiding depolarization as the trajectory approaches the upper knee of the right branch. During the spike downstroke, the conductance $\overline{{g}_{\rm K}}{w}^{4}$ for ${I}_{\rm K}$ is maximal and inactivation of ${I}_{\rm Na}$ is most severe. As the trajectory crosses the w-nullcline, the current is strongly outward; ${I}_{\rm Na}$ is almost totally deactivated for $V<-50$ mV and the flow is strongly leftward.

Fig. 10

Ionic current and total current during the recovery and depolarized phase. A. On the left red segment (Fig. 3A), $dV/dt$ is small; slow flow during recovery. Top panel: $V$ vs $t$ (blue) shows slow increase and $w$ vs $t$ (orange), shows w undershooting its steady-state level. Middle panel: the V-gated ionic currents ${I}_{\rm Na}$ (blue), ${I}_{\rm K}$ (orange) are on the order of 50 $\mu A/{cm}^{2};$ ${I}_{L}$ (green) is small and nearly constant; the steady injected current ${I}_{\rm app}$ is purple. Bottom panel: the total current (black), switches from outward to inward and it is small ( $\sim 10 \mu A/{\rm cm}^{2}$) corresponding to small $dV/dt$ and slowly rising $V(t)$. B. $dV/dt$ is of large amplitude on the trajectory’s right red segment (Fig. 3A). In the top panel $V(t)$ decays from its peak value, while ${I}_{\rm K}$ activates and ${I}_{\rm Na}$ inactivates during the spike’s broad depolarized phase. ${I}_{\rm Na}$ and ${I}_{\rm K}$ are large, on the order of 500 $\mu A/{\rm cm}^{2}$, and the total current (lower panel) is $\sim$ 100 $\mu A/{\rm cm}^{2}$, corresponding to $dV/dt$ very negative, responsible for the drop in V of order 100 mV in 1 ms over the red segment. Note that the vertical and horizontal axes are not uniform in panels A and B. The right segment is faster than the left segment due to stronger total current (and consequently larger $dV/dt$)

During a spike’s recovery phase, the trajectory is near the left branch of the V-nullcline, but not hugging the nullcline as during the depolarized phase. With ${I}_{\rm Na}$ negligible, there is a near balance among the small outward current ${I}_{\rm K}$ (small driving force, $V-{E}_{\rm K}$) and small inward current, ${I}_{L}+{I}_{\rm app}$, each of order 10 $\mu A/{\rm cm}^{2}$ (Fig. 10 left middle). Their small sum leads to a small denominator, strong downward flow direction but slow drift of increasing $V$ (Fig. 10 left bottom). The dramatic reduction in $dV/dt$ overwhelms $dw/dt$, creating a near-vertical flow.

Appendix E: Morris–Lecar model

Parameters, critical values and DoB for original Type II ML (Rinzel and Ermentrout 1998) and adjusted ML model are listed in Table.

Table 2 Parameters and critical values for Morris–Lecar model

2.

The model equations for ML model (Morris and Lecar 1981) are given by

$$\begin{aligned}{C}_{m}\frac{dV}{dt}=&{I}_{\rm app}-{g}_{L}\left(V-{E}_{L}\right)-{\overline{g} }_{\rm K}\cdot w\left(V-{E}_{\rm K}\right)\\&-{\overline{g} }_{Ca}\cdot {m}_{\infty }(V)(V-{E}_{Ca})\end{aligned}$$

$$\frac{dw}{dt}=\phi \cdot \frac{{w}_{\infty }\left(V\right)-w}{{\tau }_{w}\left(V\right)}$$

where the gating functions and time constant for K⁺ conductance are

$${m}_{\infty }\left(V\right)=\frac{1}{2}\left[1+\mathrm{tanh}\left(\frac{V-{V}_{1}}{{V}_{2}}\right)\right]$$

$${w}_{\infty }\left(V\right)=\frac{1}{2}\left[1+\mathrm{tanh}\left(\frac{V-{V}_{3}}{{V}_{4}}\right)\right]$$

$${\tau }_{w}\left(V\right)=\frac{1}{\mathrm{cosh}\left(\frac{V-{V}_{3}}{2{V}_{4}}\right)}$$

We use the usual notation $w$ (Rinzel and Ermentrout 1998), for activation of the potassium conductance, but notice that this is not identical with the recovery variable in our reduced HH model.