# Izhikevich model

Computational neuroscience utilizes mathematical models to understand the complex dynamics of neuronal activity. Among various neuron models, the Izhikevich model stands out for its ability to combine biological fidelity with computational efficiency. Developed by Eugene Izhikevich in 2003ꜛ, this model simulates the spiking and bursting behavior of neurons with a remarkable balance between simplicity and biological relevance. In this post, we explore the properties of the Izhikevich model, examining its application and adaptability in simulating single neuron behaviors.

## Mathematical foundation

The Izhikevich model bridges the gap between detailed biophysical models like the Hodgkin-Huxley model and more abstract models like the Integrate-and-Fire model. While the former is biologically realistic but computationally expensive, the latter is computationally efficient but biologically unrealistic and lacks the ability to reproduce the rich dynamics of real neurons. In contrast, the Izhikevich model offers a compromise by capturing essential neuronal dynamics with a reduced set of equations. This reduction, based on bifurcation methodologies, allows for faster simulations while retaining essential features of neuronal dynamics.

The model is defined by a two-dimensional system of ordinary differential equations (ODE),

\[\begin{align} \frac{dv}{dt} &= 0.04v^2 + 5v + 140 - u + I \label{eq:model1} \\ \frac{du}{dt} &= a(bv - u) \label{eq:model2} \end{align}\]along with an after-spike reset condition:

\[\begin{align} \text{if } v \geq 30 \text{ mV, then } & \begin{cases} v \leftarrow c \\ u \leftarrow u + d \end{cases} \label{eq:reset} \end{align}\]The variable $v$ represents the **membrane potential**, and $u$ is a **recovery variable** that provides negative feedback to $v$. Both variables are dimensionless. $I$ is the synaptic or injected current. $a$, $b$, $c$, and $d$ are dimensionless parameters that define neuron dynamics, and they can be tuned to replicate various types of neuronal behaviors. In particular,

- $a$ controls the
**time scale of the recovery variable $u$**, with smaller values leading to slower recovery. - $b$ determines the
**sensitivity of the recovery variable $u$**to the subthreshold fluctuations of membrane potential $v$. Larger values result in a stronger coupling between $v$ and $u$, which possibly leads to subthreshold oscillations and low-threshold spiking behavior. - $c$ is the
**after-spike reset value of the membrane potential $v$**. - $d$ is the
**increment of the recovery variable $u$ after a spike**. It mimics the slow recovery of the membrane potential after an action potential.

The following plot visualizes the effect of the four parameters on the dynamics of the Izhikevich model:

The specific numbers in Eq. ($\ref{eq:model1}$), $0.04v^2 + 5v + 140$, result from fitting the spike generation dynamics to experimental data of cortical neurons so that the membrane potential $v$ has the scale of mV and the time $t$ has the scale of ms. Other choices of the parameters in Eq. ($\ref{eq:model1}$) are also possible to model different types of neurons.

Eq. ($\ref{eq:model1}$) and ($\ref{eq:model2}$) describe the evolution of the membrane potential and recovery variable over time. The reset condition Eq. ($\ref{eq:reset}$) ensures the model generates realistic action potentials by resetting the membrane potential and adjusting the recovery variable whenever the potential reaches a peak of 30 mV. The Izhikevich model is particularly well-suited for exploring the diverse behaviors of different types of neurons, including regular spiking, fast spiking, and bursting neurons as we will explore in the following sections.

## Python code

To simulate the model, we choose the Euler method for numerical integration similar to Izhikevich’s original paper. However, we modify the code in such a way that the integration is performed with an adjustable step size `dt`

. Choosing a smaller step size $\lt$1 will increase the accuracy of the simulation but also increase the computational cost. A value of `dt=0.1`

is a good starting point for most simulations with a good balance between accuracy and computational efficiency.

In the code, different sets of parameters are pre-defined for various neuron types, that we will discuss in more detail in the next section. By changing the parameter set, you can simulate different types of neurons. The code also allows you to adjust the input current `I_baseline`

over time as well as the start and end time of the input current.

For plotting the membrane potential $v$, we clip the values to a maximum of 30 mV to visualize realistic action potentials.

```
import numpy as np
import matplotlib.pyplot as plt
# set simulation time and time step size:
T = 400 # total simulation time in ms
dt = 0.1 # time step size in ms
steps = int(T / dt) # number of simulation steps
t_start = 50 # start time for the input current
t_end = T # end time for the input current
# initialize parameters for one excitatory neuron:
p_RS = [0.02, 0.2, -65, 8, "regular spiking (RS)"] # regular spiking settings for excitatory neurons (RS)
p_IB = [0.02, 0.2, -55, 4, "intrinsically bursting (IB)"] # intrinsically bursting (IB)
p_CH = [0.02, 0.2, -51, 2, "chattering (CH)"] # chattering (CH)
p_FS = [0.1, 0.2, -65, 2, "fast spiking (FS)"] # fast spiking (FS)
p_TC = [0.02, 0.25, -65, 0.05, "thalamic-cortical (TC)"] # thalamic-cortical (TC) (doesn't work well)
p_LTS = [0.02, 0.25, -65, 2, "low-threshold spiking (LTS)"] # low-threshold spiking (LTS)
p_RZ = [0.1, 0.26, -65, 2, "resonator (RZ)"] # resonator (RZ)
a, b, c, d, type = p_RS # just change the parameter set here to simulate different neuron types
# initial values of v and u:
v = -65 # mV
u = b * v
# initialize array to store the u, v, I and t values over time:
u_values = np.zeros(steps)
v_values = np.zeros(steps)
I_values = np.zeros(steps)
t_values = np.zeros(steps)
# set the baseline current:
I_baseline = 10 # nA
# simulation:
for t in range(steps):
t_ms = t * dt # current time in ms
if t_ms >= t_start and t_ms <= t_end:
I = I_baseline
else:
I = 0
# check for spike and reset if v >= 30 mV (reset-condition):
if v >= 30:
v = c # reset membrane potential v to c
u += d # increase recovery variable u by d
# Euler's method for numerical integration:
v += dt * 0.5 * (0.04 * v**2 + 5 * v + 140 - u + I)
v += dt * 0.5 * (0.04 * v**2 + 5 * v + 140 - u + I)
u += dt * a * (b * v - u)
# store values for plotting:
u_values[t] = u
v_values[t] = v
I_values[t] = I
t_values[t] = t_ms
# ensure v_values do not exceed 30 mV in the plot:
v_values = np.clip(v_values, None, 30)
# plotting:
fig, ax1 = plt.subplots(figsize=(8,3.85))
# plot v_values on the left y-axis:
ax1.plot(t_values, v_values, label='Membrane potential v(t)', color='k', lw=1.3)
ax1.set_xlabel('Time (ms)')
ax1.set_ylabel('membrane potential $v$ [mV]', color='k')
ax1.tick_params(axis='y', colors='k')
# create a second y-axis for u_values:
ax2 = ax1.twinx()
ax2.plot(t_values, u_values, label='Recovery variable u(t)', color='r', lw=2, alpha=1.0)
ax2.set_ylabel('recovery variable $u$ [a.u.]', color='r')
ax2.tick_params(axis='y', colors='r')
ax2.set_ylim(-20, 10)
# create a third y-axis for I_values:
ax3 = ax1.twinx()
ax3.spines['right'].set_position(('outward', 60))
ax3.plot(t_values, I_values, label='Input Current I(t)', color='b', lw=2, alpha=0.75)
ax3.set_ylabel('input current $I$ [nA]', color='b')
ax3.tick_params(axis='y', colors='b')
ax3.set_ylim(-1,60)
ax3.set_frame_on(True)
ax3.patch.set_visible(False)
plt.title(f'Membrane potential, recovery variable, and input current, {type}\n'
f"Parameters: a={a}, b={b}, c={c}, d={d}", fontsize=12)
plt.tight_layout()
plt.show()
```

The code above shows the core concept of the Izhikevich model simulation. The entire code can be found in this Github repositoryꜛ.

## Simulating different neuron types

The Izhikevich model can simulate various types of neurons by adjusting the parameters $a$, $b$, $c$, and $d$. Here are some examples of different neuron types and their corresponding parameter sets that can be found in the mammalian neocortex:

### Regular spiking neurons

The most common type of excitatory neuron in the neocortex are regular spiking (RS) neurons. They are characterized by a regular firing pattern when exposed to a constant input current, showing a short inter-spike interval at the beginning which increases over time. This behavior is also called **spike frequency adaptation**. An increase in the input current will increase the the inter-spike frequency. However, RS neurons will never fire at high frequencies due to the the large spike-afterhyperpolarization. The parameters for RS neurons are $a=0.02$, $b=0.2$, $c=-65$, and $d=8$. $c=-65$ corresponds to a large voltage reset after a spike, and $d=8$ corresponds to a large after-spike increase of the recovery variable $u$.

### Intrinsically bursting neurons

Intrinsically bursting (IB) neurons are characterized by a burst of action potentials followed by repetitive single spikes. The parameters for IB neurons are $a=0.02$, $b=0.2$, $c=-55$, and $d=4$. During the initial burst, the recovery variable $u$ increases rapidly and then switches to (regular) spiking dynamics.

### Chattering neurons

Chattering (CH) neurons are characterized by bursts of closely spaced action potentials followed by a hyperpolarization. The inter-burst frequency can reach high values up to 40 Hz. The parameters for CH neurons are $a=0.02$, $b=0.2$, $c=-51$ tp $-50$, and $d=2$. The lower value of $d$ compared to IB and RS neurons results in a slower recovery of the membrane potential after a burst.

### Fast spiking neurons

Fast spiking (FS) neurons are one of two types of inhibitory neurons in the cortex. They fire periodic trains of action potentials at high frequencies with almost no adaptation, i.e., no slowing down of the firing rate. The parameters for FS neurons are $a=0.1$, $b=0.2$, $c=-65$, and $d=2$. The higher value of $a$ compared to RS neurons results in a faster recovery of the membrane potential.

### Low-threshold spiking neurons

Low-threshold spiking (LTS) neurons are the second type of inhibitory neurons in the cortex. Similar to FS neurons, they fire periodic trains of action potentials at high frequencies. However, LTS neurons exhibit spike-frequency adaptation, leading to a decrease in the firing rate over time. These neurons also exhbit a low-threshold for spiking and can be simulated with the parameters $a=0.02$, $b=0.25$, $c=-65$, and $d=2$. $b=0.25$ accounts for the low-threshold spiking behavior.

### Thalamic-cortical neurons

The model is also able to simulate thalamic-cortical (TC) neurons, which are found in the thalamus and project to the cortex. TC neurons provide the major input to the cortex and are involved in the generation of sleep rhythms. They have two firing modes. First, when exposed to a constant positive input current, they exhibit a tonic firing pattern. Second, when exposed to a negative input current which abruptly switches to 0, the TC neurons fire a rebound of action potentials. The parameters for TC neurons are $a=0.02$, $b=0.25$, $c=-65$, and $d=0.05$.

### Resonator neurons

In his original paper, Izhikevich shows that the model is able to simulate another interesting type of neurons: Resonator (RZ) neurons. These neurons are able to resonate to rhythmic inputs that have an appropriate frequency. As far as I understand, the model mimics this behavior by exhibiting subthreshold oscillations in response to a constant input current. Furthermore, the neuron would switch to repetitive spiking when exposed to a short-term input current pulse. The corresponding parameters for RZ neurons are $a=0.1$, $b=0.26$, $c=-65$, and $d=2$. However, I was not able to reproduce the resonator behavior with these parameters. I tried different values and different input currents, without success. If you have any suggestions on how to simulate the resonator behavior, please let me in the comments below.

### Summary of parameters for different neuron types

The following plot summarizes the different parameter sets for the various neuron types described above:

## Conclusion

The Izhikevich model is a powerful tool for simulating the spiking and bursting behavior of neurons with a remarkable balance between simplicity and biological relevance. By adjusting the parameters $a$, $b$, $c$, and $d$, the model can simulate various types of neurons found in the mammalian neocortex, including regular spiking, fast spiking, and bursting neurons, while being computationally efficient.

Despite its advantages, the Izhikevich model is not without limitations. The model can oversimplify certain neuronal behaviors, particularly those involving complex ion channel dynamics and second messenger systems. Furthermore, the model’s reliance on parameter tuning for different neuron types can make it less predictive compared to more detailed models.

Nevertheless, the Izhikevich model serves as a bridge between biologically detailed, computationally demanding models and more abstract, simplified neuronal models. It provides a versatile platform for exploring neuronal behavior and network dynamics with considerable ease and efficiency. In the next post we will discover, how the Izhikevich model can be used to efficiently simulate networks of spiking neurons.

The complete code used in this blog post is available in this Github repositoryꜛ. Feel free to experiment with it, modify the parameters, and explore the dynamics of the Izhikevich model further. And for any ideas regarding the not yet solved resonator behavior, please leave a comment below.

## References and further reading

- Izhikevich,
*Simple model of spiking neurons*, 2003, IEEE Transactions on Neural Networks, Vol. 14, Issue 6, pages 1569-1572, doi: 10.1109/TNN.2003.820440ꜛ, PDFꜛ - Izhikevich, Eugene M., (2010),
*Dynamical systems in neuroscience: The geometry of excitability and bursting (First MIT Press paperback edition)*, The MIT Press, ISBN: 978-0-262-51420-0, PDFꜛ

## Comments

Comment on this post by publicly replying to this Mastodon post using a Mastodon or other ActivityPub/Fediverse account.

Comments on this website are based on a Mastodon-powered comment system. Learn more about it here.

There are no known comments, yet. Be the first to write a reply.