EpileptorStep

EpileptorStep#

class brainmass.EpileptorStep(in_size, a=1.0, b=3.0, c=1.0, d=5.0, r=0.00035, x0=-1.6, Iext=3.1, slope=0.0, Iext2=0.45, tau=10.0, aa=6.0, bb=2.0, Kvf=0.0, Kf=0.0, Ks=0.0, tt=1.0, modification=0.0, init_x1=Constant(value=-1.5), init_y1=Constant(value=-10.0), init_z=Constant(value=3.5), init_x2=Constant(value=-1.0), init_y2=Constant(value=0.0), init_g=Constant(value=0.0), noise_x1=None, noise_x2=None, method='exp_euler')#

Epileptor model of seizure dynamics (Jirsa et al., 2014).

The Hindmarsh-Rose-Jirsa Epileptor is a six-dimensional neural-mass model that reproduces the full taxonomy of epileptic seizures: the autonomous transition from interictal (between-seizure) to ictal (seizure) activity and back, driven by a slow permittivity variable [1] [2]. It couples three subsystems on separate time scales:

  • a fast subsystem \((x_1, y_1)\) generating ictal fast discharges,

  • an ultra-slow subsystem \((x_2, y_2)\) for spike-and-wave events,

  • a slow permittivity variable \(z\) that ramps seizures on and off,

plus a low-pass filter \(g\) linking the two subsystems.

\[\begin{split}\begin{aligned} \dot x_1 &= t_t\,(y_1 - z + I_{\mathrm{ext}} + K_{vf} c_1 + f_1(x_1, x_2)\,x_1), \\ \dot y_1 &= t_t\,(c - d\,x_1^2 - y_1), \\ \dot z &= t_t\,r\,(h(x_1, z) - z + K_s c_1), \\ \dot x_2 &= t_t\,(-y_2 + x_2 - x_2^3 + I_{\mathrm{ext2}} + b_b\,g - 0.3(z - 3.5) + K_f c_2), \\ \dot y_2 &= t_t\,(-y_2 + f_2(x_2))/\tau, \\ \dot g &= t_t\,(-0.01\,(g - 0.1\,x_1)), \end{aligned}\end{split}\]

with the piecewise nonlinearities

\[\begin{split}f_1 = \begin{cases} -a x_1^2 + b x_1 & x_1 < 0 \\ \mathrm{slope} - x_2 + 0.6 (z - 4)^2 & x_1 \ge 0 \end{cases}, \quad f_2 = \begin{cases} 0 & x_2 < -0.25 \\ a_a (x_2 + 0.25) & x_2 \ge -0.25 \end{cases},\end{split}\]

and the permittivity coupling \(h\), a blend (via modification \(\in [0, 1]\)) of a linear and a sigmoidal form,

\[\begin{split}h = \mathrm{mod}\,\bigl(x_0 + \tfrac{3}{1 + e^{-(x_1 + 0.5)/0.1}}\bigr) + (1 - \mathrm{mod})\,\bigl(4(x_1 - x_0) + z_{\mathrm{nl}}\bigr), \quad z_{\mathrm{nl}} = \begin{cases} -0.1 z^7 & z < 0 \\ 0 & z \ge 0 \end{cases}.\end{split}\]

Here \(c_1\) and \(c_2\) are external/coupling inputs to populations 1 and 2 (x1_inp and x2_inp); \(c_1\) additionally drives the permittivity through \(K_s\).

Parameters:

  • in_size (int | Sequence[int] | integer | Sequence[integer]) – Spatial shape of the population. All parameters broadcast to this shape.

  • a (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Quadratic / linear coefficients of \(f_1\). Defaults 1.0, 3.0.

  • b (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Quadratic / linear coefficients of \(f_1\). Defaults 1.0, 3.0.

  • c (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Additive / quadratic coefficients of \(\dot y_1\). Defaults 1.0, 5.0.

  • d (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Additive / quadratic coefficients of \(\dot y_1\). Defaults 1.0, 5.0.

  • r (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Permittivity rate (inverse slow time constant). Default 0.00035 — the seizure-cycle time scale; larger values speed seizures up.

  • x0 (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Epileptogenicity parameter (seizure threshold). Default -1.6 (an epileptogenic node). More negative values render the node healthy.

  • Iext (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – External input to population 1. Default 3.1.

  • slope (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Linear coefficient in the \(x_1 \ge 0\) branch of \(f_1\). Default 0.0.

  • Iext2 (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – External input to population 2. Default 0.45.

  • tau (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Time scale of \(y_2\). Default 10.0.

  • aa (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Slope of \(f_2\). Default 6.0.

  • bb (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Coupling from the filter \(g\) into \(x_2\). Default 2.0.

  • Kvf (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Very-fast (\(x_1\)), fast (\(x_2\)) and slow (\(z\)) coupling scales. Defaults 0.0.

  • Kf (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Very-fast (\(x_1\)), fast (\(x_2\)) and slow (\(z\)) coupling scales. Defaults 0.0.

  • Ks (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Very-fast (\(x_1\)), fast (\(x_2\)) and slow (\(z\)) coupling scales. Defaults 0.0.

  • tt (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Global time-scale factor. Default 1.0.

  • modification (Callable | Array | ndarray | bool | number | bool | int | float | complex | Quantity | Param) – Blend in [0, 1] selecting the nonlinear (1) vs linear (0) permittivity influence on \(z\). Default 0.0.

  • init_x1 (Callable) – State initializers. Defaults reproduce the canonical initial condition (-1.5, -10.0, 3.5, -1.0, 0.0, 0.0).

  • init_y1 (Callable) – State initializers. Defaults reproduce the canonical initial condition (-1.5, -10.0, 3.5, -1.0, 0.0, 0.0).

  • init_z (Callable) – State initializers. Defaults reproduce the canonical initial condition (-1.5, -10.0, 3.5, -1.0, 0.0, 0.0).

  • init_x2 (Callable) – State initializers. Defaults reproduce the canonical initial condition (-1.5, -10.0, 3.5, -1.0, 0.0, 0.0).

  • init_y2 (Callable) – State initializers. Defaults reproduce the canonical initial condition (-1.5, -10.0, 3.5, -1.0, 0.0, 0.0).

  • init_g (Callable) – State initializers. Defaults reproduce the canonical initial condition (-1.5, -10.0, 3.5, -1.0, 0.0, 0.0).

  • noise_x1 (Noise) – Additive noise processes for populations 1 and 2. Default None.

  • noise_x2 (Noise) – Additive noise processes for populations 1 and 2. Default None.

  • method (str) – Integration method, 'exp_euler' (default; well-suited to the stiff slow z) or any braintools.quad method (e.g. 'rk4').

Return type:

Any

x1, y1

Fast subsystem (membrane-like activity and recovery).

Type:

brainstate.HiddenState

z#

Slow permittivity variable driving seizure onset/offset.

Type:

brainstate.HiddenState

x2, y2

Ultra-slow subsystem (spike-and-wave).

Type:

brainstate.HiddenState

g#

Low-pass filter of x1.

Type:

brainstate.HiddenState

Notes

  • State variables are dimensionless; every right-hand side carries unit 1/ms so an exponential-Euler step with dt in milliseconds is consistent.

  • The z time scale is stiff and slow (r = 3.5e-4): a full interictal-ictal cycle spans \(\mathcal{O}(1/r)\) time units. Raising r rescales time and makes seizures observable over fewer steps.

  • lfp() returns the standard local-field-potential proxy \(x_2 - x_1\), which update() returns.

References

Examples

>>> import brainmass
>>> import brainstate
>>> import brainunit as u
>>> model = brainmass.EpileptorStep(in_size=1)
>>> _ = brainstate.nn.init_all_states(model)
>>> with brainstate.environ.context(dt=0.1 * u.ms):
...     lfp = model.update()
>>> lfp.shape
(1,)
__init__(in_size, a=1.0, b=3.0, c=1.0, d=5.0, r=0.00035, x0=-1.6, Iext=3.1, slope=0.0, Iext2=0.45, tau=10.0, aa=6.0, bb=2.0, Kvf=0.0, Kf=0.0, Ks=0.0, tt=1.0, modification=0.0, init_x1=Constant(value=-1.5), init_y1=Constant(value=-10.0), init_z=Constant(value=3.5), init_x2=Constant(value=-1.0), init_y2=Constant(value=0.0), init_g=Constant(value=0.0), noise_x1=None, noise_x2=None, method='exp_euler')[source]#
Parameters:
derivative(state, t, x1_inp, x2_inp, add1=0.0, add2=0.0)[source]#
dg(g, x1)[source]#

Right-hand side for the low-pass filter g (unit 1/ms).

dx1(x1, y1, z, x2, x1_inp, add=0.0)[source]#

Right-hand side for the fast activity x1 (unit 1/ms).

x1_inp is the population-1 coupling current (scaled by Kvf); add is a direct additive perturbation (e.g. noise) applied at the derivative level, independent of the coupling gain.

dx2(x2, y2, z, g, x2_inp, add=0.0)[source]#

Right-hand side for the ultra-slow activity x2 (unit 1/ms).

x2_inp is the population-2 coupling current (scaled by Kf); add is a direct additive perturbation (e.g. noise) applied at the derivative level, independent of the coupling gain.

dy1(y1, x1)[source]#

Right-hand side for the fast recovery y1 (unit 1/ms).

dy2(y2, x2)[source]#

Right-hand side for the ultra-slow recovery y2 (unit 1/ms).

dz(z, x1, x1_inp)[source]#

Right-hand side for the slow permittivity z (unit 1/ms).

init_state(batch_size=None, **kwargs)[source]#

Allocate the six Epileptor states at their canonical initial values.

lfp()[source]#

Local-field-potential proxy x2 - x1 (the standard Epileptor output).

update(x1_inp=None, x2_inp=None)[source]#

Advance the six Epileptor states by one time step.

Parameters:

  • x1_inp (array-like or scalar or None, optional) – Coupling input to population 1 (also drives z through Ks). If None, treated as zero; noise_x1 is added when set.

  • x2_inp (array-like or scalar or None, optional) – Coupling input to population 2. If None, treated as zero; noise_x2 is added when set.

Returns:

The local-field-potential proxy lfp() (x2 - x1).

Return type:

array-like