In [1]:
%matplotlib widget

import matplotlib.pyplot as plt
import numpy as np
import time as tm

from tvb.simulator.lab import *


Reduction 2-dimentional Epileptor model¶

The Epileptor 2D model is the reduction of the 5 dimensions Epileptor model to only 2 dimensions, by applying averaging methods (see Proix et al., 2014). Then, the Epileptor equations become :

\begin{eqnarray} \dot{x}_{1,i} &=& - x_{1,i}^{3} - 2x_{1,i}^{2} + 1 + I_{ext_{1}}\\ \dot{z}_{i} &=& \dfrac{1}{\tau_{0}} \ (4(x_{1,i} - x_{0, i}) - z_{i} - K{s}\sum_{j=1}^{N}C_{ij}(x_{1,j} - x_{1,i}))) \end{eqnarray}

where $\tau_{0}=2857$, $I_{ext_{1}}=3.1$, and the degree of epileptogenicity is represented through the value $x_{0}$. $C_{ij}$ are the entries of the anatomical structural connectivity matrix reweigthed by the global coupling parameter $K_{s}$.

1. Exploring the Epileptor 2D model¶

Before launching any simulations, we will have a look at the phase space of the model in order to better understand its dynamics. We will use the phase plane interactive tool.

We plot the two nullclines that are defined as zero flux in either the $x_{1}$ or $z$ direction: the linear nullcline in green and the cubic nullcline in red. The intersections of the nullclines identify the fixed points of the system. The interictal state and the ictal state correspond to the left and the right branches of the cubic nullcline, respectively.

Have a look to the phase space for the 2D system (Fig. 8C in Proix et al. 2014):

In [2]:
from tvb.simulator.plot.phase_plane_interactive import PhasePlaneInteractive

# Create an Epileptor model instance.
epi = models.Epileptor2D(x0=np.array([-1]))

# Initialise a stochastic Integrator scheme.
dt = 0.5                                           #integration steps [ms]
heundetint = integrators.HeunDeterministic(dt=dt)
#heunstochint = integrators.HeunStochastic(dt=dt)
#heunstochint.noise.nsig = 0.01                      #standard deviation of the noise

# Open the phase plane tool with the Epileptor model and (stochastic) Integrator.
ppi_fig = PhasePlaneInteractive(model=epi, integrator=heundetint)
ppi_fig.show()


According to the value of the parameter $\pmb{x_{0}}$, the linear nullcline (green curve) moves left and right, changing the stability of the fixed point. Two typical trajectories are possible:

• For $\pmb{x_{0} < -2.05}$, the trajectory is attracted to the stable fixed point in the interictal state on the left cubic nullcline and the Epileptor is said not epileptogenic, meaning not triggering epileptic seizure without external input.

• For $\pmb{x_{0} > -2.05}$, the stable point disappears through a Hopf bifurcation, and there is only one unstable fixed point left. The Epileptor enters an oscillatory regime and the Epileptor is said to be epileptogenic, triggering seizures autonomously.

2. Comparison with complete 5D Epileptor¶

In order to compare the results obtained from the complete and reduced system, we will model a patient with a temporal lobe epilepsy (TLE).

To this end, we define a spatial map of epileptogenicity where each network's node $i$ is characterized by an excitability value $x_{0,i}$, which quantifies its ability to trigger a seizure.

We set:

• the right limbic areas (right hippocampus (rHC, region 47), parahippocampus (rPHC, region 62) and amygdala (rAMYG, region 40)) as Epileptic Zone (EZ) with an epileptogenicity parameter value equal to - 1.6,
• and two lesser epileptogenic regions: the inferior temporal cortex (rTCI, region 69) and the ventral temporal cortex (rTCV, region 72) as Propagation Zone (PZ), with an epileptogenicity parameter value equal to - 1.8,
• and all the other regions (or Non-Epileptogenic Zone, NEZ) are set to -2.4.
In [4]:
# Initialise a Connectivity.
con = connectivity.Connectivity.from_file()
N = con.weights.shape[0]
con.weights - con.weights * np.eye(N, N)
con.weights = con.weights / np.abs(con.weights.max())
con.tract_lengths = np.zeros((con.tract_lengths.shape))             # no time-delays
con.configure()

WARNING  File 'hemispheres' not found in ZIP.

In [8]:
# Initialise the complete model.
Epileptor5D = models.epileptor.Epileptor(Ks=np.array([1]), r=np.array([0.00015]))
Epileptor5D.x0 = np.ones((76))*-2.4                #NEZ
Epileptor5D.x0[[62, 47, 40]] = np.ones((3))*-1.6   #EZ
Epileptor5D.x0[[69, 72]] = np.ones((2))*-1.8       #PZ

#Initial conditions.
Epileptor5D.state_variable_range["x1"] = np.array([-1.8, -1.8])
Epileptor5D.state_variable_range["y1"] = np.array([-15, -15])
Epileptor5D.state_variable_range["x2"] = np.array([-1, -1])
Epileptor5D.state_variable_range["y2"] = np.array([0.01, 0.01])
Epileptor5D.state_variable_range["z"] = np.array([3.6, 3.6])

Epileptor5D.variables_of_interest = ['x2 - x1', 'x1', 'z']

In [10]:
# Initialise Simulator.
sim = simulator.Simulator(model=Epileptor5D,
connectivity=con,
coupling=coupling.Difference(a=np.array([-0.2])),
integrator=integrators.HeunDeterministic(dt=0.05),
monitors=(monitors.TemporalAverage(period=1.),))

sim.configure()

Out[10]:

Simulator

value
Type
Simulator
conduction_speed
3.0
connectivity
Connectivity gid: aabe80d5-6c97-4b40-a40f-0a2c6d9037a6
coupling
Difference gid: 89dba22a-e018-4f9d-868e-77653fd74f40
gid
UUID('ae4b6454-5c93-4e44-97b7-7855bd866f40')
initial_conditions
None
integrator
HeunDeterministic gid: ea689eef-cc91-4dbe-aa69-966eb7be6249
model
Epileptor gid: 427ca1ad-0821-4b0c-b9c3-87de9092d873
monitors
(,)
simulation_length
1000.0
stimulus
None
surface
None
title
Simulator gid: ae4b6454-5c93-4e44-97b7-7855bd866f40
In [11]:
print("Starting simulation...")
tic = tm.time()
(t, y), = sim.run(simulation_length=10000)
print("Finished simulation.")
print('execute for ' + str(tm.time()-tic))

Starting simulation...
WARNING  random_state supplied for non-stochastic integration
Finished simulation.
execute for 51.79368233680725

In [13]:
# Initialise the reduced model.
Epileptor2D = models.epileptor.Epileptor2D(Ks=np.array([1]), r=np.array([0.00015]))
Epileptor2D.x0 = np.ones((76))*-2.4               #NEZ
Epileptor2D.x0[[62, 47, 40]] = np.ones((3))*-1.6  #EZ
Epileptor2D.x0[[69, 72]] = np.ones((2))*-1.8      #PZ

#Initial conditions.
Epileptor2D.state_variable_range["x1"] = np.array([-1.8, -1.8])
Epileptor2D.state_variable_range["z"] = np.array([3.6, 3.6])

In [15]:
# Initialise Simulator.
sim = simulator.Simulator(model=Epileptor2D,
connectivity=con,
coupling=coupling.Difference(a=np.array([-0.2])),
integrator=integrators.HeunDeterministic(dt=0.05),
monitors=(monitors.TemporalAverage(period=1.),))

sim.configure()

Out[15]:

Simulator

value
Type
Simulator
conduction_speed
3.0
connectivity
Connectivity gid: aabe80d5-6c97-4b40-a40f-0a2c6d9037a6
coupling
Difference gid: ac6aa2b5-9134-435e-a473-7dcce7dcf0a1
gid
UUID('dd4b53ae-0e74-4de6-927f-c5dd6b871df1')
initial_conditions
None
integrator
HeunDeterministic gid: 6b35a241-fce4-49ac-a478-9b778e7008ed
model
Epileptor2D gid: 4d80a2f9-144a-4e69-abb7-d3a1bf1162aa
monitors
(,)
simulation_length
1000.0
stimulus
None
surface
None
title
Simulator gid: dd4b53ae-0e74-4de6-927f-c5dd6b871df1
In [16]:
print("Starting simulation...")
tic = tm.time()
(s, z), = sim.run(simulation_length=10000)
print("Finished simulation.")
print('execute for ' + str(tm.time()-tic))

Starting simulation...
WARNING  random_state supplied for non-stochastic integration
Finished simulation.
execute for 42.40776085853577

In [17]:
# Normalize the time series to have nice plots
y /= (np.max(y, 0) - np.min(y, 0))
y -= np.mean(y, 0)
z /= (np.max(z, 0) - np.min(z, 0))
z -= np.mean(z, 0)

In [18]:
# Plot time series
plt.figure(figsize=(10,10))
plt.plot(t[:], y[:, 1, :, 0] + np.r_[:76], 'C3')
plt.plot(s[:], z[:, 0, :, 0] + np.r_[:76], 'C0')
plt.title("Epileptors 5D (red) vs. 2D (blue) time series", fontsize=15)
plt.xlabel('Time [ms]', fontsize=15)
plt.yticks(np.arange(len(sim.connectivity.region_labels)), sim.connectivity.region_labels, fontsize=6)
plt.show()


Both systems show good correspondance with slight differences in the intrinsic frequencies.

3. Modification of the slow permittivity variable $z$ dynamics¶

The duration of the ictal and interictal state is approximately the same in the Epileptor, which is not the case in clinical situations, where the interictal period is typically longer. To this end, we use a slight modification of the Epileptor equations and replace the linear influence of $x_{1}$ on the slow permittivity variable $z$ by the nonlinear function $\pmb{h}$ causing a symmetry breaking between ictal and interictal period with an increase of the latter. The Epileptor equations with the here used modification read then as follows:

\begin{eqnarray} \dot{x}_{1,i} &=& - x_{1,i}^{3} - 2x_{1,i}^{2} + 1 + I_{ext_{1}}\\ \dot{z}_{i} &=& \dfrac{1}{\tau_{0}} \ (h(x_{1,i}) - z_{i} - K{s}\sum_{j=1}^{N}C_{ij}(x_{1,j} - x_{1,i}))) \end{eqnarray}

where

\begin{eqnarray} h(x_{1,i}) = x_{0,i} + 3. \ / \left(1 + \exp\left(\dfrac{-x_{1,i} - 0.5}{0.1}\right)\right) \end{eqnarray}

First, we have a look to the phase space of the 2D system for the modified Epileptor equations (see Fig. 8B in Proix et al., 2014).

In [20]:
from tvb.simulator.plot.phase_plane_interactive import PhasePlaneInteractive

# Create an Epileptor model instance.
epim = models.epileptor.Epileptor2D(x0=np.array([2]), modification=np.array([True]))
# Initialise a stochastic Integrator scheme.
dt = 0.75                                           #integration steps [ms]
heundetint = integrators.HeunDeterministic(dt=dt)
#heunstochint = integrators.HeunStochastic(dt=dt)
#heunstochint.noise.nsig = 0.01                      #standard deviation of the noise

# Open the phase plane tool with the Epileptor model and (stochastic) Integrator.
ppi_fig = PhasePlaneInteractive(model=epim, integrator=heundetint, exclude_sliders=["x0"])
ppi_fig.show()

c:\work\tvb\tvb-root\scientific_library\tvb\simulator\plot\phase_plane_interactive.py:325: UserWarning: Attempting to set identical left == right == -1.8 results in singular transformations; automatically expanding.
valinit=self.model.state_variable_range[self.svx][0])
c:\work\tvb\tvb-root\scientific_library\tvb\simulator\plot\phase_plane_interactive.py:330: UserWarning: Attempting to set identical left == right == -1.8 results in singular transformations; automatically expanding.
valinit=self.model.state_variable_range[self.svx][1])
c:\work\tvb\tvb-root\scientific_library\tvb\simulator\plot\phase_plane_interactive.py:336: UserWarning: Attempting to set identical left == right == 3.6 results in singular transformations; automatically expanding.
valinit=self.model.state_variable_range[self.svy][0])
c:\work\tvb\tvb-root\scientific_library\tvb\simulator\plot\phase_plane_interactive.py:341: UserWarning: Attempting to set identical left == right == 3.6 results in singular transformations; automatically expanding.
valinit=self.model.state_variable_range[self.svy][1])
c:\work\tvb\tvb-root\scientific_library\tvb\simulator\plot\phase_plane_interactive.py:365: UserWarning: Attempting to set identical left == right == -1.8 results in singular transformations; automatically expanding.
valinit=self.default_sv[sv, 0, 0])
c:\work\tvb\tvb-root\scientific_library\tvb\simulator\plot\phase_plane_interactive.py:365: UserWarning: Attempting to set identical left == right == 3.6 results in singular transformations; automatically expanding.
valinit=self.default_sv[sv, 0, 0])
c:\work\tvb\tvb-root\scientific_library\tvb\simulator\plot\phase_plane_interactive.py:742: UserWarning: No contour levels were found within the data range.
[0], colors="r")
c:\work\tvb\tvb-root\scientific_library\tvb\simulator\plot\phase_plane_interactive.py:745: UserWarning: No contour levels were found within the data range.
[0], colors="g")


According to the value of the parameter $\pmb{x_{0}}$, the sigmoid nullcline (green curve) moves up and down, changing the number and stability of the fixed points. Two typical trajectories are possible:

• For $\pmb{x_{0} > 2.91}$, the trajectory is attracted to the stable fixed point in the interictal state on the left cubic nullcline and the Epileptor is said not epileptogenic, meaning not triggering epileptic seizure without external input.
• For $\pmb{x_{0} < 2.91}$, the two left most fixed points disappear through a SNIC bifurcation, and there is only one unstable fixed point left. The Epileptor enters an oscillatory regime and the Epileptor is said to be epileptogenic, triggering seizures autonomously.
In [23]:
# Initialise the Model.
Epileptormod = models.epileptor.Epileptor2D(Ks=np.array([1.]), r=np.array([0.00025]), modification=np.array([True]))
Epileptormod.x0 = np.ones((76)) * 3.5
Epileptormod.x0[[62, 47, 40]] = np.ones((3)) * 2.
Epileptormod.x0[[69, 72]] = np.ones((2)) * 2.3

Epileptormod.state_variable_range["x1"] = np.array([-1.8, -1.8])
Epileptormod.state_variable_range["z"] = np.array([3.6, 3.6])

In [25]:
sim = simulator.Simulator(model=Epileptormod,
connectivity=con,
coupling=coupling.Difference(a=np.array([-0.2])),
integrator=integrators.HeunDeterministic(dt=0.05),
monitors=(monitors.TemporalAverage(period=1.),))

sim.configure()

Out[25]:

Simulator

value
Type
Simulator
conduction_speed
3.0
connectivity
Connectivity gid: aabe80d5-6c97-4b40-a40f-0a2c6d9037a6
coupling
Difference gid: 46be1ba0-faea-4f81-ac79-c90ddc8008d9
gid
UUID('4d90e73e-8af5-40d4-a130-e1b2381d2e21')
initial_conditions
None
integrator
HeunDeterministic gid: 755f1337-5d54-4d7f-9490-f4afebca2fd9
model
Epileptor2D gid: 8ee9e158-e1fc-4101-a530-1164d9e35024
monitors
(,)
simulation_length
1000.0
stimulus
None
surface
None
title
Simulator gid: 4d90e73e-8af5-40d4-a130-e1b2381d2e21
In [26]:
print("Starting simulation...")
tic = tm.time()
(s_, z_),  = sim.run(simulation_length=10000)
print("Finished simulation.")
print('execute for' + str(tm.time()-tic))

Starting simulation...
WARNING  random_state supplied for non-stochastic integration
Finished simulation.
execute for39.08166813850403

In [27]:
# Normalize the time series to have nice plots
z_ /= (np.max(z_,0) - np.min(z_,0 ))
z_ -= np.mean(z_, 0)

In [28]:
# Plot time series
plt.figure(figsize=(10,10))
#plot(s[:], z[:, 0, :, 0] + np.r_[:76], 'C4')
plt.plot(s_[:], z_[:, 0, :, 0] + np.r_[:76], 'C0')
plt.title("Modified Epileptors time series", fontsize=15)
plt.xlabel('Time [ms]', fontsize=15)
plt.yticks(np.arange(len(sim.connectivity.region_labels)), sim.connectivity.region_labels, fontsize=6)
plt.show()