In [1]:
%pylab nbagg
from tvb.simulator.lab import *
from tvb.simulator.plot.phase_plane_interactive import PhasePlaneInteractive
from tvb.simulator.models.epileptorcodim3 import EpileptorCodim3
from tvb.simulator.models.epileptorcodim3 import EpileptorCodim3SlowMod
from mpl_toolkits.mplot3d import Axes3D
Populating the interactive namespace from numpy and matplotlib
   INFO  log level set to INFO
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.wilson_cowan.WilsonCowan.state_variable_range = Final(field_type=<class 'dict'>, default={'E': array([0., 1.]), 'I': array([0., 1.])}, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.stefanescu_jirsa.ReducedSetFitzHughNagumo.state_variable_range = Final(field_type=<class 'dict'>, default={'xi': array([-4.,  4.]), 'eta': array([-3.,  3.]), 'alpha': array([-4.,  4.]), 'beta': array([-3.,  3.])}, required=True)
WARNING  default contains values out of the declared domain. Ex 1.0 
   attribute  tvb.simulator.models.stefanescu_jirsa.ReducedSetHindmarshRose.a = NArray(label=':math:`a`', dtype=float64, default=array([1.]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 3.0 
   attribute  tvb.simulator.models.stefanescu_jirsa.ReducedSetHindmarshRose.b = NArray(label=':math:`b`', dtype=float64, default=array([3.]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 1.0 
   attribute  tvb.simulator.models.stefanescu_jirsa.ReducedSetHindmarshRose.c = NArray(label=':math:`c`', dtype=float64, default=array([1.]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 3.3 
   attribute  tvb.simulator.models.stefanescu_jirsa.ReducedSetHindmarshRose.mu = NArray(label=':math:`\\mu`', dtype=float64, default=array([3.3]), dim_names=(), ndim=None, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.stefanescu_jirsa.ReducedSetHindmarshRose.state_variable_range = Final(field_type=<class 'dict'>, default={'xi': array([-4.,  4.]), 'eta': array([-25.,  20.]), 'tau': array([ 2., 10.]), 'alpha': array([-4.,  4.]), 'beta': array([-20.,  20.]), 'gamma': array([ 2., 10.])}, required=True)
WARNING  default contains values out of the declared domain. Ex 0.12 
   attribute  tvb.simulator.models.jansen_rit.JansenRit.p_min = NArray(label=':math:`p_{min}`', dtype=float64, default=array([0.12]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 0.32 
   attribute  tvb.simulator.models.jansen_rit.JansenRit.p_max = NArray(label=':math:`p_{max}`', dtype=float64, default=array([0.32]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 0.22 
   attribute  tvb.simulator.models.jansen_rit.JansenRit.mu = NArray(label=':math:`\\mu_{max}`', dtype=float64, default=array([0.22]), dim_names=(), ndim=None, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.jansen_rit.JansenRit.state_variable_range = Final(field_type=<class 'dict'>, default={'y0': array([-1.,  1.]), 'y1': array([-500.,  500.]), 'y2': array([-50.,  50.]), 'y3': array([-6.,  6.]), 'y4': array([-20.,  20.]), 'y5': array([-500.,  500.])}, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.jansen_rit.ZetterbergJansen.state_variable_range = Final(field_type=<class 'dict'>, default={'v1': array([-100.,  100.]), 'y1': array([-500.,  500.]), 'v2': array([-100.,   50.]), 'y2': array([-100.,    6.]), 'v3': array([-100.,    6.]), 'y3': array([-100.,    6.]), 'v4': array([-100.,   20.]), 'y4': array([-100.,   20.]), 'v5': array([-100.,   20.]), 'y5': array([-500.,  500.]), 'v6': array([-100.,   20.]), 'v7': array([-100.,   20.])}, required=True)
WARNING  default contains values out of the declared domain. Ex 1.0 
   attribute  tvb.simulator.models.oscillator.Generic2dOscillator.gamma = NArray(label=':math:`\\gamma`', dtype=float64, default=array([1.]), dim_names=(), ndim=None, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.oscillator.Generic2dOscillator.state_variable_range = Final(field_type=<class 'dict'>, default={'V': array([-2.,  4.]), 'W': array([-6.,  6.])}, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.oscillator.Kuramoto.state_variable_range = Final(field_type=<class 'dict'>, default={'theta': array([0.        , 6.28318531])}, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.oscillator.supHopf.state_variable_range = Final(field_type=<class 'dict'>, default={'x': array([-5.,  5.]), 'y': array([-5.,  5.])}, required=True)
WARNING  default contains values out of the declared domain. Ex -0.01 
   attribute  tvb.simulator.models.larter_breakspear.LarterBreakspear.TCa = NArray(label=':math:`T_{Ca}`', dtype=float64, default=array([-0.01]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 0.3 
   attribute  tvb.simulator.models.larter_breakspear.LarterBreakspear.TNa = NArray(label=':math:`T_{Na}`', dtype=float64, default=array([0.3]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 2.0 
   attribute  tvb.simulator.models.larter_breakspear.LarterBreakspear.aei = NArray(label=':math:`a_{ei}`', dtype=float64, default=array([2.]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 2.0 
   attribute  tvb.simulator.models.larter_breakspear.LarterBreakspear.aie = NArray(label=':math:`a_{ie}`', dtype=float64, default=array([2.]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 1.0 
   attribute  tvb.simulator.models.larter_breakspear.LarterBreakspear.ane = NArray(label=':math:`a_{ne}`', dtype=float64, default=array([1.]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 0.3 
   attribute  tvb.simulator.models.larter_breakspear.LarterBreakspear.Iext = NArray(label=':math:`I_{ext}`', dtype=float64, default=array([0.3]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 1.0 
   attribute  tvb.simulator.models.larter_breakspear.LarterBreakspear.QV_max = NArray(label=':math:`Q_{max}`', dtype=float64, default=array([1.]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 1.0 
   attribute  tvb.simulator.models.larter_breakspear.LarterBreakspear.QZ_max = NArray(label=':math:`Q_{max}`', dtype=float64, default=array([1.]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 1.0 
   attribute  tvb.simulator.models.larter_breakspear.LarterBreakspear.t_scale = NArray(label=':math:`t_{scale}`', dtype=float64, default=array([1.]), dim_names=(), ndim=None, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.larter_breakspear.LarterBreakspear.state_variable_range = Final(field_type=<class 'dict'>, default={'V': array([-1.5,  1.5]), 'W': array([-1.5,  1.5]), 'Z': array([-1.5,  1.5])}, required=True)
WARNING  default contains values out of the declared domain. Ex 0.27 
   attribute  tvb.simulator.models.wong_wang.ReducedWongWang.a = NArray(label=':math:`a`', dtype=float64, default=array([0.27]), dim_names=(), ndim=None, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.wong_wang.ReducedWongWang.state_variable_range = Final(field_type=<class 'dict'>, default={'S': array([0., 1.])}, required=True)
WARNING  default contains values out of the declared domain. Ex 10.0 
   attribute  tvb.simulator.models.wong_wang_exc_io_inh_i.ReducedWongWangExcIOInhI.tau_i = NArray(label=':math:`\\tau_i`', dtype=float64, default=array([10.]), dim_names=(), ndim=None, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.wong_wang_exc_io_inh_i.ReducedWongWangExcIOInhI.state_variable_range = Final(field_type=<class 'dict'>, default={'S_e': array([0., 1.]), 'S_i': array([0., 1.])}, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.linear.Linear.state_variable_range = Final(field_type=<class 'dict'>, default={'x': array([-1,  1])}, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.hopfield.Hopfield.state_variable_range = Final(field_type=<class 'dict'>, default={'x': array([-1.,  2.]), 'theta': array([0., 1.])}, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.epileptor.Epileptor.state_variable_range = Final(field_type=<class 'dict'>, default={'x1': array([-2.,  1.]), 'y1': array([-20.,   2.]), 'z': array([2., 5.]), 'x2': array([-2.,  0.]), 'y2': array([0., 2.]), 'g': array([-1.,  1.])}, required=True)
WARNING  default contains values out of the declared domain. Ex 1.0 
   attribute  tvb.simulator.models.epileptor.Epileptor2D.tt = NArray(label='tt', dtype=float64, default=array([1.]), dim_names=(), ndim=None, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.epileptor.Epileptor2D.state_variable_range = Final(field_type=<class 'dict'>, default={'x1': array([-2.,  1.]), 'z': array([2., 5.])}, required=True)
WARNING  default contains values out of the declared domain. Ex 1.0 
   attribute  tvb.simulator.models.JCepileptor.JC_Epileptor.gamma_rs = NArray(label=":math:'\\gamma_rs'", dtype=float64, default=array([1.]), dim_names=(), ndim=None, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.JCepileptor.JC_Epileptor.state_variable_range = Final(field_type=<class 'dict'>, default={'x1': array([-1.8, -1.4]), 'y1': array([-15, -10]), 'z': array([3.6, 4. ]), 'x2': array([-1.1, -0.9]), 'y2': array([0.001, 0.01 ]), 'g': array([-1.,  1.]), 'x_rs': array([-2.,  4.]), 'y_rs': array([-6.,  6.])}, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.epileptorcodim3.EpileptorCodim3.state_variable_range = Final(field_type=<class 'dict'>, default={'x': array([0.4, 0.6]), 'y': array([-0.1,  0.1]), 'z': array([0.  , 0.15])}, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.epileptorcodim3.EpileptorCodim3SlowMod.state_variable_range = Final(field_type=<class 'dict'>, default={'x': array([0.4, 0.6]), 'y': array([-0.1,  0.1]), 'z': array([0. , 0.1]), 'uA': array([0., 0.]), 'uB': array([0., 0.])}, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.Zerlaut.Zerlaut_adaptation_first_order.state_variable_range = Final(field_type=<class 'dict'>, default={'E': array([0. , 0.1]), 'I': array([0. , 0.1]), 'W': array([  0., 100.])}, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.simulator.models.Zerlaut.Zerlaut_adaptation_second_order.state_variable_range = Final(field_type=<class 'dict'>, default={'E': array([0. , 0.1]), 'I': array([0. , 0.1]), 'C_ee': array([0., 0.]), 'C_ei': array([0., 0.]), 'C_ii': array([0., 0.]), 'W': array([  0., 100.])}, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.datatypes.time_series.TimeSeries.labels_dimensions = Attr(field_type=<class 'dict'>, default={}, required=True)
WARNING  Field seems mutable and has a default value. Consider using a lambda as a value factory 
   attribute tvb.datatypes.projections.ProjectionMatrix.conductances = Attr(field_type=<class 'dict'>, default={'air': 0.0, 'skin': 1.0, 'skull': 0.01, 'brain': 1.0}, required=False)
WARNING  default contains values out of the declared domain. Ex 1.0 
   attribute  tvb.simulator.coupling.HyperbolicTangent.b = NArray(label=':math:`b`', dtype=float64, default=array([1.]), dim_names=(), ndim=None, required=True)
WARNING  default contains values out of the declared domain. Ex 1.0 
   attribute  tvb.simulator.coupling.Kuramoto.a = NArray(label=':math:`a`', dtype=float64, default=array([1.]), dim_names=(), ndim=None, required=True)

Epileptor codim 3 model

The Epileptor codim 3 model, developed by Saggio et al. 2017, is a neural mass model which has a rich set of possible bursting patterns. These bursters correspond to a seizure-like state as recorded in SEEG-recordings of epileptic patients. In this tutorial we will explore a number of classes of bursting which are present in this model.

The model consists of a fast and a slow subsystem. The fast subsystem is based on a spherical unfolding of the degenerate Takens-Bogdanov bifurcation. The slow subsystem controls the parameters of the fast subsystem along a path on the spherical unfolding. The class of a burster is determined by which bifurcations happen at onset and offset of the oscillation, which is in turn determined by the path on the spherical unfolding.

In [2]:
Epileptorcd3 = EpileptorCodim3()

First we look at the phase plane. The default parameters of the model are setup so we have a Saddle-Node bifurcation at onset and a Saddle-Homoclinic bifucation at offset. This corresponds to bifurcations in the original Epileptor model. By adjusting the variable z, we can travel on the spherical unfolding in a predetermined path. The Saddle-Node bifurcation happens at z=0 and the Saddle Homoclinic bifurcation happens at z=0.14.

In [3]:
Epileptorcd3.state_variable_range["x"] = array([-2.0, 2.0])
Epileptorcd3.state_variable_range["y"] = array([-2.0, 2.0])
Epileptorcd3.state_variable_range["z"] = array([-0.1, 0.3])
ppi_fig = PhasePlaneInteractive(model=Epileptorcd3)
ppi_fig.show()
Epileptorcd3.state_variable_range["x"] = numpy.array([0.4, 0.6])
Epileptorcd3.state_variable_range["y"] = numpy.array([-0.1, 0.1])
Epileptorcd3.state_variable_range["z"] = numpy.array([0.0, 0.1])

Next we simulate the model to validate that the desired burster happens.

In [4]:
Epileptorcd3.variables_of_interest=['x', 'y', 'z']
sim = simulator.Simulator(
    model=Epileptorcd3,
    connectivity=connectivity.Connectivity.from_file(),
    coupling=coupling.Linear(a=numpy.array([0.0152])),
    integrator=integrators.HeunDeterministic(dt=2 ** -4),
    monitors=(monitors.TemporalAverage(period=2 ** -2),),
    simulation_length=2 ** 10,
).configure()

(tavg_time, tavg_data), = sim.run()

figure()
plot(tavg_time, tavg_data[:, 0, 0, 0], label='x')
plot(tavg_time, tavg_data[:, 2, 0, 0], label='z')
legend()
grid(True)
xlabel('Time (ms)')
ylabel("Temporal average")
show()

fig2 = figure()
ax = fig2.add_subplot(111, projection='3d')
ax.plot(tavg_data[:, 0, 0, 0], tavg_data[:, 1, 0, 0], tavg_data[:, 2, 0, 0])
plt.grid(True)
plt.xlabel('x')
plt.ylabel('y')
plt.show()
WARNING  File 'hemispheres' not found in ZIP.

Classes of bursters

As stated earlier, the class of burster is affected by which path the bifurcation parameters take on the spherical unfolding of the degenerate Takens-Bogdanov bifurcation. These figures, from Saggio et al. 2017, show the bifurcation diagram of the deg. TB-bifurcation of the focus type. We parametrise theses paths using an arc of a great circle.

Unfolding_r0p4.png Unfolding_planar.PNG

Here are the parameter settings for the bursting classes which are listed above. They are constructed by choosing the vectors A and B as the starting and end point of the arc of the great circle respectively. Note that the bursters happen at different timescales.

This is a simulation of the 'c0' class which is a Saddle-Node/Saddle-Node burster. While this class behaves like a buster, it does not contain a limit cycle during the oscillatory phase.

In [5]:
A = [0.2649, -0.05246, 0.2951]
B = [0.2688, 0.05363, 0.2914]
c = 0.001

sim = simulator.Simulator(
    model= EpileptorCodim3(variables_of_interest=['x', 'y', 'z'], 
                           mu1_start=numpy.array([-A[1]]), 
                           mu2_start=numpy.array([A[0]]), 
                           nu_start=numpy.array([A[2]]),
                           mu1_stop=numpy.array([-B[1]]), 
                           mu2_stop=numpy.array([B[0]]), 
                           nu_stop=numpy.array([B[2]]), 
                           c=numpy.array([c])),
    connectivity=connectivity.Connectivity.from_file(),
    coupling=coupling.Linear(a=numpy.array([0.0152])),
    integrator=integrators.HeunDeterministic(dt=2 ** -4),
    monitors=(monitors.TemporalAverage(period=2 ** -2),),
    simulation_length=2 ** 11,
).configure()
(tavg_time, tavg_data), = sim.run()
figure()
plot(tavg_time, tavg_data[:, 0, 0, 0], label='x')
plot(tavg_time, tavg_data[:, 2, 0, 0], label='z')
legend()
grid(True)
xlabel('Time (ms)')
ylabel("Temporal average")
fig2 = figure()
ax = fig2.add_subplot(111, projection='3d')
ax.plot(tavg_data[:, 0, 0, 0], tavg_data[:, 1, 0, 0], tavg_data[:, 2, 0, 0])
plt.grid(True)
plt.xlabel('x')
plt.ylabel('y')
show()
WARNING  File 'hemispheres' not found in ZIP.

This is a simulation of the 'c2s' class which is a Saddle-Node/Saddle-Homoclinic burster. This is the default setting for the Epileptor codim 3 model and is the same burster which is the basis of the original Epileptor model.

In [6]:
A = [0.3448,0.02285,0.2014]
B = [0.3351,0.07465,0.2053]
c = 0.001 

sim = simulator.Simulator(
    model= EpileptorCodim3(variables_of_interest=['x', 'y', 'z'], 
                           mu1_start=numpy.array([-A[1]]), 
                           mu2_start=numpy.array([A[0]]), 
                           nu_start=numpy.array([A[2]]),
                           mu1_stop=numpy.array([-B[1]]), 
                           mu2_stop=numpy.array([B[0]]), 
                           nu_stop=numpy.array([B[2]]), 
                           c=numpy.array([c])),
    connectivity=connectivity.Connectivity.from_file(),
    coupling=coupling.Linear(a=numpy.array([0.0152])),
    integrator=integrators.HeunDeterministic(dt=2 ** -4),
    monitors=(monitors.TemporalAverage(period=2 ** -2),),
    simulation_length=2 ** 10,
).configure()
(tavg_time, tavg_data), = sim.run()
figure()
plot(tavg_time, tavg_data[:, 0, 0, 0], label='x')
plot(tavg_time, tavg_data[:, 2, 0, 0], label='z')
legend()
grid(True)
xlabel('Time (ms)')
ylabel("Temporal average")
fig2 = figure()
ax = fig2.add_subplot(111, projection='3d')
ax.plot(tavg_data[:, 0, 0, 0], tavg_data[:, 1, 0, 0], tavg_data[:, 2, 0, 0])
plt.grid(True)
plt.xlabel('x')
plt.ylabel('y')
show()
WARNING  File 'hemispheres' not found in ZIP.

This is a simulation of the 'c3s' class which is a Saddle-Node/Supercritical Hopf burster.

In [7]:
A = [0.2552,-0.0637,0.3014]
B = [0.3496,0.0795,0.1774]
c = 0.0004

sim = simulator.Simulator(
    model= EpileptorCodim3(variables_of_interest=['x', 'y', 'z'], 
                           mu1_start=numpy.array([-A[1]]), 
                           mu2_start=numpy.array([A[0]]), 
                           nu_start=numpy.array([A[2]]),
                           mu1_stop=numpy.array([-B[1]]), 
                           mu2_stop=numpy.array([B[0]]), 
                           nu_stop=numpy.array([B[2]]), 
                           c=numpy.array([c])),
    connectivity=connectivity.Connectivity.from_file(),
    coupling=coupling.Linear(a=numpy.array([0.0152])),
    integrator=integrators.HeunDeterministic(dt=2 ** -4),
    monitors=(monitors.TemporalAverage(period=2 ** -2),),
    simulation_length=2 ** 13,
).configure()
(tavg_time, tavg_data), = sim.run()
figure()
plot(tavg_time, tavg_data[:, 0, 0, 0], label='x')
plot(tavg_time, tavg_data[:, 2, 0, 0], label='z')
legend()
grid(True)
xlabel('Time (ms)')
ylabel("Temporal average")
fig2 = figure()
ax = fig2.add_subplot(111, projection='3d')
ax.plot(tavg_data[:, 0, 0, 0], tavg_data[:, 1, 0, 0], tavg_data[:, 2, 0, 0])
plt.grid(True)
plt.xlabel('x')
plt.ylabel('y'