Source code for tvb.simulator.models.wong_wang

# -*- 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
#
#

"""
Models based on Wong-Wang's work.

"""

import numpy
from .base import ModelNumbaDfun
from numba import guvectorize, float64
from tvb.basic.neotraits.api import NArray, Final, List, Range


@guvectorize([(float64[:],)*11], '(n),(m)' + ',()'*8 + '->(n)', nopython=True)
def _numba_dfun(S, c, a, b, d, g, ts, w, j, io, dx):
    "Gufunc for reduced Wong-Wang model equations."
    x = w[0]*j[0]*S[0] + io[0] + j[0]*c[0]
    h = (a[0]*x - b[0]) / (1 - numpy.exp(-d[0]*(a[0]*x - b[0])))
    dx[0] = - (S[0] / ts[0]) + (1.0 - S[0]) * h * g[0]


[docs] class ReducedWongWang(ModelNumbaDfun): r""" .. [WW_2006] Kong-Fatt Wong and Xiao-Jing Wang, *A Recurrent Network Mechanism of Time Integration in Perceptual Decisions*. Journal of Neuroscience 26(4), 1314-1328, 2006. .. [DPA_2013] Deco Gustavo, Ponce Alvarez Adrian, Dante Mantini, Gian Luca Romani, Patric Hagmann and Maurizio Corbetta. *Resting-State Functional Connectivity Emerges from Structurally and Dynamically Shaped Slow Linear Fluctuations*. The Journal of Neuroscience 32(27), 11239-11252, 2013. Equations taken from [DPA_2013]_ , page 11242 .. math:: x_k &= w\,J_N \, S_k + I_o + J_N \mathbf\Gamma(S_k, S_j, u_{kj})\\ H(x_k) &= \dfrac{ax_k - b}{1 - \exp(-d(ax_k -b))}\\ \dot{S}_k &= -\dfrac{S_k}{\tau_s} + (1 - S_k) \, H(x_k) \, \gamma """ # Define traited attributes for this model, these represent possible kwargs. a = NArray( label=":math:`a`", default=numpy.array([0.270, ]), domain=Range(lo=0.0, hi=0.270, step=0.01), doc="[n/C]. Input gain parameter, chosen to fit numerical solutions.") b = NArray( label=":math:`b`", default=numpy.array([0.108, ]), domain=Range(lo=0.0, hi=1.0, step=0.01), doc="[kHz]. Input shift parameter chosen to fit numerical solutions.") d = NArray( label=":math:`d`", default=numpy.array([154., ]), domain=Range(lo=0.0, hi=200.0, step=0.01), doc="""[ms]. Parameter chosen to fit numerical solutions.""") gamma = NArray( label=r":math:`\gamma`", default=numpy.array([0.641, ]), domain=Range(lo=0.0, hi=1.0, step=0.01), doc="""Kinetic parameter""") tau_s = NArray( label=r":math:`\tau_S`", default=numpy.array([100., ]), domain=Range(lo=50.0, hi=150.0, step=1.0), doc="""Kinetic parameter. NMDA decay time constant.""") w = NArray( label=r":math:`w`", default=numpy.array([0.6, ]), domain=Range(lo=0.0, hi=1.0, step=0.01), doc="""Excitatory recurrence""") J_N = NArray( label=r":math:`J_{N}`", default=numpy.array([0.2609, ]), domain=Range(lo=0.2609, hi=0.5, step=0.001), doc="""Excitatory recurrence""") I_o = NArray( label=":math:`I_{o}`", default=numpy.array([0.33, ]), domain=Range(lo=0.0, hi=1.0, step=0.01), doc="""[nA] Effective external input""") sigma_noise = NArray( label=r":math:`\sigma_{noise}`", default=numpy.array([0.000000001, ]), domain=Range(lo=0.0, hi=0.005, step=0.0001), doc="""[nA] Noise amplitude. Take this value into account for stochatic integration schemes.""") state_variable_range = Final( label="State variable ranges [lo, hi]", default={"S": numpy.array([0.0, 1.0])}, doc="Population firing rate") state_variable_boundaries = Final( label="State Variable boundaries [lo, hi]", default={"S": numpy.array([0.0, 1.0])}, doc="""The values for each state-variable should be set to encompass the boundaries of the dynamic range of that state-variable. Set None for one-sided boundaries""") variables_of_interest = List( of=str, label="Variables watched by Monitors", choices=("S",), default=("S",), doc="""default state variables to be monitored""") state_variables = ['S'] _nvar = 1 cvar = numpy.array([0], dtype=numpy.int32)
[docs] def configure(self): """ """ super(ReducedWongWang, self).configure() self.update_derived_parameters()
def _numpy_dfun(self, state_variables, coupling, local_coupling=0.0): S = state_variables[0, :] c_0 = coupling[0, :] # if applicable lc_0 = local_coupling * S x = self.w * self.J_N * S + self.I_o + self.J_N * c_0 + self.J_N * lc_0 H = (self.a * x - self.b) / (1 - numpy.exp(-self.d * (self.a * x - self.b))) dS = - (S / self.tau_s) + (1 - S) * H * self.gamma derivative = numpy.array([dS]) return derivative
[docs] def dfun(self, x, c, local_coupling=0.0): r""" Equations taken from [DPA_2013]_ , page 11242 .. math:: x_k &= w\,J_N \, S_k + I_o + J_N \mathbf\Gamma(S_k, S_j, u_{kj})\\ H(x_k) &= \dfrac{ax_k - b}{1 - \exp(-d(ax_k -b))}\\ \dot{S}_k &= -\dfrac{S_k}{\tau_s} + (1 - S_k) \, H(x_k) \, \gamma """ x_ = x.reshape(x.shape[:-1]).T c_ = c.reshape(c.shape[:-1]).T + local_coupling * x[0] deriv = _numba_dfun(x_, c_, self.a, self.b, self.d, self.gamma, self.tau_s, self.w, self.J_N, self.I_o) return deriv.T[..., numpy.newaxis]