# NaK model analysis

Here we will show you the neurodynamics analysis of a two-dimensional system model with the example of the \(I_{\rm{Na,p+}}-I_K\) Model.

The dynamical system is given by:

where

This model specifies a leak current \(I_L\), persistent sodium current \(I_{\rm{Na, p}}\) with instantaneous activation kinetic, and a relatively slower persistent potassium current \(I_K\) with either high or low threshold (the two choices result in fundamentally different dynamics).

```
[1]:
```

```
from collections import OrderedDict
import brainpy as bp
bp.backend.set(dt=0.02)
```

```
[2]:
```

```
C = 1
E_L = -78 # different from high-threshold model
g_L = 8
g_Na = 20
g_K = 10
E_K = -90
E_Na = 60
Vm_half = -20
k_m = 15
Vn_half = -45 # different from high-threshold model
k_n = 5
tau = 1
@bp.odeint
def integral(V, n, t, Iext):
m_inf = 1 / (1 + bp.ops.exp((Vm_half - V) / k_m))
I_leak = g_L * (V - E_L)
I_Na = g_Na * m_inf * (V - E_Na)
I_K = g_K * n * (V - E_K)
dvdt = (-I_leak - I_Na - I_K + Iext) / C
n_inf = 1 / (1 + bp.ops.exp((Vn_half - V) / k_n))
dndt = (n_inf - n) / tau
return dvdt, dndt
```

## Phase plane analysis

```
[3]:
```

```
analyzer = bp.analysis.PhasePlane(
integral,
target_vars=OrderedDict(V=[-90, 20], n=[0., 1.]),
pars_update={'Iext': 50.})
analyzer.plot_nullcline()
analyzer.plot_vector_field()
analyzer.plot_fixed_point()
analyzer.plot_trajectory([{'V': -10, 'n': 0.2}, {'V': -80, 'n': 0.4}],
duration=100.,
show=True)
```

```
plot nullcline ...
SymPy solve "integral(V, n) = 0" to "n = f(V, )", success.
SymPy solve "integral(V, n) = 0" to "n = f(V, )", success.
plot vector field ...
```

```
:2: RuntimeWarning: divide by zero encountered in true_divide
```

```
plot fixed point ...
SymPy solve derivative of "integral(V, n)" by "V", success.
SymPy solve derivative of "integral(V, n)" by "n", success.
SymPy solve derivative of "integral(V, n)" by "V", success.
SymPy solve derivative of "integral(V, n)" by "n", success.
Fixed point #1 at V=-51.60868767199546, n=0.2105293582848025 is a unstable focus.
plot trajectory ...
```

## Codimension 1 bifurcation analysis

Here we show the codimension 1 bifurcation analysis of the \(I_{\rm{Na,p+}}-I_K\) Model, in which \(I_{ext}\) is varied in [0., 50.].

```
[4]:
```

```
analyzer = bp.analysis.Bifurcation(
integral,
target_pars={'Iext': [0, 50.]},
target_vars={"V": [-90., 20.], 'n': [0., 1.]},
numerical_resolution=0.1)
analyzer.plot_bifurcation(show=True)
```

```
plot bifurcation ...
SymPy solve "integral(V, n, Iext) = 0" to "n = f(V, Iext)", success.
SymPy solve derivative of "integral(V, n, Iext)" by "V", success.
SymPy solve derivative of "integral(V, n, Iext)" by "n", success.
SymPy solve derivative of "integral(V, n, Iext)" by "V", success.
SymPy solve derivative of "integral(V, n, Iext)" by "n", success.
```

## Codimension 2 bifurcation analysis

Codimension 2 bifurcation analysis of the \(I_{\rm{Na,p+}}-I_K\) Model, in which \(I_{ext}\) is varied in [0., 50.], and “Vn_half” is varied in [-50, -40].

```
[5]:
```

```
analyzer = bp.analysis.Bifurcation(
integral,
target_pars={'Iext': [0, 50.], 'Vn_half': [-50, -40]},
target_vars={"V": [-90., 20.], 'n': [0., 1.]},
numerical_resolution=0.1)
analyzer.plot_bifurcation(show=True)
```

```
plot bifurcation ...
SymPy solve "integral(V, n, Iext, Vn_half) = 0" to "n = f(V, Iext,Vn_half)", success.
SymPy solve derivative of "integral(V, n, Iext, Vn_half)" by "V", success.
SymPy solve derivative of "integral(V, n, Iext, Vn_half)" by "n", success.
SymPy solve derivative of "integral(V, n, Iext, Vn_half)" by "V", success.
SymPy solve derivative of "integral(V, n, Iext, Vn_half)" by "n", success.
```

## Reference

[1] Izhikevich, Eugene M. Dynamical systems in neuroscience (Chapter 4). MIT press, 2007.