Source code for tvb.simulator.models.oscillator

# -*- coding: utf-8 -*-
#
#
# TheVirtualBrain-Scientific Package. This package holds all simulators, and
# analysers necessary to run brain-simulations. You can use it stand alone or
# in conjunction with TheVirtualBrain-Framework Package. See content of the
# documentation-folder for more details. See also http://www.thevirtualbrain.org
#
# (c) 2012-2023, Baycrest Centre for Geriatric Care ("Baycrest") and others
#
# This program is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software Foundation,
# either version 3 of the License, or (at your option) any later version.
# This program is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
# PARTICULAR PURPOSE.  See the GNU General Public License for more details.
# You should have received a copy of the GNU General Public License along with this
# program.  If not, see <http://www.gnu.org/licenses/>.
#
#
#   CITATION:
# When using The Virtual Brain for scientific publications, please cite it as explained here:
# https://www.thevirtualbrain.org/tvb/zwei/neuroscience-publications
#
#

"""
Oscillator models.

"""

from .base import Model, ModelNumbaDfun
import numpy
from numba import guvectorize, float64
from tvb.basic.neotraits.api import NArray, Final, List, Range
from tvb.simulator.backend.ref import RefBase


[docs] class Generic2dOscillator(ModelNumbaDfun): r""" The Generic2dOscillator model is a generic dynamic system with two state variables. The dynamic equations of this model are composed of two ordinary differential equations comprising two nullclines. The first nullcline is a cubic function as it is found in most neuron and population models; the second nullcline is arbitrarily configurable as a polynomial function up to second order. The manipulation of the latter nullcline's parameters allows to generate a wide range of different behaviours. Equations: .. math:: \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) See: .. [FH_1961] FitzHugh, R., *Impulses and physiological states in theoretical models of nerve membrane*, Biophysical Journal 1: 445, 1961. .. [Nagumo_1962] Nagumo et.al, *An Active Pulse Transmission Line Simulating Nerve Axon*, Proceedings of the IRE 50: 2061, 1962. .. [SJ_2011] Stefanescu, R., Jirsa, V.K. *Reduced representations of heterogeneous mixed neural networks with synaptic coupling*. Physical Review E, 83, 2011. .. [SJ_2010] Jirsa VK, Stefanescu R. *Neural population modes capture biologically realistic large-scale network dynamics*. Bulletin of Mathematical Biology, 2010. .. [SJ_2008_a] Stefanescu, R., Jirsa, V.K. *A low dimensional description of globally coupled heterogeneous neural networks of excitatory and inhibitory neurons*. PLoS Computational Biology, 4(11), 2008). The model's (:math:`V`, :math:`W`) time series and phase-plane its nullclines can be seen in the figure below. The model with its default parameters exhibits FitzHugh-Nagumo like dynamics. +---------------------------+ | 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 | +---------------------------+ +---------------------------+ | 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 | +---------------------------+ +---------------------------+ | 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 | +---------------------------+ +---------------------------+ | 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 | | | +---------------------------+ +---------------------------+ | 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. .. figure :: img/Generic2dOscillator_01_mode_0_pplane.svg .. _phase-plane-Generic2D: :alt: Phase plane of the generic 2D population model with (V, W) The (:math:`V`, :math:`W`) phase-plane for the generic 2D population model for default parameters. The dynamical system has an equilibrium point. .. automethod:: Generic2dOscillator.dfun """ # Define traited attributes for this model, these represent possible kwargs. tau = NArray( label=r":math:`\tau`", default=numpy.array([1.0]), domain=Range(lo=1.0, hi=5.0, step=0.01), doc="""A time-scale hierarchy can be introduced for the state variables :math:`V` and :math:`W`. Default parameter is 1, which means no time-scale hierarchy.""") I = NArray( label=":math:`I_{ext}`", default=numpy.array([0.0]), domain=Range(lo=-5.0, hi=5.0, step=0.01), doc="""Baseline shift of the cubic nullcline""") a = NArray( label=":math:`a`", default=numpy.array([-2.0]), domain=Range(lo=-5.0, hi=5.0, step=0.01), doc="""Vertical shift of the configurable nullcline""") b = NArray( label=":math:`b`", default=numpy.array([-10.0]), domain=Range(lo=-20.0, hi=15.0, step=0.01), doc="""Linear slope of the configurable nullcline""") c = NArray( label=":math:`c`", default=numpy.array([0.0]), domain=Range(lo=-10.0, hi=10.0, step=0.01), doc="""Parabolic term of the configurable nullcline""") d = NArray( label=":math:`d`", default=numpy.array([0.02]), domain=Range(lo=0.0001, hi=1.0, step=0.0001), doc="""Temporal scale factor. Warning: do not use it unless you know what you are doing and know about time tides.""") e = NArray( label=":math:`e`", default=numpy.array([3.0]), domain=Range(lo=-5.0, hi=5.0, step=0.0001), doc="""Coefficient of the quadratic term of the cubic nullcline.""") f = NArray( label=":math:`f`", default=numpy.array([1.0]), domain=Range(lo=-5.0, hi=5.0, step=0.0001), doc="""Coefficient of the cubic term of the cubic nullcline.""") g = NArray( label=":math:`g`", default=numpy.array([0.0]), domain=Range(lo=-5.0, hi=5.0, step=0.5), doc="""Coefficient of the linear term of the cubic nullcline.""") alpha = NArray( label=r":math:`\alpha`", default=numpy.array([1.0]), domain=Range(lo=-5.0, hi=5.0, step=0.0001), doc="""Constant parameter to scale the rate of feedback from the slow variable to the fast variable.""") beta = NArray( label=r":math:`\beta`", default=numpy.array([1.0]), domain=Range(lo=-5.0, hi=5.0, step=0.0001), doc="""Constant parameter to scale the rate of feedback from the slow variable to itself""") # This parameter is basically a hack to avoid having a negative lower boundary in the global coupling strength. gamma = NArray( label=r":math:`\gamma`", default=numpy.array([1.0]), domain=Range(lo=-1.0, hi=1.0, step=0.1), doc="""Constant parameter to reproduce FHN dynamics where excitatory input currents are negative. It scales both I and the long range coupling term.""") state_variable_range = Final( label="State Variable ranges [lo, hi]", default={"V": numpy.array([-2.0, 4.0]), "W": numpy.array([-6.0, 6.0])}, doc="""The values for each state-variable should be set to encompass the expected dynamic range of that state-variable for the current parameters, it is used as a mechanism for bounding random initial conditions when the simulation isn't started from an explicit history, it is also provides the default range of phase-plane plots.""") variables_of_interest = List( of=str, label="Variables or quantities available to Monitors", choices=("V", "W", "V + W", "V - W"), default=("V",), doc="The quantities of interest for monitoring for the generic 2D oscillator.") state_variables = ('V', 'W') _nvar = 2 cvar = numpy.array([0], dtype=numpy.int32) def _numpy_dfun(self, state_variables, coupling, local_coupling=0.0): V = state_variables[0, :] W = state_variables[1, :] # [State_variables, nodes] c_0 = coupling[0, :] tau = self.tau I = self.I a = self.a b = self.b c = self.c d = self.d e = self.e f = self.f g = self.g beta = self.beta alpha = self.alpha gamma = self.gamma lc_0 = local_coupling * V # Pre-allocate the result array then instruct numexpr to use it as output. # This avoids an expensive array concatenation derivative = numpy.empty_like(state_variables) ev = RefBase.evaluate ev('d * tau * (alpha * W - f * V**3 + e * V**2 + g * V + gamma * I + gamma *c_0 + lc_0)', out=derivative[0]) ev('d * (a + b * V + c * V**2 - beta * W) / tau', out=derivative[1]) return derivative
[docs] def dfun(self, vw, c, local_coupling=0.0): r""" The two state variables :math:`V` and :math:`W` are typically considered to represent a function of the neuron's membrane potential, such as the firing rate or dendritic currents, and a recovery variable, respectively. If there is a time scale hierarchy, then typically :math:`V` is faster than :math:`W` corresponding to a value of :math:`\tau` greater than 1. The equations of the generic 2D population model read .. math:: \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) where external currents :math:`I` provide the entry point for local, long-range connectivity and stimulation. """ lc_0 = local_coupling * vw[0, :, 0] vw_ = vw.reshape(vw.shape[:-1]).T c_ = c.reshape(c.shape[:-1]).T deriv = _numba_dfun_g2d(vw_, c_, self.tau, self.I, self.a, self.b, self.c, self.d, self.e, self.f, self.g, self.beta, self.alpha, self.gamma, lc_0) return deriv.T[..., numpy.newaxis]
@guvectorize([(float64[:],) * 16], '(n),(m)' + ',()' * 13 + '->(n)', nopython=True) def _numba_dfun_g2d(vw, c_0, tau, I, a, b, c, d, e, f, g, beta, alpha, gamma, lc_0, dx): "Gufunc for Generic2dOscillator model equations." V = vw[0] V2 = V * V W = vw[1] dx[0] = d[0] * tau[0] * ( alpha[0] * W - f[0] * V2 * V + e[0] * V2 + g[0] * V + gamma[0] * I[0] + gamma[0] * c_0[0] + lc_0[0]) dx[1] = d[0] * (a[0] + b[0] * V + c[0] * V2 - beta[0] * W) / tau[0]
[docs] class Kuramoto(Model): r""" The Kuramoto model is a model of synchronization phenomena derived by Yoshiki Kuramoto in 1975 which has since been applied to diverse domains including the study of neuronal oscillations and synchronization. See: .. [YK_1975] Y. Kuramoto, in: H. Arakai (Ed.), International Symposium on Mathematical Problems in Theoretical Physics, *Lecture Notes in Physics*, page 420, vol. 39, 1975. .. [SS_2000] S. H. Strogatz. *From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators*. Physica D, 143, 2000. .. [JC_2011] J. Cabral, E. Hugues, O. Sporns, G. Deco. *Role of local network oscillations in resting-state functional connectivity*. NeuroImage, 57, 1, 2011. The :math:`\theta` variable is the phase angle of the oscillation. Dynamic equations: .. math:: \dot{\theta}_{k} = \omega_{k} + \mathbf{\Gamma}(\theta_k, \theta_j, u_{kj}) + \sin(W_{\zeta}\theta) """ # Define traited attributes for this model, these represent possible kwargs. omega = NArray( label=r":math:`\omega`", default=numpy.array([1.0]), domain=Range(lo=0.01, hi=200.0, step=0.1), doc=r""":math:`\omega` sets the base line frequency for the Kuramoto oscillator in [rad/ms]""") state_variable_range = Final( label="State Variable ranges [lo, hi]", default={"theta": numpy.array([0.0, numpy.pi * 2.0]), }, doc="""The values for each state-variable should be set to encompass the expected dynamic range of that state-variable for the current parameters, it is used as a mechanism for bounding random initial conditions when the simulation isn't started from an explicit history, it is also provides the default range of phase-plane plots.""") variables_of_interest = List( of=str, label="Variables watched by Monitors", choices=("theta",), default=("theta",), doc="""This represents the default state-variables of this Model to be monitored. It can be overridden for each Monitor if desired. The Kuramoto model, however, only has one state variable with and index of 0, so it is not necessary to change the default here.""") state_variables = ['theta'] _nvar = 1 cvar = numpy.array([0], dtype=numpy.int32)
[docs] def dfun(self, state_variables, coupling, local_coupling=0.0, ev=RefBase.evaluate, sin=numpy.sin, pi2=numpy.pi * 2): r""" The :math:`\theta` variable is the phase angle of the oscillation. .. math:: \dot{\theta}_{k} = \omega_{k} + \mathbf{\Gamma}(\theta_k, \theta_j, u_{kj}) + \sin(W_{\zeta}\theta) where :math:`I` is the input via local and long range connectivity, passing first through the Kuramoto coupling function, """ theta = state_variables[0, :] # import pdb; pdb.set_trace() # A) Distribution of phases according to the local connectivity kernel local_range_coupling = numpy.sin(local_coupling * theta) # NOTE: To evaluate. # B) Strength of the interactions # local_range_coupling = local_coupling * numpy.sin(theta) I = coupling[0, :] + local_range_coupling if not hasattr(self, 'derivative'): self.derivative = numpy.empty((1,) + theta.shape) # phase update self.derivative[0] = self.omega + I # all this pi makeh me have great hungary, can has sum NaN? return self.derivative
[docs] class SupHopf(ModelNumbaDfun): r""" The supHopf model describes the normal form of a supercritical Hopf bifurcation in Cartesian coordinates. This normal form has a supercritical bifurcation at a=0 with a the bifurcation parameter in the model. So for a < 0, the local dynamics has a stable fixed point and the system corresponds to a damped oscillatory state, whereas for a > 0, the local dynamics enters in a stable limit cycle and the system switches to an oscillatory state. See for examples: .. [Kuznetsov_2013] Kuznetsov, Y.A. *Elements of applied bifurcation theory.* Springer Sci & Business Media, 2013, vol. 112. .. [Deco_2017a] Deco, G., Kringelbach, M.L., Jirsa, V.K., Ritter, P. *The dynamics of resting fluctuations in the brain: metastability and its dynamical cortical core* Sci Reports, 2017, 7: 3095. The equations of the supHopf equations read as follows: .. math:: \dot{x}_{i} &= (a_{i} - x_{i}^{2} - y_{i}^{2})x_{i} - {\omega}{i}y_{i} \\ \dot{y}_{i} &= (a_{i} - x_{i}^{2} - y_{i}^{2})y_{i} + {\omega}{i}x_{i} where a is the local bifurcation parameter and omega the angular frequency. """ a = NArray( label=r":math:`a`", default=numpy.array([-0.5]), domain=Range(lo=-10.0, hi=10.0, step=0.01), doc="""Local bifurcation parameter.""") omega = NArray( label=r":math:`\omega`", default=numpy.array([1.]), domain=Range(lo=0.05, hi=630.0, step=0.01), doc="""Angular frequency.""") # Initialization. state_variable_range = Final( label="State Variable ranges [lo, hi]", default={"x": numpy.array([-5.0, 5.0]), "y": numpy.array([-5.0, 5.0])}, doc="""The values for each state-variable should be set to encompass the expected dynamic range of that state-variable for the current parameters, it is used as a mechanism for bounding random initial conditions when the simulation isn't started from an explicit history, it is also provides the default range of phase-plane plots.""") variables_of_interest = List( of=str, label="Variables watched by Monitors", choices=("x", "y"), default=("x",), doc="Quantities of supHopf available to monitor.") state_variables = ["x", "y"] _nvar = 2 # number of state-variables cvar = numpy.array([0, 1], dtype=numpy.int32) # coupling variables def _numpy_dfun(self, state_variables, coupling, local_coupling=0.0, array=numpy.array, where=numpy.where, concat=numpy.concatenate): y = state_variables ydot = numpy.empty_like(state_variables) # long-range coupling c_0 = coupling[0] c_1 = coupling[1] # short-range (local) coupling lc_0 = local_coupling * y[0] # supHopf's equations in Cartesian coordinates: ydot[0] = (self.a - y[0] ** 2 - y[1] ** 2) * y[0] - self.omega * y[1] + c_0 + lc_0 ydot[1] = (self.a - y[0] ** 2 - y[1] ** 2) * y[1] + self.omega * y[0] + c_1 return ydot
[docs] def dfun(self, x, c, local_coupling=0.0): r""" Computes the derivatives of the state-variables of supHopf with respect to time. The equations of the supHopf equations read as follows: .. math:: \dot{x}_{i} &= (a_{i} - x_{i}^{2} - y_{i}^{2})x_{i} - {\omega}{i}y_{i} \\ \dot{y}_{i} &= (a_{i} - x_{i}^{2} - y_{i}^{2})y_{i} + {\omega}{i}x_{i} where a is the local bifurcation parameter and omega the angular frequency. """ x_ = x.reshape(x.shape[:-1]).T c_ = c.reshape(c.shape[:-1]).T lc_0 = local_coupling * x[0, :, 0] deriv = _numba_dfun_supHopf(x_, c_, self.a, self.omega, lc_0) return deriv.T[..., numpy.newaxis]
@guvectorize([(float64[:],) * 6], '(n),(m)' + ',()' * 3 + '->(n)', nopython=True) def _numba_dfun_supHopf(y, c, a, omega, lc_0, ydot): "Gufunc for supHopf model equations." # long-range coupling c_0 = c[0] c_1 = c[1] # supHopf equations ydot[0] = (a[0] - y[0] ** 2 - y[1] ** 2) * y[0] - omega[0] * y[1] + c_0 + lc_0[0] ydot[1] = (a[0] - y[0] ** 2 - y[1] ** 2) * y[1] + omega[0] * y[0] + c_1