Simulation of the Izhikevich model

Izhikevic developed a two-dimensional model that includes a variety of response patterns seen in real neurons. The subthreshold phase is modelled by the following system of nonlinear differential equations:

$$\begin{aligned}\frac{{dV}}{{dt}} &= 0.04\,{V^2} + 5\,V + 140 - U + {I_{{\rm{stim}}}}\\ \frac{{dU}}{{dt}} &= a \cdot (b \cdot V - U)\end{aligned}$$ where $V$ the membrane potential, $U$ is the recovery variable, and $I_\mathrm{stim}$ is the stimulus. The parameter $a$ describes the time scale of the recovery variable and the parameter $b$ describes the senitivity of the recovery variable to the subthreshold fluctuations of the membrane potential. When the membrane potential $V$ reaches a threshold value it is reset to the value $c$ and the recovery variable $U$ is increased with the value of $d$: $$V(t + \delta ) \rightarrow c\quad {\rm{and}}\quad U(t + \delta ) \rightarrow U + d\quad {\rm{when}}\quad \delta \rightarrow 0$$