Defining new models in TVB using the TVB DSL

In this demo we show how to use the TVB DSL in order to automatically generate model code and execute it in the TVB simulator. The TVB DSL is an extension of LEMS and is defined using XML syntax.

First let us load the necessary functions to generate code

In [1]:
import tvb.dsl.LEMS2python as templating
from IPython.display import Code

In this demo we will define a new model to TVB called the Montbrio model. This model is based on: Montbri├│, Ernest, Diego Paz├│, and Alex Roxin. "Macroscopic description for networks of spiking neurons." Physical Review X 5.2 (2015): 021028. This is a model with two state variables derived from on a network of all-to-all coupled heterogeneous quadratic integrate-and-fire QIF neurons, which is exact in the thermodynamic limit.

Getting to know the TVB DSL

In the example of the montbrio model, we can see different elements which allow us to express the model variables and dynamics. We can now take a deeper look into each component.

Lets start with the definition of the model name and description documentation. This is done by defining a new Component using the keyword ComponentType.

<ComponentType name="MontbrioT" description="This is a demo describing the Montbrio model" value="">

Constants for a model are defined using the Constant keyword, such as in:

<Constant name="I" domain="lo=-10.0, hi=10.0, step=0.01" default="0.0" description="A constant"/>

The domain specification allows TVB to check if a value set by the user to this constant before simulation is within the expected range or not.

The model dynamics are specified in a special section tagged as Dynamics. Inside this section we can define state variables using the keyword StateVariable, other variables derived from combinations of constants and state variables using the keyword DerivedVariable, and the time derivatives of the model with the keyword TimeDerivative. Some examples are illustrated below:

<StateVariable name="r" default="0., 2.0" boundaries="0.0, inf"/>

The default term defines the range within which the default value of this state variable can be defined before the simulation takes place. The boundaries specify hard limits for the values of these variables in order to avoid numerical overflows.

<DerivedVariable name="Coupling_global" expression="alpha * coupling[0]"/>

The expression term in the DerivedVariable is used to state the combination of constanst / variables that define the new variable. The experssion is expected to be in Python-like syntax

<TimeDerivative name="dx" expression="Delta / pi + 2 * V * r - k * r**2 + Gamma * r / pi"/>

As in the derived variable, the expression term in the TimeDerivative is used to mathematically describe the update in the state variable that will be performed on each integration step and should be expressed in Python-like syntax.

It is important to highlight that the order of the definition of the state variables must be preserved for the definition of their time derivatives in order for the DSL to correctly match them.

Finally, outside of the dynamics section we can use the keyword Exposure to define the state variables we are interested on monitoring from the outside e.g.:

<Exposure name="r" choices="r, V" default="r, V" description="The quantities of interest for monitoring for the Infinite QIF 2D oscillator."/>

We can now take a look at the whole XML code which is used to define the Montbrio model. The model in this example has been defined using the DSL in the file /tvb/dsl/NeuroML/XMLmodels.

In [2]:
import tvb.dsl
import os

xmlfile = open(os.path.join(os.path.dirname(tvb.dsl.__file__), "NeuroML/XMLmodels/montbriot.xml"),"r")
model =
display(Code(model, language='xml'))
      xsi:schemaLocation=" ../../LEMS/Schemas/LEMS/LEMS_v0.7.4.xsd"
      description="Rate based/population models translated using LEMS.">

    <ComponentType name="MontbrioT"
                   description="2D model describing the Ott-Antonsen reduction of infinitely all-to-all coupled QIF neurons (Theta-neurons)."

        <!-- If empty then none -->
        <Constant name="I" domain="lo=-10.0, hi=10.0, step=0.01" default="0.0" description="???"/>
        <Constant name="Delta" domain="lo=0.0, hi=10.0, step=0.01" default="1.0" description="Vertical shift of the configurable nullcline."/>
        <Constant name="alpha" domain="lo=0.0, hi=1.0, step=0.1" default="1.0" description=":math:`\alpha` ratio of effect between long-range and local connectivity."/>
        <Constant name="s" domain="lo=-15.0, hi=15.0, step=0.01" default="0.0" description="QIF membrane reversal potential."/>
        <Constant name="k" domain="lo=-15.0, hi=15.0, step=0.01" default="0.0" description="Switch for the terms specific to Coombes model."/>
        <Constant name="J" domain="lo=-25.0, hi=25.0, step=0.0001" default="15.0" description="Constant parameter to scale the rate of feedback from the slow variable to the firing rate variable."/>
        <Constant name="eta" domain="lo=-10.0, hi=10.0, step=0.0001" default="-5.0" description="Constant parameter to scale the rate of feedback from the firing rate variable to itself"/>
        <Constant name="Gamma" domain="lo=0., hi=10.0, step=0.1" default="0.0" description="Derived from eterogeneous currents and synaptic weights (see Montbrio p.12)."/>
        <Constant name="gamma" domain="lo=-2.0, hi=2.0, step=0.1" default="1.0" description="Constant parameter to reproduce FHN dynamics where excitatory input currents are negative. It scales both I and the long range coupling term."/>

            <!-- "State variable ranges [lo, hi]" values are entered with keyword "default" -->
            <!-- For each state variable a set of boundaries can be added to encompass the boundaries of the dynamic range -->
            <!-- Leave empty "" for no boundaries. Set None for one-sided boundaries, i.e. "1.0, None" -->
            <StateVariable name="r" default="0., 2.0" boundaries="0.0, inf"/>
            <StateVariable name="V" default="-2.0, 1.5" boundaries=""/>

            <!-- Derived variables can be used to simplify the time derivatives, enter the local coupling formulas or any formula -->
            <!-- syntax: [name]=[expression] -->
            <!-- Define for ex. global and local coupling: c_0 = coupling[0, ] and lc_0 = local_coupling -->
            <DerivedVariable name="Coupling_global" expression="alpha * coupling[0]"/>
            <DerivedVariable name="Coupling_local" expression="(1-alpha) * local_coupling * r"/>
            <DerivedVariable name="Coupling_Term" expression="Coupling_global + Coupling_local"/>

            <!-- For conditionals use &lt(=); or &gt;(=) for less- or greater then (equal to)  -->
            <!-- Conditional used for if statement, syntax: if {condition} -> {cases[0]} else {cases[1]}. Cases are separated by (,) -->
            <TimeDerivative name="dx" expression="Delta / pi + 2 * V * r - k * r**2 + Gamma * r / pi"/>
            <TimeDerivative name="dy" expression="V**2 - pi**2 * r**2 + eta + (k * s + J) * r - k * V * r + gamma * I + Coupling_Term"/>


        <!-- Exposures are used for observables, for the name enter variable to be observed (usually states)
        and for dimension enter the reduction functionality. Will be represented as variables_of_interest.
        Choices and default list with a (,) separator-->
        <Exposure name="r" choices="r, V" default="r, V" description="The quantities of interest for monitoring for the Infinite QIF 2D oscillator."/>


Automatically generating the model code and integrating it into TVB

We will call the templating function in order to automatically generate the model code. This will be directly moved to the right path in the tvb library so it will be recognized by the simulator when we call it later. Please notice that the name provided to the templating function needs to be the exact name given to the model ComponentType in xml.

In [3]:
2020-08-10 20:02:34,664 - INFO - tvb.dsl.LEMS2python - model file generated MontbrioT

Testing the new model in a TVB simulation

In [13]:
from tvb.simulator.lab import *
import numpy as np
from tvb.simulator import models

sim_length = 400
white_matter = connectivity.Connectivity.from_file(source_file="")
white_matter_coupling = coupling.Linear(a=np.array([1.0]))
SC = white_matter.weights
integrator = integrators.EulerDeterministic(dt=0.1)
populations = models.MontbrioT()
monitorsen = (monitors.TemporalAverage(period=10.0))

sim = simulator.Simulator(model=populations, 

(_,tavg_data) =[0]
WARNING  File 'cortical' not found in ZIP.
WARNING  File 'hemispheres' not found in ZIP.
WARNING  File 'areas' not found in ZIP.
[[[[ 0.23047686]
   [ 0.18730459]
   [ 0.10174776]
   [ 0.11056869]
   [ 0.13605601]
   [ 0.21458262]]

  [[-1.5721524 ]
   [-1.3126255 ]

 [[[ 0.08186933]
   [ 0.08153367]
   [ 0.08115131]
   [ 0.08207059]
   [ 0.08121951]
   [ 0.08235958]]


 [[[ 0.08143515]
   [ 0.08126934]
   [ 0.0811502 ]
   [ 0.08170015]
   [ 0.08119331]
   [ 0.08169049]]



 [[[ 0.08128816]
   [ 0.08120127]
   [ 0.08115008]
   [ 0.08128428]
   [ 0.08118454]
   [ 0.08139337]]


 [[[ 0.08128816]
   [ 0.08120127]
   [ 0.08115008]
   [ 0.08128428]
   [ 0.08118454]
   [ 0.08139337]]


 [[[ 0.08128816]
   [ 0.08120127]
   [ 0.08115008]
   [ 0.08128428]
   [ 0.08118454]
   [ 0.08139337]]

In [10]:
import tvb.simulator.models

In [11]:
from tvb.simulator.models.montbriot import MontbrioT
In [ ]: