(initially assembled for McIntosh Lab theory group meeting in August 2015)
Define some variables
tvb_lib_folder = '../../../../libraries_of_others/github/tvb-library'
tvb_dat_folder = '../../../../libraries_of_others/github/tvb-data'
eg_rois = ['lV1', 'rV1', 'lA1', 'rA1']
second_y_rois = ['lV1', 'rA1']
Importage
import os,sys,glob,numpy as np,pandas as pd
%matplotlib inline
from matplotlib import pyplot as plt
tvb_lib_folder = os.path.abspath(tvb_lib_folder)
tvb_dat_folder = os.path.abspath(tvb_dat_folder)
sys.path+=[tvb_lib_folder,tvb_dat_folder]
from tvb.simulator.lab import *
from mne.io import RawArray
from mne import create_info
2019-06-20 23:53:40,188 - WARNING - tvb.simulator.common - psutil module not available: no warnings will be issued when a
simulation may require more memory than available
INFO log level set to INFO
$\dot{V} = d \, \tau (-f V^3 + e V^2 + g V + \alpha W + \gamma I)$
$\dot{W} = \dfrac{d}{\tau}\,\,(c V^2 + b V - \beta W + a)$
Parameter table for the Generic 2D Oscillator model (generated automatically from traits class structure)
G = models.Generic2dOscillator
kks = ['description', 'default', 'minValue', 'maxValue']
thing = []
for k in G.trait:
intf = getattr(G,k).interface
thingy = [k]
for kk in kks:
if kk in intf:
thingy.append(intf[kk])
thing.append(thingy)
df = pd.DataFrame(thing, columns = ['parameter'] + kks)
df.set_index('parameter', inplace=True)
myparams = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'alpha',
'beta', 'gamma', 'I', 'tau']
df.ix[myparams]#.to_pickle(f)
description | default | minValue | maxValue | |
---|---|---|---|---|
parameter | ||||
a | Vertical shift of the configurable nullcline | [-2.0] | -5.0000 | 5.0 |
b | Linear slope of the configurable nullcline | [-10.0] | -20.0000 | 15.0 |
c | Parabolic term of the configurable nullcline | [0.0] | -10.0000 | 10.0 |
d | Temporal scale factor. Warning: do not use it ... | [0.02] | 0.0001 | 1.0 |
e | Coefficient of the quadratic term of the cubic... | [3.0] | -5.0000 | 5.0 |
f | Coefficient of the cubic term of the cubic nul... | [1.0] | -5.0000 | 5.0 |
g | Coefficient of the linear term of the cubic nu... | [0.0] | -5.0000 | 5.0 |
alpha | Constant parameter to scale the rate of feedba... | [1.0] | -5.0000 | 5.0 |
beta | Constant parameter to scale the rate of feedba... | [1.0] | -5.0000 | 5.0 |
gamma | Constant parameter to reproduce FHN dynamics w... | [1.0] | -1.0000 | 1.0 |
I | Baseline shift of the cubic nullcline | [0.0] | -5.0000 | 5.0 |
tau | A time-scale hierarchy can be introduced for t... | [1.0] | 1.0000 | 5.0 |
G2dO Configuration 1 - Excitable
+---------------------------+
| Table 1 |
+--------------+------------+
| EXCITABLE CONFIGURATION |
+--------------+------------+
|Parameter | Value |
+==============+============+
| a | -2.0 |
+--------------+------------+
| b | -10.0 |
+--------------+------------+
| c | 0.0 |
+--------------+------------+
| d | 0.02 |
+--------------+------------+
| I | 0.0 |
+--------------+------------+
| limit cycle if a is 2.0 |
+---------------------------+
G2dO Configuration 2 - Bistable
+---------------------------+
| Table 2 |
+--------------+------------+
| BISTABLE CONFIGURATION |
+--------------+------------+
|Parameter | Value |
+==============+============+
| a | 1.0 |
+--------------+------------+
| b | 0.0 |
+--------------+------------+
| c | -5.0 |
+--------------+------------+
| d | 0.02 |
+--------------+------------+
| I | 0.0 |
+--------------+------------+
| monostable regime: |
| fixed point if Iext=-2.0 |
| limit cycle if Iext=-1.0 |
+---------------------------+
G2dO Configuration 3 - Excitable (Morris-Lecar-like)
+---------------------------+
| Table 3 |
+--------------+------------+
| EXCITABLE CONFIGURATION |
+--------------+------------+
| (similar to Morris-Lecar)|
+--------------+------------+
|Parameter | Value |
+==============+============+
| a | 0.5 |
+--------------+------------+
| b | 0.6 |
+--------------+------------+
| c | -4.0 |
+--------------+------------+
| d | 0.02 |
+--------------+------------+
| I | 0.0 |
+--------------+------------+
| excitable regime if b=0.6 |
| oscillatory if b=0.4 |
+---------------------------+
G2dO Configuration 4 - 10 Hz Peak
+---------------------------+
| Table 4 |
+--------------+------------+
| GhoshetAl, 2008 |
| KnocketAl, 2009 |
+--------------+------------+
|Parameter | Value |
+==============+============+
| a | 1.05 |
+--------------+------------+
| b | -1.00 |
+--------------+------------+
| c | 0.0 |
+--------------+------------+
| d | 0.1 |
+--------------+------------+
| I | 0.0 |
+--------------+------------+
| alpha | 1.0 |
+--------------+------------+
| beta | 0.2 |
+--------------+------------+
| gamma | -1.0 |
+--------------+------------+
| e | 0.0 |
+--------------+------------+
| g | 1.0 |
+--------------+------------+
| f | 1/3 |
+--------------+------------+
| tau | 1.25 |
+--------------+------------+
| |
| frequency peak at 10Hz |
| |
+---------------------------+
G2dO Configuration 5 - 10 Hz Freq Peak
+---------------------------+
| Table 5 |
+--------------+------------+
| SanzLeonetAl 2013 |
+--------------+------------+
|Parameter | Value |
+==============+============+
| a | - 0.5 |
+--------------+------------+
| b | -10.0 |
+--------------+------------+
| c | 0.0 |
+--------------+------------+
| d | 0.02 |
+--------------+------------+
| I | 0.0 |
+--------------+------------+
| |
| intrinsic frequency is |
| approx 10 Hz |
| |
+---------------------------+
NOTE: This regime, if I = 2.1, is called subthreshold regime.
Unstable oscillations appear through a subcritical Hopf bifurcation.
Simple run simulation function
def run_smallsim(g2do_params,n_step=50000,dt=0.1):
mod = models.Generic2dOscillator(**g2do_params)
t,y = mod.stationary_trajectory(n_step=n_step,dt=dt)
y = np.squeeze(y)
df_VW = pd.DataFrame(y,index=t,columns=['V', 'W'])
df_VW.columns.names = ['svar']
df_VW.index.names = ['t']
info = create_info(ch_names = ['V', 'W'], ch_types = ['eeg', 'eeg'], sfreq=100./dt)
mne_VW = RawArray(info=info,data=y.T)
return df_VW,mne_VW
Specify regimes
regime_1a = {'a': -2.0, 'b': -10.0, 'c': 0.0, 'd': 0.02, 'I': 0.0}
regime_1b = {'a': 2.0, 'b': -10.0, 'c': 0.0, 'd': 0.02, 'I': 0.0}
regime_2a = {'a': 1.0, 'b': 0.0, 'c': -5.0, 'd': 0.02, 'I': 0.0}
regime_2b = {'a': 1.0, 'b': 0.0, 'c': -5.0, 'd': 0.02, 'I': -2.0}
regime_2c = {'a': 1.0, 'b': 0.0, 'c': -5.0, 'd': 0.02, 'I': -1.0}
regime_3a = {'a': 0.5, 'b': 0.6, 'c': -4.0, 'd': 0.02, 'I': 0.0}
regime_3b = {'a': 0.5, 'b': 0.4, 'c': -4.0, 'd': 0.02, 'I': 0.0}
regime_4 = {'a': 1.05, 'b': -1.00, 'c': 0.0, 'd': 0.1,
'I': 0.0, 'alpha': 1.0, 'beta': 0.2, 'gamma': -1.0,
'e': 0.0, 'g': 1.0, 'f': 1/3, 'tau': 1.25}
regime_5a = {'a': -0.5, 'b': -10.0, 'c': 0.0, 'd': 0.02, 'I': 0.0}
regime_5b = {'a': -0.5, 'b': -10.0, 'c': 0.0, 'd': 0.02, 'I': 2.1}
%%time
VW_1a,mne_1a = run_smallsim(regime_1a,n_step=10000,dt=0.5)
VW_1b,mne_1b = run_smallsim(regime_1b,n_step=10000,dt=0.5)
VW_2a,mne_2a = run_smallsim(regime_2a,n_step=10000,dt=0.5)
VW_2b,mne_2b = run_smallsim(regime_2b,n_step=10000,dt=0.5)
VW_3a,mne_3a = run_smallsim(regime_3a,n_step=10000,dt=0.5)
VW_3b,mne_3b = run_smallsim(regime_3b,n_step=10000,dt=0.5)
VW_4,mne_4 = run_smallsim(regime_4)
VW_5a,mne_5a = run_smallsim(regime_5a,n_step=10000,dt=0.5)
VW_5b,mne_5b = run_smallsim(regime_5b,n_step=10000,dt=0.5)
Creating RawArray with float64 data, n_channels=2, n_times=1001
Range : 0 ... 1000 = 0.000 ... 5.000 secs
Ready.
Creating RawArray with float64 data, n_channels=2, n_times=1001
Range : 0 ... 1000 = 0.000 ... 5.000 secs
Ready.
Creating RawArray with float64 data, n_channels=2, n_times=1001
Range : 0 ... 1000 = 0.000 ... 5.000 secs
Ready.
Creating RawArray with float64 data, n_channels=2, n_times=1001
Range : 0 ... 1000 = 0.000 ... 5.000 secs
Ready.
Creating RawArray with float64 data, n_channels=2, n_times=1001
Range : 0 ... 1000 = 0.000 ... 5.000 secs
Ready.
Creating RawArray with float64 data, n_channels=2, n_times=1001
Range : 0 ... 1000 = 0.000 ... 5.000 secs
Ready.
Creating RawArray with float64 data, n_channels=2, n_times=5001
Range : 0 ... 5000 = 0.000 ... 5.000 secs
Ready.
Creating RawArray with float64 data, n_channels=2, n_times=1001
Range : 0 ... 1000 = 0.000 ... 5.000 secs
Ready.
Creating RawArray with float64 data, n_channels=2, n_times=1001
Range : 0 ... 1000 = 0.000 ... 5.000 secs
Ready.
CPU times: user 8.45 s, sys: 203 ms, total: 8.66 s
Wall time: 8.95 s
fig, ax = plt.subplots(ncols=2,nrows=2, figsize=(12,6))
VW_1a.loc[1000:3000].plot(ax=ax[0][0])
VW_1a.plot(x='V', y='W',ax=ax[1][0])#,xlim=[-5,5], ylim=[-5,5])
VW_1b.loc[1000:3000].plot(ax=ax[0][1])
VW_1b.plot(x='V', y='W',ax=ax[1][1])#,xlim=[-5,5], ylim=[-5,5])
fig, ax = plt.subplots(ncols=2,nrows=1, figsize=(12,3))
mne_1a.plot_psd(ax=ax[0],show=False);
mne_1b.plot_psd(ax=ax[1],show=False);
Effective window size : 5.005 (s)
Effective window size : 5.005 (s)
fig, ax = plt.subplots(ncols=2,nrows=2, figsize=(12,6))
VW_2a.loc[1000:3000].plot(ax=ax[0][0]);
VW_2a.plot(x='V', y='W',ax=ax[1][0]);
VW_2b.loc[1000:3000].plot(ax=ax[0][1]);
VW_2b.plot(x='V', y='W',ax=ax[1][1]);
fig, ax = plt.subplots(ncols=2,nrows=1, figsize=(12,3))
_ = mne_2a.plot_psd(ax=ax[0],show=False,fmin=0.,fmax=60.);
_ = mne_2b.plot_psd(ax=ax[1],show=False,fmin=0.,fmax=60.);
Effective window size : 5.005 (s)
Effective window size : 5.005 (s)
fig, ax = plt.subplots(ncols=2,nrows=2, figsize=(12,6))
VW_3a.loc[1000:3000].plot(ax=ax[0][0]);
VW_3a.plot(x='V', y='W',ax=ax[1][0]);
VW_3b.loc[1000:3000].plot(ax=ax[0][1]);
VW_3b.plot(x='V', y='W',ax=ax[1][1]);
fig, ax = plt.subplots(ncols=2,nrows=1, figsize=(12,3));
_ = mne_3a.plot_psd(ax=ax[0],show=False,fmin=0.,fmax=60.);
_ = mne_3b.plot_psd(ax=ax[1],show=False,fmin=0.,fmax=60.);
Effective window size : 5.005 (s)
Effective window size : 5.005 (s)
fig, ax = plt.subplots(ncols=2,nrows=2, figsize=(12,6))
VW_5a.loc[1000:3000].plot(ax=ax[0][0])
VW_5a.plot(x='V', y='W',ax=ax[1][0])#,xlim=[-5,5], ylim=[-5,5])
VW_5b.loc[1000:3000].plot(ax=ax[0][1])
VW_5b.plot(x='V', y='W',ax=ax[1][1])#,xlim=[-5,5], ylim=[-5,5])
fig, ax = plt.subplots(ncols=2,nrows=1, figsize=(12,3))
_ = mne_5a.plot_psd(ax=ax[0],show=False,fmin=0.,fmax=50.);
_ = mne_5b.plot_psd(ax=ax[1],show=False,fmin=0.,fmax=50.);
Effective window size : 5.005 (s)
Effective window size : 5.005 (s)
fig, ax = plt.subplots(ncols=2,nrows=2, figsize=(12,6))
VW_4.loc[1000:3000].plot(ax=ax[0][0]);
VW_4.plot(x='V', y='W',ax=ax[1][0]);
ax[0][1].axis('off');
ax[1][1].axis('off');
fig, ax = plt.subplots(ncols=2,nrows=1, figsize=(12,3));
_ = mne_4.plot_psd(ax=ax[0],show=False,fmin=0.,fmax=60.);
ax[1].axis('off');
#_ = mne_2b.plot_psd(ax=ax[1],show=False,fmin=0.,fmax=60.);
Effective window size : 2.048 (s)
(0.0, 1.0, 0.0, 1.0)