%pylab nbagg from tvb.simulator.lab import *
Populating the interactive namespace from numpy and matplotlib INFO log level set to INFO
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.
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()
C:\Users\mw\Downloads\TVB_Distribution\tvb_data\Lib\site-packages\matplotlib\cbook.py:137: MatplotlibDeprecationWarning: The set_color_cycle attribute was deprecated in version 1.5. Use set_prop_cycle instead. warnings.warn(message, mplDeprecation, stacklevel=1)
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$
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=-0.2, Kf=0.1, r=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:
Lets load the connectivity matrix and choose the coupling function
con = connectivity.Connectivity(load_default=True)
WARNING File 'hemispheres' not found in ZIP.
We choose a difference coupling function
coupl = coupling.Difference(a=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(load_default=True) #Initialise some Monitors with period in physical time mon_tavg = monitors.TemporalAverage(period=1.) mon_EEG = monitors.EEG(load_default=True, region_mapping=rm, period=1.) mon_SEEG = monitors.iEEG(load_default=True, region_mapping=rm, 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()
C:\Users\mw\Downloads\TVB_Distribution\tvb_data\Lib\site-packages\IPython\kernel\__main__.py:3: RuntimeWarning: invalid value encountered in divide app.launch_new_instance()
figure(figsize=(10,10)) plot(teeg[:], 10*eeg[:, 0, :, 0] + np.r_[:65]) yticks(np.r_[:75], mon_SEEG.sensors.labels[:65]) title("EEG")