%matplotlib widget
import matplotlib.pyplot as plt
import numpy as np
import time as tm
from tvb.simulator.lab import *
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}$.
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):
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.
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:
# 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()
# 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']
# 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()
print("Starting simulation...")
tic = tm.time()
(t, y), = sim.run(simulation_length=10000)
print("Finished simulation.")
print('execute for ' + str(tm.time()-tic))
# 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])
# 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()
print("Starting simulation...")
tic = tm.time()
(s, z), = sim.run(simulation_length=10000)
print("Finished simulation.")
print('execute for ' + str(tm.time()-tic))
# 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)
# 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.
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).
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()
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:
# 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])
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()
print("Starting simulation...")
tic = tm.time()
(s_, z_), = sim.run(simulation_length=10000)
print("Finished simulation.")
print('execute for' + str(tm.time()-tic))
# Normalize the time series to have nice plots
z_ /= (np.max(z_,0) - np.min(z_,0 ))
z_ -= np.mean(z_, 0)
# 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()