1628-Problem 1.sagews

Computational Modeling of Neuronal Systems

Homework 1

 

Kenneth L. Ho

[email protected]

version()

'Sage Version 4.3.2, Release Date: 2010-02-06'

import pylab
from scipy import fftpack
from scipy import integrate
from scipy import optimize

Choose two out of the following three exercises on the Hodgkin-Huxley (HH) model. You may use the XPP or Matlab codes on our course website. Choose a or for numerical integration so that tightening the criterion gives comparable results. State the error criteria and integration method that you use.

 

Problem 1

Do some simulations with the standard HH model; set the temperature to $18.5$ $^{\circ}$C.

1. Use an $I_{\mathrm{app}}$ amplitude for a $1$-ms duration that produces a spike, and a weaker stimulus that produces a just-subthreshold response. Show time course plots of $V$ and of $m$, $h$, and $n$ (these three on one set of axes), and show phase plane projections: $V$ vs. $n$ and $h$ vs. $n$. Describe the features of the response as viewed in the phase plane.
2. Use a long duration $I_{\mathrm{app}}$, say, $100$ ms, to see repetitive firing for $I_{\mathrm{app}} = 20$ $\mu$A/cm$^{2}$. Produce corresponding plots as for (a) and describe the features.
3. Develop a three-variable HH model by setting $m = m_{\infty} (V)$. Redo (a) and (b), and compare your results to the full model. Describe any differences.
4. Make a two-variable HH model by setting $m = m_{\infty} (V)$ and $h = 1 - n$. Redo (a) and (b), and compare your results to the full model. Describe any differences.
5. If you use XPP, use the command to compute the steady state and stability properties (including the eigenvalues) of the three versions of the model for $I_{\mathrm{app}} = 0$ and $20$ $\mu$A/cm$^{2}$. Compare and discuss the results.

 

Solution

We implement the core functions for the required HH models below in Python, making use of Cython for speed, following the Matlab codes provided.

%cython

import numpy
cimport numpy
from numpy.linalg import norm

cdouble = numpy.double
ctypedef numpy.double_t cdouble_t

cpdef tuple m_dynamics(double V):
    cdef double alpha = ((V + 45)/10) / (1 - numpy.exp(-(V + 45)/10))
    cdef double beta = 4*numpy.exp(-(V + 70)/18)
    cdef double m_inf = alpha / (alpha + beta)
    cdef double tau_m = 1 / (alpha + beta)
    return m_inf, tau_m

cpdef tuple h_dynamics(double V):
    cdef double alpha = 0.07*numpy.exp(-(V + 70)/20)
    cdef double beta = 1 / (1 + numpy.exp(-(V + 40)/10))
    cdef double h_inf = alpha / (alpha + beta)
    cdef double tau_h = 1 / (alpha + beta)
    return h_inf, tau_h

cpdef tuple n_dynamics(double V):
    cdef double alpha = 0.1*((V + 60)/10) / (1 - numpy.exp(-(V + 60)/10))
    cdef double beta = 0.125*numpy.exp(-(V + 70)/80)
    cdef double n_inf = alpha / (alpha + beta)
    cdef double tau_n = 1 / (alpha + beta)
    return n_inf, tau_n

cpdef double I_inst(double V, double m, double h, double n, \
                    double g_Na, double g_K, double g_L,    \
                    double E_Na, double E_K, double E_L):
    return g_Na*m**3*h*(V - E_Na) + g_K*n**4*(V - E_K) + g_L*(V - E_L)

cpdef numpy.ndarray hh_ode(double t, numpy.ndarray y,                     \
                           double C, double g_Na, double g_K, double g_L, \
                           double E_Na, double E_K, double E_L, I):
    cdef double V = y[0]
    cdef double m = y[1]
    cdef double h = y[2]
    cdef double n = y[3]
    cdef double m_inf, tau_m, h_inf, tau_h, n_inf, tau_n
    m_inf, tau_m = m_dynamics(V)
    h_inf, tau_h = h_dynamics(V)
    n_inf, tau_n = n_dynamics(V)
    return numpy.array([(I(t) - I_inst(V, m, h, n, g_Na, g_K, g_L, E_Na, E_K, E_L))/C, \
                        (m_inf - m)/tau_m,                                             \
                        (h_inf - h)/tau_h,                                             \
                        (n_inf - n)/tau_n])

cpdef numpy.ndarray hh2_ode(double t, numpy.ndarray y,                     \
                            double C, double g_Na, double g_K, double g_L, \
                            double E_Na, double E_K, double E_L, I):
    cdef double V = y[0]
    cdef double n = y[1]
    cdef double m = m_dynamics(V)[0]
    cdef double h = 1 - n
    cdef double n_inf, tau_n
    n_inf, tau_n = n_dynamics(V)
    return numpy.array([(I(t) - I_inst(V, m, h, n, g_Na, g_K, g_L, E_Na, E_K, E_L))/C, \
                        (n_inf - n)/tau_n])

cpdef numpy.ndarray hh3_ode(double t, numpy.ndarray y,                     \
                            double C, double g_Na, double g_K, double g_L, \
                            double E_Na, double E_K, double E_L, I):
    cdef double V = y[0]
    cdef double h = y[1]
    cdef double n = y[2]
    cdef double m = m_dynamics(V)[0]
    cdef double h_inf, tau_h, n_inf, tau_n
    h_inf, tau_h = h_dynamics(V)
    n_inf, tau_n = n_dynamics(V)
    return numpy.array([(I(t) - I_inst(V, m, h, n, g_Na, g_K, g_L, E_Na, E_K, E_L))/C, \
                        (h_inf - h)/tau_h,                                             \
                        (n_inf - n)/tau_n])

cpdef dt_rel(deriv, double t, numpy.ndarray y, double dt_min, double dt_max, \
             double C, double g_Na, double g_K, double g_L,                  \
             double E_Na, double E_K, double E_L, I):
    cdef numpy.ndarray[cdouble_t] dy = deriv(t, y, C, g_Na, g_K, g_L, E_Na, E_K, E_L, I)
    cdef double a = norm(dy / y)
    return dt_min*(a/(1 + a)) + dt_max*(1/(1 + a))

This is followed by a general driver, which can be used to run simulations of the full, three-, and two-variable models, as desired.

def hh(deriv, y0, tspan, dt_min, dt_max, *params):
    y = []
    t = []
    y.append(y0)
    t.append(tspan[0])
    r = integrate.ode(deriv)
    r.set_initial_value(y[0], t[0])
    r.set_f_params(*params)
    while r.successful() and r.t < tspan[1]:
        dt = dt_rel(deriv, t[-1], pylab.array(y[-1]), dt_min, dt_max, *params)
        r.integrate(r.t + dt)
        t.append(r.t)
        y.append(r.y)
    return pylab.array(t), pylab.array(y)

def hh_plot(t, y):
    pylab.rc('figure', figsize=(6,4))
    pylab.rc('legend', fontsize='x-small')
    pylab.rc('xtick', labelsize='xx-small')
    pylab.rc('ytick', labelsize='xx-small')
    pylab.figure()
    pylab.subplot(221)
    pylab.plot(t, y[:,0])
    pylab.xlim(tspan)
    pylab.xlabel(r'Time ($t$) [ms]', size='small')
    pylab.ylabel(r'Voltage ($V$) [mV]', size='small')
    pylab.subplot(222)
    pylab.plot(t, y[:,1], 'b', label=r'Na + ($m$)')
    pylab.plot(t, y[:,2], 'g', label=r'Na - ($h$)')
    pylab.plot(t, y[:,3], 'r', label=r'K + ($n$)')
    pylab.xlim(tspan)
    pylab.legend(loc='best')
    pylab.xlabel(r'Time ($t$) [ms]', size='small')
    pylab.ylabel(r'Gating variables', size='small')
    pylab.subplot(223)
    pylab.plot(y[:,0], y[:,3])
    pylab.xlabel(r'Voltage ($V$) [mV]', size='small')
    pylab.ylabel(r'K activation ($n$)', size='small')
    pylab.subplot(224)
    pylab.plot(y[:,2], y[:,3])
    pylab.xlabel(r'Na inhibition ($h$)', size='small')
    pylab.ylabel(r'K activation ($n$)', size='small')
    pylab.subplots_adjust(wspace=0.3, hspace=0.3)
    pylab.savefig('sage.png')
    pylab.close()
    pylab.rcdefaults()

def proc_hh2(y):
    n = len(y)
    m = pylab.empty(n)
    for i in (0..n-1): m[i] = m_dynamics(y[i,0])[0]
    y = pylab.hstack((y, pylab.transpose([m]), pylab.transpose([1-y[:,1]])))
    y[:,1:4] = y[:,[2,3,1]]
    return y

def proc_hh3(y):
    n = len(y)
    m = pylab.empty(n)
    for i in (0..n-1): m[i] = m_dynamics(y[i,0])[0]
    y = pylab.hstack((y, pylab.transpose([m])))
    y[:,1:4] = y[:,[3,1,2]]
    return y

A note on the numerics: we use the solver , as wrapped in SciPy, for numerical integration. As performs automatic timestep adjustment with regard to numerical stability, we need only specify the timestep with regard to time resolution. To do this, we use the criterion $$\Delta t = \left( \frac{\left\| \Delta \right\|}{1 + \left\| \Delta \right\|} \right) \Delta t_{\min} + \left( \frac{1}{1 + \left\| \Delta \right\|} \right) \Delta t_{\max},$$ where $\Delta_{i} = \dot{x}_{i} / x_{i}$, i.e., $\Delta$ is the vector of relative derivatives; and where $\Delta t_{\min}$ and $\Delta t_{\max}$ are parameters.

To initialize the models, we use the parameters from the Matlab codes; the initial condition is set by assuming a resting potential of $V_{\mathrm{rest}} = -70$ mV and setting $m$, $h$, and $n$ to their corresponding equilibrium values at $V_{\mathrm{rest}}$. For the applied current, we use the wave form $$I_{\mathrm{app}} \left( t; A, \tau \right) = \begin{cases} A, & t \leq \tau,\\ 0, & t > \tau. \end{cases}$$

C = 1

g_Na = 120
g_K  =  36
g_L  =   0.3

E_Na =  45
E_K  = -82
E_L  = -59

V_rest  = -70
m_rest = m_dynamics(V_rest)[0]
h_rest = h_dynamics(V_rest)[0]
n_rest = n_dynamics(V_rest)[0]
y0 = [V_rest, m_rest, h_rest, n_rest]

def I(t, A, tau):
    if t <= tau: return A
    else:        return 0

The first simulation (a) is of the full model for a spike in response to a short stimulus. To do this, we use $A = 20$ $\mu$A/cm$^{2}$ and $\tau = 1$ ms. For resolution, we set $\Delta t_{\min} = 0.01$ and $\Delta t_{\max} = 0.1$.

tspan = [0, 30]
dt_min = 1e-2
dt_max = 1e-1

A = 20
tau = 1

t, y = hh(hh_ode, y0, tspan, dt_min, dt_max, \
          C, g_Na, g_K, g_L, E_Na, E_K, E_L, \
          lambda t: I(t, A, tau))
hh_plot(t, y)

As desired, we obtain a spike, with characteristic time courses for both the membrane potential $V$ and the gating variables $m$, $h$, and $n$. The spiking behavior is especially evident in the $V$-$n$ phase plane, which produces a full cycle of the relaxation dynamics. Futhermore, interestingly, we observe a strong linear dependence between the two variables in the $n$-$h$ phase plane, suggesting a possible model reduction as we will see later.

The next simulation is of the full model for a short subthreshold stimulation ($A = 4.5$ $\mu$A/cm$^{2}$, $\tau = 1$ ms).

tspan = [0, 30]
dt_min = 1e-2
dt_max = 1e-1

A = 4.5
tau = 1

t, y = hh(hh_ode, y0, tspan, dt_min, dt_max, \
          C, g_Na, g_K, g_L, E_Na, E_K, E_L, \
          lambda t: I(t, A, tau))
hh_plot(t, y)

Indeed, we do not produce a spike, with the membrane voltage remaining around $V_{\mathrm{rest}} \approx -70$ mV, and the gating variables more or less constant throughout. In the $n$-$V$ phase plane, we see an attempt by the system to leave the rest state, as in (a), but it fails to cross the threshold and quickly restabilizes. The $n$-$h$ phase plane is similar to that previously obtained.

We now run a simulation (b) of the full model in response to a long stimulus ($A = 20$ $\mu$A/cm$^{2}$, $\tau = 100$ ms).

tspan = [0, 100]
dt_min = 1e-2
dt_max = 1e-1

A = 20
tau = 100

t, y = hh(hh_ode, y0, tspan, dt_min, dt_max, \
          C, g_Na, g_K, g_L, E_Na, E_K, E_L, \
          lambda t: I(t, A, tau))
hh_plot(t, y)

Repetitive firing is observed, though, interestingly, the first spike appears larger than those ensuing. Biochemically, this could be due to, e.g., synaptic depression, though this, of course, does not appear explicitly in the model. The phase planes show, generically, multiple passes over the trajectories in (a), with the system locking onto a limit cycle sustained by the elevated $I_{\mathrm{app}}$ as evident in the $n$-$V$ phase plane.

The fourth simulation (c) is of the three-variable model, obtained by setting $m = m_{\infty} (V)$. For the remainder of this problem, we consider only the short and long spike-generating stimuli ($A = 20$ $\mu$A/cm$^{2}$ with $\tau = 1$ ms and $100$ ms, respectively). For the short stimulus, we obtain the following:

tspan = [0, 30]
dt_min = 1e-2
dt_max = 1e-1

A = 20
tau = 1

t, y = hh(hh3_ode, [y0[0], y0[2], y0[3]], tspan, dt_min, dt_max, \
          C, g_Na, g_K, g_L, E_Na, E_K, E_L,                     \
          lambda t: I(t, A, tau))
y = proc_hh3(y)
hh_plot(t, y)

The plots are generally similar to those obtained in (a), suggesting a justified model reduction, though $m$ opens and hence allows the cell to fire a bit sooner. The time courses of $h$ and $n$ are, in response, a bit sharper to account for the instantaneous relaxation introduced in $m$. Furthermore, the $n$-$V$ phase plane appears rotated relative to that in (a), no doubt a result of the instantaneous dynamics.

Likewise, the long stimulus gives:

tspan = [0, 100]
dt_min = 1e-2
dt_max = 1e-1

A = 20
tau = 100

t, y = hh(hh3_ode, [y0[0], y0[2], y0[3]], tspan, dt_min, dt_max, \
          C, g_Na, g_K, g_L, E_Na, E_K, E_L,                     \
          lambda t: I(t, A, tau))
y = proc_hh3(y)
hh_plot(t, y)

which, again, is quite similar to that obtained in (b). Similar comments apply as in the short stimulus case. Moreover, effective synaptic depression appears to be weaker.

Finally, we run simulations (d) of the two-variable model, obtained by further setting $h = 1 - n$, as suggested by essentially all of the above $n$-$h$ phase planes.

tspan = [0, 30]
dt_min = 1e-2
dt_max = 1e-1

A = 20
tau = 1

t, y = hh(hh2_ode, [y0[0], y0[3]], tspan, dt_min, dt_max, \
          C, g_Na, g_K, g_L, E_Na, E_K, E_L,              \
          lambda t: I(t, A, tau))
y = proc_hh2(y)
hh_plot(t, y)

For the short stimulus, we obtain fairly good agreement with previous results, the primary point of contention (apart from the obvious collapse of the $n$-$h$ phase plane) being that the $V$ time course appears segmented on the downstroke. Note that the $n$-$h$ phase plane produced here is an artifact of the matplotlib graphics system; inspection of the actual data reveals a sustained line much like that to appear below.

tspan = [0, 100]
dt_min = 1e-2
dt_max = 1e-1

A = 20
tau = 100

t, y = hh(hh2_ode, [y0[0], y0[3]], tspan, dt_min, dt_max, \
          C, g_Na, g_K, g_L, E_Na, E_K, E_L,              \
          lambda t: I(t, A, tau))
y = proc_hh2(y)
hh_plot(t, y)

The long stimulus produces results very similar to those in (c). Most striking, though, is the loss of variability in the limit cycle in the $n$-$V$ phase plane.

We study the steady states (e) of the three models above by numerically solving $I_{\mathrm{inst}} \left( V, m, h, n \right) = I_{\mathrm{app}}$ for $V$, where $$I_{\mathrm{inst}} \left( V, m, h, n \right) = g_{\mathrm{Na}} m^{3} h \left( V - E_{\mathrm{Na}} \right) + g_{\mathrm{K}} n^{4} \left( V - E_{\mathrm{K}} \right) + g_{\mathrm{L}} \left( V - E_{\mathrm{L}} \right)$$ is the instantaneous current through the membrane, on setting $m = m_{\infty} (V)$ and $n = n_{\infty} (V)$, and $h = h_{\infty} (V)$ for the full and three-variable models, and $h = 1 - n_{\infty} (V)$ for the two-variable model. The stabilities of the steady states are evaluated by symbolically constructing the Jacobian $J$ of the system $$\begin{align*} \dot{V} &= \frac{I_{\mathrm{app}} - I_{\mathrm{inst}} \left( V, m, h, n \right)}{C},\\ \dot{m} &= \alpha_{m} \left( V \right) \left( 1 - m \right) - \beta_{m} \left( V \right) m,\\ \dot{h} &= \alpha_{h} \left( V \right) \left( 1 - h \right) - \beta_{h} \left( V \right) h,\\ \dot{n} &= \alpha_{n} \left( V \right) \left( 1 - n \right) - \beta_{n} \left( V \right) n, \end{align*}$$ and computing its eigenvalues at the steady states.

The steady-state potential $V_{\infty}$ satisfying $I_{\mathrm{inst}} = I_{\mathrm{app}}$ is obtained using MINPACK, and the eigenvalues of $J$ computed using LAPACK, both as wrapped by NumPy/SciPy.

def V_ss_eqn(nvar, g_Na, g_K, g_L, E_Na, E_K, E_L, I):
    def V_ss(V):
        m = m_dynamics(V)[0]
        n = n_dynamics(V)[0]
        if nvar == 3 or nvar == 4: h = h_dynamics(V)[0]
        elif nvar == 2:            h = 1 - n
        return I_inst(V, m, h, n, g_Na, g_K, g_L, E_Na, E_K, E_L) - I
    return V_ss

def hh_ss(nvar, x0, g_Na, g_K, g_L, E_Na, E_K, E_L, I):
    V = optimize.fsolve(V_ss_eqn(nvar, g_Na, g_K, g_L, E_Na, E_K, E_L, I), x0)
    m = m_dynamics(V)[0]
    n = n_dynamics(V)[0]
    if nvar == 3 or nvar == 4: h = h_dynamics(V)[0]
    elif nvar == 2:            h = 1 - n
    return V, m, h, n

V, m, h, n = var('V m h n')
gNa, gK, gL = var('gNa gK gL')
ENa, EK, EL = var('ENa EK EL')
I_sym = gNa*m^3*h*(V - ENa) + gK*n^4*(V - EK) + gL*(V - EL)

def J_ss(nvar, V_ss, C, g_Na, g_K, g_L, E_Na, E_K, E_L):
    alpha_m = ((V + 45)/10) / (1 - exp(-(V + 45)/10))
    beta_m = 4*exp(-(V + 70)/18)
    m_inf = alpha_m / (alpha_m + beta_m)
    dm = alpha_m*(1 - m) - beta_m*m
    alpha_n = 0.1*((V + 60)/10) / (1 - exp(-(V + 60)/10))
    beta_n = 0.125*exp(-(V + 70)/80)
    n_inf = alpha_n / (alpha_n + beta_n)
    dn = alpha_n*(1 - n) - beta_n*n
    if nvar == 3 or nvar == 4:
        alpha_h = 0.07*exp(-(V + 70)/20)
        beta_h = 1 / (1 + exp(-(V + 40)/10))
        h_inf = alpha_h / (alpha_h + beta_h)
        dh = alpha_h*(1 - h) - beta_h*h
    elif nvar == 2: h_inf = 1 - n
    ss = {V  : V_ss,               \
          m  : m_inf.subs(V=V_ss), \
          h  : h_inf.subs(V=V_ss), \
          n  : n_inf.subs(V=V_ss), \
          gNa: g_Na,               \
          gK : g_K,                \
          gL : g_L,                \
          ENa: E_Na,               \
          EK : E_K,                \
          EL : E_L}
    if nvar == 4:
        dV = -I_sym / C
        J = [[f.diff(x).subs(ss) for x in [V, m, h, n]] for f in [dV, dm, dh, dn]]
    elif nvar == 3:
        dV = -I_sym.subs(m=m_inf) / C
        J = [[f.diff(x).subs(ss) for x in [V, h, n]] for f in [dV, dh, dn]]
    elif nvar == 2:
        dV = -I_sym.subs({m: m_inf, h: 1-n}) / C
        J = [[f.diff(x).subs(ss) for x in [V, n]] for f in [dV, dn]]
    return pylab.array(J, dtype=float)

We compute the steady states and their associated eigenvalues for $I_{\mathrm{app}} = 0$ and $20$ $\mu$A/cm$^{2}$. For the full model, these are:

html('full model:')
print
for I in [0, 20]:
    V_ss, m_ss, h_ss, n_ss = hh_ss(4, 0, g_Na, g_K, g_L, E_Na, E_K, E_L, I)
    html(r'$(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty})$ at $I_{\mathrm{app}} = %i$:' % I)
    pretty_print((V_ss, m_ss, h_ss, n_ss))
    print
    html('eigenvalues of $J$ at $I_{\mathrm{app}} = %i$:' % I)
    pretty_print(pylab.eigvals(J_ss(4, V_ss, C, g_Na, g_K, g_L, E_Na, E_K, E_L)))
    print
    print

full model:

(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty}) at I_{\mathrm{app}} = 0:

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\hbox{-69.8976728964}, 0.0535746092323, 0.592537659007, 0.319246167222\right)

eigenvalues of J at I_{\mathrm{app}} = 0:

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{[-4.66655603+0.j -0.20044489+0.3876745j -0.20044489-0.3876745j -0.12082016+0.j ]}

(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty}) at I_{\mathrm{app}} = 20:

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\hbox{-61.5643150194}, 0.13497383385, 0.306805365841, 0.450925338254\right)

eigenvalues of J at I_{\mathrm{app}} = 20:

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{[-5.28245933+0.j 0.15640342+0.64195044j 0.15640342-0.64195044j -0.15780123+0.j ]}

i.e., the steady state at $I_{\mathrm{app}} = 0$ $\mu$A/cm$^{2}$ is stable, while that at $I_{\mathrm{app}} = 20$ $\mu$A/cm$^{2}$ is unstable. There is evidently only one steady state as a later plot will show. The results for the three-variable model are similar:

html('three-variable model:')
print
for I in [0, 20]:
    V_ss, m_ss, h_ss, n_ss = hh_ss(3, 0, g_Na, g_K, g_L, E_Na, E_K, E_L, I)
    html(r'$(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty})$ at $I_{\mathrm{app}} = %i$:' % I)
    pretty_print((V_ss, m_dynamics(V_ss)[0], h_ss, n_ss))
    print
    html('eigenvalues of $J$ at $I_{\mathrm{app}} = %i$:' % I)
    pretty_print(pylab.eigvals(J_ss(3, V_ss, C, g_Na, g_K, g_L, E_Na, E_K, E_L)))
    print
    print

three-variable model:

(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty}) at I_{\mathrm{app}} = 0:

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\hbox{-69.8976728964}, 0.0535746092323, 0.592537659007, 0.319246167222\right)

eigenvalues of J at I_{\mathrm{app}} = 0:

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{[-0.21057224+0.40888856j -0.21057224-0.40888856j -0.12080474+0.j ]}

(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty}) at I_{\mathrm{app}} = 20:

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\hbox{-61.5643150194}, 0.13497383385, 0.306805365841, 0.450925338254\right)

eigenvalues of J at I_{\mathrm{app}} = 20:

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{[ 0.44882718+0.77185623j 0.44882718-0.77185623j -0.15773290+0.j ]}

For both of these models, therefore, the steady applied current induces a bifurcation from a steady rest state to an oscillating limit cycle, producing repetitive firing as has already been seen earlier.

For the two-variable model, however, we get something a little different. In particular, there are now three steady states: first,

html(r'two-variable model near $V_{\infty} \approx -70$ mV:')
print
for I in [0, 20]:
    V_ss, m_ss, h_ss, n_ss = hh_ss(2, -70, g_Na, g_K, g_L, E_Na, E_K, E_L, I)
    html(r'$(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty})$ at $I_{\mathrm{app}} = %i$:' % I)
    pretty_print((V_ss, m_dynamics(V_ss)[0], 1 - n_ss, n_ss))
    print
    html('eigenvalues of $J$ at $I_{\mathrm{app}} = %i$:' % I)
    pretty_print(pylab.eigvals(J_ss(2, V_ss, C, g_Na, g_K, g_L, E_Na, E_K, E_L)))
    print
    print

two-variable model near V_{\infty} \approx -70 mV:

(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty}) at I_{\mathrm{app}} = 0:

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\hbox{-69.7262906627}, 0.0546657944747, 0.678120453084, 0.321879546916\right)

eigenvalues of J at I_{\mathrm{app}} = 0:

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{[-0.17284854+0.41570498j -0.17284854-0.41570498j]}

Warning: The iteration is not making good progress, as measured by the improvement from the last ten iterations.

(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty}) at I_{\mathrm{app}} = 20:

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\hbox{-58.9351936326}, 0.175377559045, 0.508015695861, 0.491984304139\right)

eigenvalues of J at I_{\mathrm{app}} = 20:

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{[ 6.92443132e+00 -3.60566769e-04]}

second:

html(r'two-variable model near $V_{\infty} \approx -50$ mV:')
print
for I in [0, 20]:
    V_ss, m_ss, h_ss, n_ss = hh_ss(2, -50, g_Na, g_K, g_L, E_Na, E_K, E_L, I)
    html(r'$(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty})$ at $I_{\mathrm{app}} = %i$:' % I)
    pretty_print((V_ss, m_dynamics(V_ss)[0], 1 - n_ss, n_ss))
    print
    html('eigenvalues of $J$ at $I_{\mathrm{app}} = %i$:' % I)
    pretty_print(pylab.eigvals(J_ss(2, V_ss, C, g_Na, g_K, g_L, E_Na, E_K, E_L)))
    print
    print

two-variable model near V_{\infty} \approx -50 mV:

(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty}) at I_{\mathrm{app}} = 0:

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\hbox{-54.0694295866}, 0.27115411034, 0.435878786546, 0.564121213454\right)

eigenvalues of J at I_{\mathrm{app}} = 0:

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{[ 20.38459968 -0.08022857]}

Warning: The iteration is not making good progress, as measured by the improvement from the last five Jacobian evaluations.

(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty}) at I_{\mathrm{app}} = 20:

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\hbox{-58.9502980712}, 0.175122534275, 0.508248314858, 0.491751685142\right)

eigenvalues of J at I_{\mathrm{app}} = 20:

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{[ 6.89629867e+00 1.62786029e-05]}

and third:

html(r'two-variable model near $V_{\infty} \approx -25$ mV:')
print
for I in [0, 20]:
    V_ss, m_ss, h_ss, n_ss = hh_ss(2, -25, g_Na, g_K, g_L, E_Na, E_K, E_L, I)
    html(r'$(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty})$ at $I_{\mathrm{app}} = %i$:' % I)
    pretty_print((V_ss, m_dynamics(V_ss)[0], 1 - n_ss, n_ss))
    print
    html('eigenvalues of $J$ at $I_{\mathrm{app}} = %i$:' % I)
    pretty_print(pylab.eigvals(J_ss(2, V_ss, C, g_Na, g_K, g_L, E_Na, E_K, E_L)))
    print
    print

two-variable model near V_{\infty} \approx -25 mV:

(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty}) at I_{\mathrm{app}} = 0:

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\hbox{-26.5156451255}, 0.859989720941, 0.17297674804, 0.82702325196\right)

eigenvalues of J at I_{\mathrm{app}} = 0:

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{[ 2.65557338+3.70630697j 2.65557338-3.70630697j]}

(V_{\infty}, m_{\infty}, h_{\infty}, n_{\infty}) at I_{\mathrm{app}} = 20:

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\hbox{-26.1174195632}, 0.864288169081, 0.170787639316, 0.829212360684\right)

eigenvalues of J at I_{\mathrm{app}} = 20:

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{[ 1.88345651+4.23974391j 1.88345651-4.23974391j]}

The first steady state, again, shares the same stable-unstable transition (Hopf bifurcation) as the other two models, while the second and third steady states are inherently unstable, generating natural thresholds for spiking behavior. Note that the first and second activated ($I_{\mathrm{app}} = 20$ $\mu$A/cm$^{2}$) steady states may be switched, judging based on the value of $V_{\infty}$, though the character of the eigenvalues may suggest otherwise.

For justification of the number of steady states that each model admits, we appeal to the following graphical analysis, as obtained by plotting $I_{\mathrm{inst}} (V)$ with $m$, $h$, and $n$ set accordingly at their steady-state values.

V_spn = pylab.linspace(-80, -20, 500)
pylab.rc('figure', figsize=(6.5,3))
pylab.rc('xtick', labelsize='x-small')
pylab.rc('ytick', labelsize='x-small')
pylab.figure()
pylab.subplot(131)
pylab.plot(V_spn, [V_ss_eqn(4, g_Na, g_K, g_L, E_Na, E_K, E_L, 0)(V) for V in V_spn])
pylab.axhline(c='k', ls='--')
pylab.ylabel(r'Current ($I_{\mathrm{inst}}$) [$\mu$A/cm$^{2}$]', size='small')
pylab.title('Full model', size='small')
pylab.subplot(132)
pylab.plot(V_spn, [V_ss_eqn(3, g_Na, g_K, g_L, E_Na, E_K, E_L, 0)(V) for V in V_spn])
pylab.axhline(c='k', ls='--')
pylab.xlabel(r'Voltage ($V$) [mV]', size='small')
pylab.title('Three-var model', size='small')
pylab.subplot(133)
pylab.plot(V_spn, [V_ss_eqn(2, g_Na, g_K, g_L, E_Na, E_K, E_L, 0)(V) for V in V_spn])
pylab.axhline(c='k', ls='--')
pylab.title('Two-var model', size='small')
pylab.subplots_adjust(bottom=0.15, wspace=0.3)
pylab.savefig('sage.png')
pylab.close()
pylab.rcdefaults()

The monotonic nature of the curves for the full and three-variable models gives clearly that only one steady state exists, whereas the more complicated shape of that of the two-variable model suggests the existence of three steady states.

 

Problem 2

Consider HH without $I_{\mathrm{K}}$ (i.e., $g_{\mathrm{K}} = 0$). Show that with adjustments in $g_{\mathrm{Na}}$ (and maybe $g_{\mathrm{L}}$), the HH model is still excitable and generates an action potential [do it with $m = m_{\infty} (V)$]. Study this two-variable ($V$-$h$) model in the phase plane: nullclines, stability of rest state, trajectories, etc. Then consider a range of $I_{\mathrm{app}}$ to see if you get repetitive firing. Compute the frequency vs. $I_{\mathrm{app}}$ relation; study in the phase plane. Do analysis to see that the rest point must be on the middle branch to get a limit cycle.

 

Solution

The model setup is exactly that of the three-variable model but with $g_{\mathrm{K}} = 0$, hence we will make use of the three-variable model in the following. To facilitate analysis, we first implement code for computing nullclines in the $h$-$V$ phase plane.

Iapp = var('Iapp')
h_nc = (Iapp - I_sym).solve(h)[0].rhs()

def hh_noK_pplot(t, y, g_Na, g_L, E_Na, E_L, I):
    def V_nc(V):
        return h_nc.subs({var('V'): V,                \
                          m       : m_dynamics(V)[0], \
                          n       : 0,                \
                          gNa     : g_Na,             \
                          gK      : 0,                \
                          gL      : g_L,              \
                          ENa     : E_Na,             \
                          EK      : 0,                \
                          EL      : E_L,              \
                          Iapp    : I}).n()
    pylab.rc('figure', figsize=(5.5,3))
    pylab.rc('legend', fontsize='x-small')
    pylab.rc('xtick', labelsize='x-small')
    pylab.rc('ytick', labelsize='x-small')
    pylab.figure()
    pylab.subplot(121)
    pylab.plot(t, y[:,0])
    pylab.xlim(tspan)
    pylab.xlabel(r'Time ($t$) [ms]', size='small')
    pylab.ylabel(r'Voltage ($V$) [mV]', size='small')
    pylab.subplot(122)
    pylab.plot(y[:,0], y[:,2])
    V_lim = pylab.xlim()
    ylim = pylab.ylim()
    V_spn = pylab.linspace(*V_lim)
    pylab.plot(V_spn, [h_dynamics(V)[0] for V in V_spn], label=r'$\dot{h} = 0$')
    pylab.plot(V_spn, [V_nc(V) for V in V_spn], label=r'$\dot{V} = 0$')
    pylab.legend(loc='best')
    pylab.ylim(ylim)
    pylab.xlabel(r'Voltage ($V$) [mV]', size='small')
    pylab.ylabel(r'Na inhibition ($h$)', size='small')
    pylab.subplots_adjust(bottom=0.2, wspace=0.3)
    pylab.savefig('sage.png')
    pylab.close()
    pylab.rcdefaults()

We then run the three-variable model with zero input ($A = 0$ $\mu$A/cm$^{2}$), using the previous values of $g_{\mathrm{Na}}$ and $g_{\mathrm{L}}$.

tspan = [0, 100]
dt_min = 1e-2
dt_max = 1e-1

A = 0
tau = 100
    
t, y = hh(hh3_ode, [y0[0], y0[2], y0[3]], tspan, dt_min, dt_max, \
          C, g_Na, 0, g_L, E_Na, E_K, E_L,                       \
          lambda t: I(t, A, tau))
y = proc_hh3(y)
hh_noK_pplot(t, y, g_Na, g_L, E_Na, E_L, A)

As observed, the system spontaneously generates one spike, then stabilizes to a depolarized state. This may clearly be attributed to the lack of potassium outflux to repolarize the membrane. Hence to restore proper spiking behavior, we must intuitively decrease $g_{\mathrm{Na}}$ and increase $g_{\mathrm{L}}$. We do this by setting $\tilde{g}_{\mathrm{Na}} = \varphi_{\mathrm{Na}} g_{\mathrm{Na}}$ and $\tilde{g}_{\mathrm{L}} = \varphi_{\mathrm{L}} g_{\mathrm{L}}$, where, for illustration, we choose $\varphi_{\mathrm{Na}} = 0.5$ and $\varphi_{\mathrm{L}} = 10$, and use $\tilde{g}_{\mathrm{Na}}$ and $\tilde{g}_{\mathrm{L}}$ in place of $g_{\mathrm{Na}}$ and $g_{\mathrm{L}}$. For example, using a stimulus with $A = 10$ $\mu$A/cm$^{2}$ and $\tau = 100$ ms gives repetitive firing:

tspan = [0, 100]
dt_min = 1e-2
dt_max = 1e-1

A = 10
tau = 100

phi_Na = 0.5
phi_L = 10

t, y = hh(hh3_ode, [y0[0], y0[2], y0[3]], tspan, dt_min, dt_max, \
          C, phi_Na*g_Na, 0, phi_L*g_L, E_Na, E_K, E_L,          \
          lambda t: I(t, A, tau))
y = proc_hh3(y)
hh_noK_pplot(t, y, phi_Na*g_Na, phi_L*g_L, E_Na, E_L, A)

The phase plane is particular illuminating, demonstrating the limit cycle about the unstable steady state. We note here that we have plotted only the nullclines as the stability of the rest state is obvious given these and the trajectory.

Furthermore, that repetitive firing requires the rest state to be on the middle downward branch of the $V$-nullcline is also immediate given the diagram. To see this, we adopt a geometric view, and suppose for simplicity that all trajectories start from rest, i.e., the upper left of the above phase plane; moreover, we assume that the dynamics in $V$ are fast, as supported by the plots. Then if the steady state lies on the rightmost branch of the $V$-nullcline, the trajectory cannot cross and so hugs the $V$-nullcline until it collides with the steady state, muting any potential oscillation. Similarly, if the steady state lies on the leftmost branch, the same argument applies after the trajectory crosses over past the bend of the middle branch. Therefore, a limit cycle is possible only if the steady state lies on the middle branch, in which case, we have, generically, oscillations between the left and right branches of the $V$-nullcline, with each escape made possible by the corresponding bends defining the extent of the middle branch.

To compute the frequency $\omega$ for a given repetitive firing trajectory, we make use of the fast Fourier transform (FFT). The following function uniformly resamples the data, then uses the FFT to pick out the dominant $\omega$.

def hh_freq(t, y, m):
    n = len(t)
    t_reg = pylab.linspace(t[0], t[-1], m)
    y_reg = pylab.interp(t_reg, t, y)
    omega = pylab.linspace(0, m-1, m) / (t[-1] - t[0])
    f = abs(fftpack.fft(y_reg))
    omega = omega[pylab.find(f[1:] == max(f[1:]))[0] + 1]
    return omega

Initial numerical explorations indicate that repetitive firing is observed from $A \approx -3$ to $20$ $\mu$A/cm$^{2}$, so we use this as our parameter range.

tspan = [0, 100]
dt_min = 1e-2
dt_max = 1e-1

A = pylab.linspace(-3, 20)
tau = 100

phi_Na = 0.5
phi_L = 10

omega = []
for a in A:
    t, y = hh(hh3_ode, [y0[0], y0[2], y0[3]], tspan, dt_min, dt_max, \
              C, phi_Na*g_Na, 0, phi_L*g_L, E_Na, E_K, E_L,          \
              lambda t: I(t, a, tau))
    omega.append(1e3 * hh_freq(t, y[:,0], 1e4))

The frequency relation is shown below, from which we see the characteristic upward trend indicative, physiologically, of a rate code (with $A$ now identified with $I_{\mathrm{app}}$).

pylab.rc('figure', figsize=(4.5,3.5))
pylab.rc('xtick', labelsize='small')
pylab.rc('ytick', labelsize='small')
pylab.figure()
pylab.subplot(111)
pylab.plot(A, omega)
pylab.xlabel(r'Applied current ($I_{\mathrm{app}}$) [$\mu$A/cm$^{2}$]', size='small')
pylab.ylabel(r'Spiking frequency ($\omega$) [Hz]', size='small')
pylab.subplots_adjust(left=0.15, bottom=0.2)
pylab.savefig('sage.png')
pylab.close()
pylab.rcdefaults()

Problem 3

Convert the HH model into "phasic" mode. By "phasic" I mean that the point-neuron does not fire repetitively for any $I_{\mathrm{app}}$ values—only one to a few spikes at onset and then it settles to a stable steady state. Many neurons in the auditory system behave phasically. Do this by, say, sliding some channel-gating dynamics along the $V$-axis (probably just for $I_{\mathrm{K}}$). [If you slide $x_{\infty} (V)$, then you must also slide $\tau_{x} (V)$.] If it can be done using $h = 1 - n$ and $m = m_{\infty} (V)$, then do the phase plane analysis.

 

Solution

We begin by implementing modified versions of the HH models in Problem 1, taking as an additional argument the voltage shift $V_{\mathrm{shift}}$, manifested in the dynamics of the potassium activation gating variable $n$ (to modify $I_{\mathrm{K}}$) as $$\dot{n} = \alpha \left( V - V_{\mathrm{shift}} \right) \left( 1 - n \right) - \beta \left( V - V_{\mathrm{shift}} \right) n,$$ i.e., we slide $n_{\infty}$ and $\tau_{n}$.

%cython

import numpy
cimport numpy
from numpy.linalg import norm

cdouble = numpy.double
ctypedef numpy.double_t cdouble_t

cpdef tuple m_dynamics(double V):
    cdef double alpha = ((V + 45)/10) / (1 - numpy.exp(-(V + 45)/10))
    cdef double beta = 4*numpy.exp(-(V + 70)/18)
    cdef double m_inf = alpha / (alpha + beta)
    cdef double tau_m = 1 / (alpha + beta)
    return m_inf, tau_m

cpdef tuple h_dynamics(double V):
    cdef double alpha = 0.07*numpy.exp(-(V + 70)/20)
    cdef double beta = 1 / (1 + numpy.exp(-(V + 40)/10))
    cdef double h_inf = alpha / (alpha + beta)
    cdef double tau_h = 1 / (alpha + beta)
    return h_inf, tau_h

cpdef tuple n_dynamics(double V):
    cdef double alpha = 0.1*((V + 60)/10) / (1 - numpy.exp(-(V + 60)/10))
    cdef double beta = 0.125*numpy.exp(-(V + 70)/80)
    cdef double n_inf = alpha / (alpha + beta)
    cdef double tau_n = 1 / (alpha + beta)
    return n_inf, tau_n

cpdef double I_inst(double V, double m, double h, double n, \
                    double g_Na, double g_K, double g_L,    \
                    double E_Na, double E_K, double E_L):
    return g_Na*m**3*h*(V - E_Na) + g_K*n**4*(V - E_K) + g_L*(V - E_L)

cpdef numpy.ndarray hh_Kmod_ode(double t, numpy.ndarray y,                               \
                                double C, double g_Na, double g_K, double g_L,           \
                                double E_Na, double E_K, double E_L, I, double V_shift):
    cdef double V = y[0]
    cdef double m = y[1]
    cdef double h = y[2]
    cdef double n = y[3]
    cdef double m_inf, tau_m, h_inf, tau_h, n_inf, tau_n
    m_inf, tau_m = m_dynamics(V)
    h_inf, tau_h = h_dynamics(V)
    n_inf, tau_n = n_dynamics(V - V_shift)
    return numpy.array([(I(t) - I_inst(V, m, h, n, g_Na, g_K, g_L, E_Na, E_K, E_L))/C, \
                        (m_inf - m)/tau_m,                                             \
                        (h_inf - h)/tau_h,                                             \
                        (n_inf - n)/tau_n])

cpdef numpy.ndarray hh2_Kmod_ode(double t, numpy.ndarray y,                               \
                                 double C, double g_Na, double g_K, double g_L,           \
                                 double E_Na, double E_K, double E_L, I, double V_shift):
    cdef double V = y[0]
    cdef double n = y[1]
    cdef double m = m_dynamics(V)[0]
    cdef double h = 1 - n
    cdef double n_inf, tau_n
    n_inf, tau_n = n_dynamics(V - V_shift)
    return numpy.array([(I(t) - I_inst(V, m, h, n, g_Na, g_K, g_L, E_Na, E_K, E_L))/C, \
                        (n_inf - n)/tau_n])

cpdef dt_rel(deriv, double t, numpy.ndarray y, double dt_min, double dt_max, \
             double C, double g_Na, double g_K, double g_L,                  \
             double E_Na, double E_K, double E_L, I, double V_shift):
    cdef numpy.ndarray[cdouble_t] dy =                             \
        deriv(t, y, C, g_Na, g_K, g_L, E_Na, E_K, E_L, I, V_shift)
    cdef double a = norm(dy / y)
    return dt_min*(a/(1 + a)) + dt_max*(1/(1 + a))

Numerical simulation then shows that taking $V_{\mathrm{shift}} = 50$ mV suffices for phasic conversion.

tspan = [0, 100]
dt_min = 1e-2
dt_max = 1e-1

A = 20
tau = 100
V_shift = 50

t, y = hh(hh_Kmod_ode, y0, tspan, dt_min, dt_max, \
          C, g_Na, g_K, g_L, E_Na, E_K, E_L,      \
          lambda t: I(t, A, tau), V_shift)

pylab.rc('figure', figsize=(5.5,3))
pylab.rc('legend', fontsize='x-small')
pylab.rc('xtick', labelsize='x-small')
pylab.rc('ytick', labelsize='x-small')
pylab.figure()
pylab.subplot(121)
pylab.plot(t, y[:,0])
pylab.xlim(tspan)
pylab.xlabel(r'Time ($t$) [ms]', size='small')
pylab.ylabel(r'Voltage ($V$) [mV]', size='small')
pylab.subplot(122)
pylab.plot(t, y[:,1], 'b', label=r'Na + ($m$)')
pylab.plot(t, y[:,2], 'g', label=r'Na - ($h$)')
pylab.plot(t, y[:,3], 'r', label=r'K + ($n$)')
pylab.xlim(tspan)
pylab.legend(loc='best')
pylab.xlabel(r'Time ($t$) [ms]', size='small')
pylab.ylabel(r'Gating variables', size='small')
pylab.subplots_adjust(bottom=0.2, wspace=0.3)
pylab.savefig('sage.png')
pylab.close()
pylab.rcdefaults()

In anticipation of the corresponding two-variable model, we implement the following code for plotting nullclines in the $n$-$V$ phase plane.

n_nc = (Iapp - I_sym).subs(h=1-n).solve(n)

def hh_Kmod_pplot(t, y, g_Na, g_K, g_L, E_Na, E_K, E_L, I, V_shift):
    def V_nc(V):
        return n_nc[1].rhs().subs({var('V'): V,                \
                                   m       : m_dynamics(V)[0], \
                                   gNa     : g_Na,             \
                                   gK      : g_K,              \
                                   gL      : g_L,              \
                                   ENa     : E_Na,             \
                                   EK      : E_K,              \
                                   EL      : E_L,              \
                                   Iapp    : I}).n()
    pylab.rc('figure', figsize=(5.5,3))
    pylab.rc('legend', fontsize='x-small')
    pylab.rc('xtick', labelsize='x-small')
    pylab.rc('ytick', labelsize='x-small')
    pylab.figure()
    pylab.subplot(121)
    pylab.plot(t, y[:,0])
    pylab.xlim(tspan)
    pylab.xlabel(r'Time ($t$) [ms]', size='small')
    pylab.ylabel(r'Voltage ($V$) [mV]', size='small')
    pylab.subplot(122)
    pylab.plot(y[:,0], y[:,3])
    V_lim = pylab.xlim()
    ylim = pylab.ylim()
    V_spn = pylab.linspace(*V_lim)
    pylab.plot(V_spn, [n_dynamics(V - V_shift)[0] for V in V_spn], label=r'$\dot{n} = 0$')
    pylab.plot(V_spn, [V_nc(V) for V in V_spn], label=r'$\dot{V} = 0$')
    pylab.legend(loc='best')
    pylab.ylim(ylim)
    pylab.xlabel(r'Voltage ($V$) [mV]', size='small')
    pylab.ylabel(r'K activation ($h$)', size='small')
    pylab.subplots_adjust(bottom=0.2, wspace=0.3)
    pylab.savefig('sage.png')
    pylab.close()
    pylab.rcdefaults()

There are four roots of $n$ to the equation $I_{\mathrm{inst}} = I_{\mathrm{app}}$ for $m = m_{\infty} (V)$ and $h = 1 - n$. To demonstrate that we have selected the correct root in the code above, we check which of the roots satisfies $0 \leq n \leq 1$ at standard resting values:

html(r'roots of $n(V_{\infty})$ for two-variable model:')
for s in n_nc:
    pretty_print(s.rhs().subs({V   : V_rest,                \
                               m   : m_dynamics(V_rest)[0], \
                               gNa : g_Na,                  \
                               gK  : g_K,                   \
                               gL  : g_L,                   \
                               ENa : E_Na,                  \
                               EK  : E_K,                   \
                               EL  : E_L,                   \
                               Iapp: 0}).n())

roots of n(V_{\infty}) for two-variable model:

\newcommand{\Bold}[1]{\mathbf{#1}}-0.344018153581324

\newcommand{\Bold}[1]{\mathbf{#1}}0.322725399339642

\newcommand{\Bold}[1]{\mathbf{#1}}0.0106463771208409 - 0.333711600086087i

\newcommand{\Bold}[1]{\mathbf{#1}}0.0106463771208409 + 0.333711600086087i

We now first study the two-variable model with $V_{\mathrm{shift}} = 0$ mV, necessary for proper interpretation by comparison later, with a sustained applied current of $A = 20$ $\mu$A/cm$^{2}$ and $\tau = 100$ ms.

tspan = [0, 100]
dt_min = 1e-2
dt_max = 1e-1

A = 20
tau = 100
V_shift = 0

t, y = hh(hh2_Kmod_ode, [y0[0], y0[3]], tspan, dt_min, dt_max, \
          C, g_Na, g_K, g_L, E_Na, E_K, E_L,                   \
          lambda t: I(t, A, tau), V_shift)
y = proc_hh2(y)
hh_Kmod_pplot(t, y, g_Na, g_K, g_L, E_Na, E_K, E_L, A, V_shift)

In this case, we obtain repetitive firing as demonstrated in Problem 1. The oscillating steady state indeed lies on the middle branch of the $V$-nullcline as suggested by Problem 2. In contrast, for $V_{\mathrm{shift}} = 50$ mV:

tspan = [0, 100]
dt_min = 1e-2
dt_max = 1e-1

A = 20
tau = 100
V_shift = 50

t, y = hh(hh2_Kmod_ode, [y0[0], y0[3]], tspan, dt_min, dt_max, \
          C, g_Na, g_K, g_L, E_Na, E_K, E_L,                   \
          lambda t: I(t, A, tau), V_shift)
y = proc_hh2(y)
hh_Kmod_pplot(t, y, g_Na, g_K, g_L, E_Na, E_K, E_L, A, V_shift)

we see that the $n$-nullcline has shifted rightward so that only one (stable) steady state now exists, which in particular does not lie on the middle branch. The trajectory hence relaxes to rest along the $n$-manifold, producing phasic firing.

We note also that phasic mode may be induced, for example, by modifying the conductances $g_{\mathrm{Na}}$, $g_{\mathrm{K}}$, and $g_{\mathrm{L}}$, as shown in Problem 2.