%pylab nbagg
from tvb.simulator.lab import *
Modeling Epilepsy¶
The main goal of this tutorial is to provide a clear understanding of how we can reproduce clinically relevant senarios such as the propagation of an epileptic seizure in the human brain, electrical stimulation of a brain region that can trigger a seizure, or surgical resection of brain regions.
Exploring the Epileptor model¶
The Epileptor is a phenomenological neural mass model able to reproduce epileptic seizure dynamics such as recorded with intracranial EEG electrodes (see Jirsa_et_al). Before launching any simulations, we will have a look at the phase space of the Epileptor model to better understand its dynamics. We will use the phase plane interactive tool.
from tvb.simulator.plot.phase_plane_interactive import PhasePlaneInteractive
# Create an Epileptor model instance
epileptor = models.Epileptor()
# Open the phase plane tool with the epileptor model
ppi_fig = PhasePlaneInteractive(model=epileptor)
ppi_fig.show()
Looking at the phase space, we have here the first population (variables $y_0$ in abscissa and $y_1$ in ordinate). The left most intersection of the nullcline defines a stable fixed point, representing the interictal state, whereas the rightmost intersection is the center of a limit cycle, being the ictal state. Both states are separated by a separatrix, as you can see by drawing different trajectories in this phase space (left click on the figure).
You can also look at other variables in the phase space, such as the second population $y_3$ & $y_4$, responsible for hte interical spikes in the Epileptor model. Change the lower and upper bound of the axis to see correctly the phase space.
You can continue to play along to explore the dynamics of this model. For instance, try changing the parameter $x_0$
Region based simulation of a temporal lobe seizure¶
We will model a patient with temporal lobe epilepsy (TLE). We will set different values of epileptogenicity $x_0$ parameter in the Epileptor according to the region positions, thereby introducing heterogeneities in the network parameters. We set the right limbic areas (right hippocampus (rHC, region 47), parahippocampus (rPHC, region 62) and amygdala (rAMYG, region 40)) as epileptic zones. We also add two lesser epileptogenic regions: the superior temporal cortex (rTCI, region 69) and the ventral temporal cortex (rTCV, region 72).
In other words, we assign to all the nodes the Dynamics for which $x_0$ has a value of value of $-2.2$. We apply the epileptogenic configuration ($-1.6$) to the right limbic areas.
Additionally, we chose which kind of coupling we want (between the fast variable (Kvf), the spike-and-wave events (Kf), or the slow permittive coupling (Ks)). Here we use Kf and Ks of them.
Finally, we also slow-down the dynamics of the Epileptor by choosing r=0.00015
epileptors = models.Epileptor(Ks=numpy.array([-0.2]), Kf=numpy.array([0.1]), r=numpy.array([0.00015]))
epileptors.x0 = np.ones((76))*-2.4
epileptors.x0[[62, 47, 40]] = np.ones((3))*-1.6
epileptors.x0[[69, 72]] = np.ones((2))*-1.8
All the other model parameters are the default ones:
| Model parameter | Value |
|---|---|
| $Iext$ | 3.1 |
| $Iext2$ | 0.45 |
| $slope$ | 0.0 |
Lets load the connectivity matrix and choose the coupling function
con = connectivity.Connectivity.from_file()
We choose a difference coupling function
coupl = coupling.Difference(a=numpy.array([1.]))
We use a stochastic integration scheme; the noise is only added on the two variables of the second population (y3, y4)
hiss = noise.Additive(nsig = numpy.array([0., 0., 0., 0.0003, 0.0003, 0.]))
heunint = integrators.HeunStochastic(dt=0.05, noise=hiss)
We will now set the monitors: a temporal average, an EEG and a SEEG. We need for this to load a region mapping, the projection matrices and the sensors.
In the Epileptor model, the LFP is define as -y0+y3. We want the projection matrices to be applied on the LFP, so we use this as a 'pre' expression. We also keep track of the slow permittivity variable y2.
# load the default region mapping
rm = region_mapping.RegionMapping.from_file()
#Initialise some Monitors with period in physical time
mon_tavg = monitors.TemporalAverage(period=1.)
mon_EEG = monitors.EEG.from_file()
mon_EEG.region_mapping=rm
mon_EEG.period=1.
mon_SEEG = monitors.iEEG.from_file()
mon_SEEG.region_mapping=rm
mon_SEEG.period=1.
#Bundle them
what_to_watch = (mon_tavg, mon_EEG, mon_SEEG)
Finally, we iniatilise and configure our Simulator object
#Initialise a Simulator -- Model, Connectivity, Integrator, and Monitors.
sim = simulator.Simulator(model=epileptors, connectivity=con,
coupling=coupl,
integrator=heunint, monitors=what_to_watch)
sim.configure()
We perform the simulation of 10.000 ms
(ttavg, tavg), (teeg, eeg), (tseeg, seeg) = sim.run(simulation_length=10000)
And we plot the results
# Normalize the time series to have nice plots
tavg /= (np.max(tavg,0) - np.min(tavg,0 ))
eeg /= (np.max(eeg,0) - np.min(eeg,0 ))
eeg -= np.mean(eeg, 0)
seeg /= (np.max(seeg,0) - np.min(seeg, 0))
seeg -= np.mean(seeg, 0)
#Plot raw time series
figure(figsize=(10,10))
plot(ttavg[:], tavg[:, 0, :, 0] + np.r_[:76], 'r')
title("Epileptors time series")
show()
figure(figsize=(10,10))
plot(teeg[:], 10*eeg[:, 0, :, 0] + np.r_[:65])
yticks(np.r_[:65], mon_EEG.sensors.labels)
title("EEG")
figure(figsize=(10,10))
plot(tseeg[:], seeg[:, 0, :75, 0] + np.r_[:75])
yticks(np.r_[:75], mon_SEEG.sensors.labels)
title("SEEG")
Modeling surgical resection¶
Surgical resection is used for around 20% of epileptic patient whose seizures are drug- resistant. We will simulate the hypothetic case of a surgical resection of the amygdala and the hippocampus, but leaving the parahippocampal cortex.
We set all the connections to the right amygdala (40) and right hippocampus (47) to 0 in the connectivity matrix. The resection of the EZ is not complete, will it be enough to prevent seizures?
con = connectivity.Connectivity.from_file()
con.weights[[ 47, 40]] = 0.
con.weights[:, [47, 40]] = 0.
# we plot only the right hemisphere
# the lines and columns set to 0 are clearly visible
figure()
imshow(con.weights[38:, 38:], interpolation='nearest', cmap='binary')
colorbar()
show()
The rest of the model is set as before, but we just use a time average monitor:
coupl = coupling.Difference(a=numpy.array([1.]))
hiss = noise.Additive(nsig = numpy.array([0., 0., 0., 0.0003, 0.0003, 0.]))
heunint = integrators.HeunStochastic(dt=0.05, noise=hiss)
mon_tavg = monitors.TemporalAverage(period=1.)
We can now relaunch our first simulation, taking care of replacing the dynamic of the EZ by a stable node, as if the region was resected.
epileptors = models.Epileptor(Ks=numpy.array([-0.2]), Kf=numpy.array([0.1]), r=numpy.array([0.00015]))
epileptors.x0 = np.ones((76))*-2.4
epileptors.x0[[69, 72]] = np.ones((2))*-1.8
sim = simulator.Simulator(model=epileptors, connectivity=con,
coupling=coupl,
integrator=heunint, monitors=(mon_tavg,))
sim.configure();
(ttavg, tavg), = sim.run(simulation_length=10000)
As you can see, no seizure is triggered anymore
# normalize the time series
tavg /= (np.max(tavg,0) - np.min(tavg,0 ))
figure(figsize=(10,10))
plot(ttavg[:], tavg[:, 0, :, 0] + np.r_[:76], 'r')
title("Epileptors time series")
show()
Triggering a seizure by stimulation¶
We are now going to model an electric stimulation and trigger a seizure. We set the whole brain to non-epileptogenic but close to the threshold
epileptors = models.Epileptor(Ks=numpy.array([-0.2]), Kf=numpy.array([0.1]), r=numpy.array([0.00035]))
epileptors.x0 = np.ones((76))*-2.1
con = connectivity.Connectivity.from_file()
coupl = coupling.Difference(a=numpy.array([0.]))
heunint = integrators.HeunDeterministic(dt=0.05)
mon_tavg = monitors.TemporalAverage(period=1.)
# Weighting for regions to receive stimuli.
weighting = numpy.zeros((76))
weighting[[69, 72]] = numpy.array([2.])
eqn_t = equations.PulseTrain()
eqn_t.parameters["T"] = 10000.0
eqn_t.parameters["onset"] = 2500.0
eqn_t.parameters["tau"] = 50.0
eqn_t.parameters["amp"] = 20.0
stimulus = patterns.StimuliRegion(temporal = eqn_t,
connectivity = con,
weight = weighting)
#Configure space and time
stimulus.configure_space()
stimulus.configure_time(numpy.arange(0., 3000., heunint.dt))
#And take a look
plot_pattern(stimulus)
show()
#Bundle them
what_to_watch = (mon_tavg, )
#Initialise Simulator -- Model, Connectivity, Integrator, Monitors, and stimulus.
sim = simulator.Simulator(model=epileptors, connectivity=con,
coupling=coupl,
integrator=heunint, monitors=what_to_watch,
stimulus=stimulus)
sim.configure()
(t, tavg), = sim.run(simulation_length=10000)
# Normalize the time series to have nice plots
tavg /= (np.max(tavg,0) - np.min(tavg,0 ))
#Plot raw time series
figure(figsize=(10,10))
plot(t[:], tavg[:, 0, :, 0] + np.r_[:76], 'r')
title("Epileptors time series")
#Show them
show()
