N-Point Surface Simulations in TVB

(by analogy with Jirsa ‘2-point’ connections)

Notebook Setup

Define some variables

outdir = '/tmp/n-point_tvb_surface_sims'
!mkdir -p $outdir

import os
tvb_folder = os.path.abspath('../../../../../Code/libraries_of_others/github')
ctx_file = tvb_folder + '/tvb-data/tvb_data/surfaceData/cortex_16384.zip'

Importage

import os,sys,glob,numpy as np,pandas as pd
from copy import deepcopy

%matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns

sys.path += [tvb_folder + '/' + t for t in ['tvb-library', 'tvb-data']]
from tvb.simulator.lab import *
from tvb.datatypes.region_mapping import RegionMapping
from tvb.datatypes.cortex import Cortex

from nilearn.plotting import plot_surf
2019-10-21 19:44:27,866 - 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
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.epileptor_rs.EpileptorRestingState.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.epileptor_rs.EpileptorRestingState.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.hopfield.Hopfield.state_variable_range = Final(field_type=<class 'dict'>, default={'x': array([-1.,  2.]), 'theta': array([0., 1.])}, 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 -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  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  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  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  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  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_inh.ReducedWongWangExcInh.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_inh.ReducedWongWangExcInh.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.zerlaut.ZerlautFirstOrder.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.ZerlautSecondOrder.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)


/nethome/kcni/jgriffiths/Software/miniconda3/envs/jupyter_py3/lib/python3.7/site-packages/sklearn/externals/joblib/__init__.py:15: DeprecationWarning: sklearn.externals.joblib is deprecated in 0.21 and will be removed in 0.23. Please import this functionality directly from joblib, which can be installed with: pip install joblib. If this warning is raised when loading pickled models, you may need to re-serialize those models with scikit-learn 0.21+.
  warnings.warn(msg, category=DeprecationWarning)

Go to output dir

os.chdir(outdir)

Ok, let’s get cracking…

origctx = cortex.Cortex.from_file(ctx_file)
origctx.configure()

origconn = connectivity.Connectivity.from_file('connectivity_76.zip')
origconn.configure()
WARNING  File 'hemispheres' not found in ZIP.
lA1_nodenum = 38
rA1_nodenum = 0
rA1lA1_nodenums = np.array([rA1_nodenum,lA1_nodenum])

#newmat = np.zeros_like(conn.weights)
#newmat[rA1_nodenum,lA1_nodenum] = 3
#newmat[lA1_nodenum,rA1_nodenum] = 3

newmat = np.array([[0,1],[1,0]])
newconn = connectivity.Connectivity()
newconn.weights = newmat
newconn.region_labels = origconn.region_labels[rA1lA1_nodenums]
newconn.centres = origconn.centres[rA1lA1_nodenums]
newconn.cortical = origconn.cortical[rA1lA1_nodenums]
newconn.hemispheres = origconn.hemispheres[rA1lA1_nodenums]
newconn.tract_lengths = origconn.tract_lengths[rA1lA1_nodenums,:][:,rA1lA1_nodenums]
newconn.configure()

lA1_rmidxs_orig = np.nonzero(origctx.region_mapping==lA1_nodenum)[0]
rA1_rmidxs_orig = np.nonzero(origctx.region_mapping==rA1_nodenum)[0]
alltris_in_lA1_idxs_orig = []
alltris_in_rA1_idxs_orig = []

for l_it,l in enumerate(origctx.triangles):
  if (l[0] in lA1_rmidxs_orig and l[1] in lA1_rmidxs_orig and l[2] in lA1_rmidxs_orig): 
    alltris_in_lA1_idxs_orig.append(l_it)
    
  if (l[0] in rA1_rmidxs_orig and l[1] in rA1_rmidxs_orig and l[2] in rA1_rmidxs_orig): 
    alltris_in_rA1_idxs_orig.append(l_it)    
    
    
# Tests: these should error if any of the triangle values are not in the 'alltris' arrays 

assert len([l for l in origctx.triangles[alltris_in_lA1_idxs_orig].ravel() if l not in lA1_rmidxs_orig]) == 0
assert len([l for l in origctx.triangles[alltris_in_rA1_idxs_orig].ravel() if l not in rA1_rmidxs_orig]) == 0


newtris = origctx.triangles[alltris_in_rA1_idxs_orig+alltris_in_lA1_idxs_orig]
newrm = np.array([rA1_nodenum for r in rA1_rmidxs_orig] + [lA1_nodenum for l in lA1_rmidxs_orig])
newrm_v2 = newrm.copy()
newrm_v2[newrm==38] = 1

newvtxids_un = np.unique(newtris)
newvtcs = origctx.vertices[newvtxids_un]
newvtcs_key = {k: k_it for k_it,k in enumerate(newvtxids_un)}
newtris_renamed = np.array([[newvtcs_key[kk] for kk in k] for k in newtris])


new_vertex_triangles = [origctx.vertex_triangles[r] for r in newvtxids_un]
# test: newtris_renamed should pick out same array from newvtcs as newtris picks out from orig vtcs
assert len(np.nonzero((newvtcs[newtris_renamed].ravel() == origctx.vertices[newtris].ravel())==False)[0]) == 0
newctx = Cortex()
#newctx = cortex.CorticalSurface()
newctx.coupling_strength = np.array([2**-10])

newctx.vertices = newvtcs
newctx.triangles = newtris_renamed
#newctx.region_mapping_data = newrm
#setattr(newctx,'region_mapping', newrm)

#tmprm = RegionMappingData()
tmprm = RegionMapping()
tmprm.array_data = newrm_v2
newctx.region_mapping_data = tmprm
#newctx.vertex_triangles = new_vertex_triangles
#newctx.vertex_normals = origctx.vertex_normals[newvtxids_un]
newctx.configure()

model = models.Generic2dOscillator()
solver = integrators.HeunStochastic()
#cpl = coupling.Linear(a=0.1)
cpl = coupling.Linear(a=np.array([0.1]))
mons = (monitors.TemporalAverage(period=2**-2),)


#Initialise Simulator -- Model, Connectivity, Integrator, Monitors, and surface.
sim = simulator.Simulator(model = model, connectivity = newconn,
                          coupling = cpl,
                          integrator = solver,monitors = mons,
                          surface = newctx)
sim.configure()

Simulator

value
Type
Simulator
conduction_speed
3.0
connectivity
Connectivity
coupling
Linear
initial_conditions
None
integrator
HeunStochastic
model
Generic2dOscillator
monitors
(,)
simulation_length
1000.0
stimulus
None
surface
Cortex
title
Simulator
res = sim.run(simulation_length=500)

tavg_time,tavg_dat = res[0]
df_tavg = pd.DataFrame(np.squeeze(tavg_dat),index=tavg_time)


df_tavg.plot(legend=False,figsize=(12,3),linewidth=0.01)


ax = plot_surf([newctx.vertices,newctx.triangles]);


fig, ax =plt.subplots(ncols=3,figsize=(12,3))

x,y,z = newctx.vertices.T

ax[0].scatter(x,y,c='b',s=50)
ax[1].scatter(x,z,c='b',s=50)
ax[2].scatter(y,z,c='b',s=50)


x,y,z = origctx.vertices.T
ax[0].scatter(x,y,linewidth=0.0001,alpha=0.05,c='k',s=5)
ax[1].scatter(x,z,linewidth=0.001,alpha=0.05,c='k',s=5)
ax[2].scatter(y,z,linewidth=0.0001,alpha=0.05,c='k',s=5)
<matplotlib.collections.PathCollection at 0x7fe94dbac748>

png

png

png