︠83f38885-a083-49e6-83d6-b5f6895b315fi︠ %html

Computational Modeling of Neuronal Systems

Homework 1

 

Kenneth L. Ho

ho@courant.nyu.edu

︡7427524b-2a99-4649-80af-058723c44578︡{"html": "

Computational Modeling of Neuronal Systems

\n

Homework 1

\n

 

\n

Kenneth L. Ho

\n

ho@courant.nyu.edu

"}︡ ︠cfaab037-6524-44da-93bb-6e68b0e098c3︠ version() ︡e1e56d57-a759-4214-a2b6-79eaab3ab4aa︡{"stdout": "'Sage Version 4.3.2, Release Date: 2010-02-06'"}︡ ︠b0c2cd7e-e652-4af5-bf46-184416daa42e︠ import pylab from scipy import fftpack from scipy import integrate from scipy import optimize ︡b46bd627-5a02-4943-b459-41babd87443f︡︡ ︠9d7685b5-d0fb-412c-9565-16433d801315i︠ %html

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 dt or tol 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 SingPt 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.

︡4f249b0b-47c0-4078-90ae-1deffdf3cdf7︡{"html": "

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 dt or tol for numerical integration so that tightening the criterion gives comparable results. State the error criteria and integration method that you use.

\n

 

\n

Problem 1

\n

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

\n
    \n
  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. \n
  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. \n
  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. \n
  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. \n
  5. If you use XPP, use the SingPt 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.
    6. \n
\n

 

\n

Solution

\n

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

"}︡ ︠09706331-f634-4cab-a5e9-9002e0a6342d︠ %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)) ︡48597572-d957-45b5-889f-8804298d6f8d︡︡ ︠eb5af20c-9efa-4c5d-ad1a-163252f7d3edi︠ %html

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

︡ea5c7cac-c089-4990-8d23-4d226f8f8aa7︡{"html": "

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

"}︡ ︠0baba189-f8ce-42ef-b940-a5fc2b95596a︠ 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 ︡dabb544c-a983-4eff-9891-59bac2fa516c︡︡ ︠05ce0246-55b2-40b9-90ec-bc2f523f9324i︠ %html

A note on the numerics: we use the solver vode, as wrapped in SciPy, for numerical integration. As vode 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}$$

︡2a544a90-afff-44a7-a0a3-e16068cf9a12︡{"html": "

A note on the numerics: we use the solver vode, as wrapped in SciPy, for numerical integration. As vode 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.

\n

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}$$

"}︡ ︠850977be-9f3f-4aaf-bf3e-6c23bc38ff67︠ 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 ︡590cf8e4-0a27-46a0-b4cd-4da6c41efbf6︡︡ ︠322444ad-ceaa-49c0-a7ea-687c9e6d9e3bi︠ %html

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

︡0c776af8-ef87-4ff3-a19a-891d29724e10︡{"html": "

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

"}︡ ︠3a22fc30-53e7-43c2-978d-fe34696087c0︠ 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) ︡96d1473c-76f8-4b3a-98f8-8ad30d269227︡︡ ︠6a372dff-ae21-4aad-bce3-101277c0f64ai︠ %html

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

︡24fa478a-71b6-40a8-b41e-822010a648fb︡{"html": "

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.

\n

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

"}︡ ︠7e5605a1-8d76-4e2c-b4a1-cbe820230528︠ 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) ︡3d332c1f-ff94-42f8-add4-f4de45fa746a︡︡ ︠369b583b-f53e-4843-8e10-308aeb10d269i︠ %html

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

︡094e73b6-c2b4-432b-aa70-db44a6d6d21c︡{"html": "

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.

\n

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

"}︡ ︠9786cb0a-6dbe-4a19-85fc-073b53391b39︠ 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) ︡addc0a7e-678e-41b0-bed2-ec41df39ac4e︡︡ ︠7930b01f-28f8-4c6a-9d61-078aa90a7832i︠ %html

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:

︡bcf74a19-f362-4424-9207-0e88560f983e︡{"html": "

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.

\n

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:

"}︡ ︠c62efebb-8875-46e5-a4f4-9fe0f21025bb︠ 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) ︡8526b2d6-2ad4-4c90-99b5-34c970b4b614︡︡ ︠3c2373d0-d7a1-4686-819c-6264bbb17412i︠ %html

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:

︡9d950add-433a-490f-961b-ac5e8328a358︡{"html": "

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.

\n

Likewise, the long stimulus gives:

"}︡ ︠f1a46b6a-a1eb-4dca-a167-8ebb485f2bc6︠ 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) ︡eb50a8ea-3cc4-43ab-be8b-f5aa49195d1a︡︡ ︠108b3d70-76e8-41b7-adc5-c17464caf88fi︠ %html

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.

︡775a0588-d793-4305-8574-ac0d837dd4ee︡{"html": "

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.

\n

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.

"}︡ ︠d55293dd-8931-42fd-bd9b-512b1f6b7be8︠ 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) ︡eab4e323-9c6c-418b-b7f4-deda32833919︡︡ ︠7c6c6d4b-0b39-4e48-abe9-58d9edc6de21i︠ %html

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.

︡f6783f4b-8cda-43f6-841e-b906b99bbf86︡{"html": "

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.

"}︡ ︠5daa2576-2b4b-4d52-bd03-e6a825efc8ce︠ 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) ︡404f437d-692c-4165-a2f4-b13969b8293e︡︡ ︠84c4902e-5b8e-4b22-9216-4d0b08ae6053i︠ %html

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.

︡3952ebec-20fd-4da1-87d2-3caf6d4bc61a︡{"html": "

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.

\n

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.

\n

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.

"}︡ ︠c8060223-0fa7-4ea6-8365-09bbb84349d0︠ 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) ︡fdb87130-8139-405f-9a24-08414149ba57︡︡ ︠2a2a9cfb-01b3-4aa0-9d1e-12ae29c7364di︠ %html

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:

︡f4baee61-31f5-4c8e-a6e6-006b02727018︡{"html": "

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:

"}︡ ︠786f6498-c15f-4d41-9c9e-3927f5b7f131︠ 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 ︡3e5298df-38a5-4ba7-ac4a-79a28d30f351︡{"html": "full model:"}︡{"html": "(V_{\\infty}, m_{\\infty}, h_{\\infty}, n_{\\infty}) at I_{\\mathrm{app}} = 0:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\hbox{-69.8976728964}, 0.0535746092323, 0.592537659007, 0.319246167222\\right)"}︡{"html": "eigenvalues of J at I_{\\mathrm{app}} = 0:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\hbox{[-4.66655603+0.j -0.20044489+0.3876745j -0.20044489-0.3876745j\n -0.12082016+0.j ]}"}︡{"html": "(V_{\\infty}, m_{\\infty}, h_{\\infty}, n_{\\infty}) at I_{\\mathrm{app}} = 20:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\hbox{-61.5643150194}, 0.13497383385, 0.306805365841, 0.450925338254\\right)"}︡{"html": "eigenvalues of J at I_{\\mathrm{app}} = 20:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\hbox{[-5.28245933+0.j 0.15640342+0.64195044j 0.15640342-0.64195044j\n -0.15780123+0.j ]}"}︡ ︠1828c6ea-d8f8-4c85-8be6-447142d26840i︠ %html

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:

︡c9e862f8-3d4f-474d-b2c5-5b267e4897ac︡{"html": "

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:

"}︡ ︠853b0e78-6fae-47d2-a189-46216de2ecbc︠ 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 ︡b41bdc50-dda4-4d6c-8bf5-965cce8d46e9︡{"html": "three-variable model:"}︡{"html": "(V_{\\infty}, m_{\\infty}, h_{\\infty}, n_{\\infty}) at I_{\\mathrm{app}} = 0:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\hbox{-69.8976728964}, 0.0535746092323, 0.592537659007, 0.319246167222\\right)"}︡{"html": "eigenvalues of J at I_{\\mathrm{app}} = 0:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\hbox{[-0.21057224+0.40888856j -0.21057224-0.40888856j -0.12080474+0.j ]}"}︡{"html": "(V_{\\infty}, m_{\\infty}, h_{\\infty}, n_{\\infty}) at I_{\\mathrm{app}} = 20:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\hbox{-61.5643150194}, 0.13497383385, 0.306805365841, 0.450925338254\\right)"}︡{"html": "eigenvalues of J at I_{\\mathrm{app}} = 20:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\hbox{[ 0.44882718+0.77185623j 0.44882718-0.77185623j -0.15773290+0.j ]}"}︡ ︠85b73135-bb39-42a3-a8a6-eb1c128dc826i︠ %html

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,

︡97688a58-1d4e-4cc6-a621-66387de1756a︡{"html": "

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.

\n

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

"}︡ ︠24712e28-bfd1-43ce-916b-71b3d158e4b9︠ 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 ︡67d1a6e3-f348-4262-8dd9-f318abc8622d︡{"html": "two-variable model near V_{\\infty} \\approx -70 mV:"}︡{"html": "(V_{\\infty}, m_{\\infty}, h_{\\infty}, n_{\\infty}) at I_{\\mathrm{app}} = 0:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\hbox{-69.7262906627}, 0.0546657944747, 0.678120453084, 0.321879546916\\right)"}︡{"html": "eigenvalues of J at I_{\\mathrm{app}} = 0:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\hbox{[-0.17284854+0.41570498j -0.17284854-0.41570498j]}"}︡{"stdout": "Warning: The iteration is not making good progress, as measured by the \n improvement from the last ten iterations."}︡{"html": "(V_{\\infty}, m_{\\infty}, h_{\\infty}, n_{\\infty}) at I_{\\mathrm{app}} = 20:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\hbox{-58.9351936326}, 0.175377559045, 0.508015695861, 0.491984304139\\right)"}︡{"html": "eigenvalues of J at I_{\\mathrm{app}} = 20:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\hbox{[ 6.92443132e+00 -3.60566769e-04]}"}︡ ︠4aba287f-f3de-4d4b-8754-887313ff2280i︠ %html

second:

︡08f9cd98-390f-45ba-a6fd-fdd65dc1ecee︡{"html": "

second:

"}︡ ︠9064ffcd-c7de-4014-a1ca-c74bd6a2c352︠ 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 ︡ad4b79aa-63c3-4b3b-86ee-8e493270c398︡{"html": "two-variable model near V_{\\infty} \\approx -50 mV:"}︡{"html": "(V_{\\infty}, m_{\\infty}, h_{\\infty}, n_{\\infty}) at I_{\\mathrm{app}} = 0:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\hbox{-54.0694295866}, 0.27115411034, 0.435878786546, 0.564121213454\\right)"}︡{"html": "eigenvalues of J at I_{\\mathrm{app}} = 0:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\hbox{[ 20.38459968 -0.08022857]}"}︡{"stdout": "Warning: The iteration is not making good progress, as measured by the \n improvement from the last five Jacobian evaluations."}︡{"html": "(V_{\\infty}, m_{\\infty}, h_{\\infty}, n_{\\infty}) at I_{\\mathrm{app}} = 20:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\hbox{-58.9502980712}, 0.175122534275, 0.508248314858, 0.491751685142\\right)"}︡{"html": "eigenvalues of J at I_{\\mathrm{app}} = 20:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\hbox{[ 6.89629867e+00 1.62786029e-05]}"}︡ ︠25062fd3-3e63-4f5a-8f7a-f2e70489f08fi︠ %html

and third:

︡449d1e00-71b8-4876-b630-cc1c5f453285︡{"html": "

and third:

"}︡ ︠cc783e24-9587-4383-bf36-cee7d075bc82︠ 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 ︡5fbb1060-d550-49db-89ee-751349543eba︡{"html": "two-variable model near V_{\\infty} \\approx -25 mV:"}︡{"html": "(V_{\\infty}, m_{\\infty}, h_{\\infty}, n_{\\infty}) at I_{\\mathrm{app}} = 0:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\hbox{-26.5156451255}, 0.859989720941, 0.17297674804, 0.82702325196\\right)"}︡{"html": "eigenvalues of J at I_{\\mathrm{app}} = 0:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\hbox{[ 2.65557338+3.70630697j 2.65557338-3.70630697j]}"}︡{"html": "(V_{\\infty}, m_{\\infty}, h_{\\infty}, n_{\\infty}) at I_{\\mathrm{app}} = 20:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\hbox{-26.1174195632}, 0.864288169081, 0.170787639316, 0.829212360684\\right)"}︡{"html": "eigenvalues of J at I_{\\mathrm{app}} = 20:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\hbox{[ 1.88345651+4.23974391j 1.88345651-4.23974391j]}"}︡ ︠e3f6041e-93fd-4507-baf9-350068c6ffafi︠ %html

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.

︡b9f4636f-b556-4042-9829-f5c72562aa9c︡{"html": "

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.

\n

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.

"}︡ ︠142c6ec6-a931-4847-af3d-931d57391991︠ 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() ︡fe476609-0816-4294-aa23-acbbdad1f3a5︡︡ ︠f8797363-eb81-4c65-8066-913a381099cbi︠ %html

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.

︡2b1e736a-82f0-4a54-a7f2-2561c5a5cedd︡{"html": "

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.

\n

 

\n

Problem 2

\n

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.

\n

 

\n

Solution

\n

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.

"}︡ ︠6707f698-4f3c-43be-a149-35963729a853︠ 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() ︡c089e728-8bb9-4a5b-89d3-395e6aedc81f︡︡ ︠90d14ccd-1206-4ea0-b84e-90514f47df96i︠ %html

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}}$.

︡abf4b909-dd50-44b9-a000-f2aadb89486f︡{"html": "

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}}$.

"}︡ ︠d8cf486d-3cc8-4e61-aff9-94e66653761e︠ 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) ︡7d735ff9-550e-4e65-b9e2-9aab6a5452a1︡︡ ︠8f5b3112-9d8b-48a7-8009-db6dafa2ff8ci︠ %html

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:

︡894e3190-a808-4df1-8098-a72ffed5404a︡{"html": "

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:

"}︡ ︠c97af2cb-c07c-4b49-b8f8-76909ef4bbad︠ 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) ︡2247e0b5-0d6a-42ec-ab1d-31bb9a6139c5︡︡ ︠499c0c7a-02cc-4bfe-8156-2449a6412023i︠ %html

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

︡139d47aa-70a4-4d21-a37c-936c6b716206︡{"html": "

\n

\n

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.

\n

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.

\n

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

\n
\n

"}︡ ︠f3520d97-a440-48da-9b55-e1e22524b08f︠ 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 ︡fee6f700-91e5-4871-81d7-1aabf68927ed︡︡ ︠ff69c36e-49d9-468b-b7af-6fa771e63928i︠ %html

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.

︡61ca3840-a311-49ba-9d69-d6ea2d951b99︡{"html": "

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.

"}︡ ︠2d7cf766-3d94-4c20-96dd-4b64a209858c︠ 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)) ︡e4ea697f-0d2d-4b55-8a68-d865ccff149a︡︡ ︠d95aa617-2d5a-411d-9e1d-21c68336a5e1i︠ %html

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}}$).

︡932cfad8-c0f4-44a9-9fce-0eadc33b38a0︡{"html": "

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}}$).

"}︡ ︠3603e236-2565-4062-bd5f-60327a65f9b9︠ 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() ︡476012b3-ff16-4750-9b7b-66727d36e90d︡︡ ︠df856da4-729a-44ea-919e-dd1db2b63e1ci︠ %html

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}$.

︡fbd48c0a-afb0-47ef-b7fb-15fbb1028e1c︡{"html": "

Problem 3

\n

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.

\n

 

\n

Solution

\n

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}$.

"}︡ ︠4235a6c2-f3d8-4ba2-b121-6e6df535f811︠ %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)) ︡4417d61c-534b-42e7-8226-44de791f0e18︡︡ ︠86e1c8f7-cd09-4785-bc5e-f37fd40ae237i︠ %html

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

︡700c07a7-4255-4e82-ac4e-9e09991471ba︡{"html": "

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

"}︡ ︠9c556706-524d-490a-ac84-0d8abf4b8fb2︠ 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() ︡0dedb172-401c-4abd-a2a4-47e0b62fb015︡︡ ︠ef73773c-c835-40b8-9414-0b916a062690i︠ %html

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

︡91d2e72b-f604-4d2c-90e2-613f51e46f40︡{"html": "

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

"}︡ ︠f0d4c96a-80a3-4258-8b0f-c8012d7aa2c0︠ 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() ︡aeee3742-76d4-4a27-ad69-49f5c9fad765︡︡ ︠ee040959-05de-4cac-a5f0-9d58e20f80b3i︠ %html

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:

︡6623a8fb-357f-4d2c-a824-6f9fb913a133︡{"html": "

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:

"}︡ ︠60b6d528-6b31-4e73-8f4c-dca0f78b6ce0︠ 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()) ︡45a63efc-cc61-4f9d-8e76-2be13ba2bc5b︡{"html": "roots of n(V_{\\infty}) for two-variable model:"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}-0.344018153581324"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}0.322725399339642"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}0.0106463771208409 - 0.333711600086087i"}︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}0.0106463771208409 + 0.333711600086087i"}︡ ︠cfd60f8b-96a5-4f49-a62f-4133678e1315i︠ %html

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.

︡a65bd40c-9dcc-4ab4-bc83-4730f3f81340︡{"html": "

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.

"}︡ ︠83befe64-2a4a-4ab5-ba59-079c8f87df92︠ 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) ︡cb65d927-b447-4fdf-ac24-3915a8b634ce︡︡ ︠16f92bec-8f99-4041-ad9f-e29a9245a23ei︠ %html

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:

︡4c4951b3-6bd0-4767-aa02-5ce1c267d316︡{"html": "

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:

"}︡ ︠89d606a3-8de9-4ae0-8d3f-a17d98db8f4a︠ 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) ︡9f296aa5-896d-49f0-880f-b0929946be3e︡︡ ︠e7d0ccd9-9df7-48ad-9e94-fb8b9f6bec44i︠ %html

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.

︡3b596186-0241-4cac-bc5a-4ee5cd93e3b6︡{"html": "

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.

\n

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.

"}︡