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