# -*- 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
#
#
"""
Larter-Breakspear model based on the Morris-Lecar equations.
"""
import numpy
from .base import Model
from tvb.basic.neotraits.api import NArray, Final, List, Range
[docs]class LarterBreakspear(Model):
r"""
A modified Morris-Lecar model that includes a third equation which simulates
the effect of a population of inhibitory interneurons synapsing on
the pyramidal cells.
.. [Larteretal_1999] Larter et.al. *A coupled ordinary differential equation
lattice model for the simulation of epileptic seizures.* Chaos. 9(3):
795, 1999.
.. [Breaksetal_2003_a] Breakspear, M.; Terry, J. R. & Friston, K. J. *Modulation of excitatory
synaptic coupling facilitates synchronization and complex dynamics in an
onlinear model of neuronal dynamics*. Neurocomputing 52–54 (2003).151–158
.. [Breaksetal_2003_b] M. J. Breakspear et.al. *Modulation of excitatory
synaptic coupling facilitates synchronization and complex dynamics in a
biophysical model of neuronal dynamics.* Network: Computation in Neural
Systems 14: 703-732, 2003.
.. [Honeyetal_2007] Honey, C.; Kötter, R.; Breakspear, M. & Sporns, O. * Network structure of
cerebral cortex shapes functional connectivity on multiple time scales*. (2007)
PNAS, 104, 10240
.. [Honeyetal_2009] Honey, C. J.; Sporns, O.; Cammoun, L.; Gigandet, X.; Thiran, J. P.; Meuli,
R. & Hagmann, P. *Predicting human resting-state functional connectivity
from structural connectivity.* (2009), PNAS, 106, 2035-2040
.. [Alstottetal_2009] Alstott, J.; Breakspear, M.; Hagmann, P.; Cammoun, L. & Sporns, O.
*Modeling the impact of lesions in the human brain*. (2009)), PLoS Comput Biol, 5, e1000408
Equations and default parameters are taken from [Breaksetal_2003_b]_.
All equations and parameters are non-dimensional and normalized.
For values of d_v < 0.55, the dynamics of a single column settles onto a
solitary fixed point attractor.
Parameters used for simulations in [Breaksetal_2003_a]_ Table 1. Page 153.
Two nodes were coupled. C=0.1
+---------------------------+
| Table 1 |
+--------------+------------+
|Parameter | Value |
+==============+============+
| I | 0.3 |
+--------------+------------+
| a_ee | 0.4 |
+--------------+------------+
| a_ei | 0.1 |
+--------------+------------+
| a_ie | 1.0 |
+--------------+------------+
| a_ne | 1.0 |
+--------------+------------+
| a_ni | 0.4 |
+--------------+------------+
| r_NMDA | 0.2 |
+--------------+------------+
| delta | 0.001 |
+--------------+------------+
| Breakspear et al. 2003 |
+---------------------------+
+---------------------------+
| Table 2 |
+--------------+------------+
|Parameter | Value |
+==============+============+
| gK | 2.0 |
+--------------+------------+
| gL | 0.5 |
+--------------+------------+
| gNa | 6.7 |
+--------------+------------+
| gCa | 1.0 |
+--------------+------------+
| a_ne | 1.0 |
+--------------+------------+
| a_ni | 0.4 |
+--------------+------------+
| a_ee | 0.36 |
+--------------+------------+
| a_ei | 2.0 |
+--------------+------------+
| a_ie | 2.0 |
+--------------+------------+
| VK | -0.7 |
+--------------+------------+
| VL | -0.5 |
+--------------+------------+
| VNa | 0.53 |
+--------------+------------+
| VCa | 1.0 |
+--------------+------------+
| phi | 0.7 |
+--------------+------------+
| b | 0.1 |
+--------------+------------+
| I | 0.3 |
+--------------+------------+
| r_NMDA | 0.25 |
+--------------+------------+
| C | 0.1 |
+--------------+------------+
| TCa | -0.01 |
+--------------+------------+
| d_Ca | 0.15 |
+--------------+------------+
| TK | 0.0 |
+--------------+------------+
| d_K | 0.3 |
+--------------+------------+
| VT | 0.0 |
+--------------+------------+
| ZT | 0.0 |
+--------------+------------+
| TNa | 0.3 |
+--------------+------------+
| d_Na | 0.15 |
+--------------+------------+
| d_V | 0.65 |
+--------------+------------+
| d_Z | d_V |
+--------------+------------+
| QV_max | 1.0 |
+--------------+------------+
| QZ_max | 1.0 |
+--------------+------------+
| Alstott et al. 2009 |
+---------------------------+
NOTES about parameters
:math:`\delta_V` : for :math:`\delta_V` < 0.55, in an uncoupled network,
the system exhibits fixed point dynamics; for 0.55 < :math:`\delta_V` < 0.59,
limit cycle attractors; and for :math:`\delta_V` > 0.59 chaotic attractors
(eg, d_V=0.6,aee=0.5,aie=0.5, gNa=0, Iext=0.165)
:math:`\delta_Z`
this parameter might be spatialized: ones(N,1).*0.65 + modn*(rand(N,1)-0.5);
:math:`C`
The long-range coupling :math:`\delta_C` is ‘weak’ in the sense that
the model is well behaved for parameter values for which C < a_ee and C << a_ie.
.. figure :: img/LarterBreakspear_01_mode_0_pplane.svg
:alt: Larter-Breaskpear phase plane (V, W)
The (:math:`V`, :math:`W`) phase-plane for the Larter-Breakspear model.
Dynamic equations:
.. math::
\dot{V}_k & = - (g_{Ca} + (1 - C) \, r_{NMDA} \, a_{ee} \, Q_V + C \, r_{NMDA} \, a_{ee} \, \langle Q_V\rangle^{k}) \, m_{Ca} \, (V - VCa) \\
& \,\,- g_K \, W \, (V - VK) - g_L \, (V - VL) \\
& \,\,- (g_{Na} \, m_{Na} + (1 - C) \, a_{ee} \, Q_V + C \, a_{ee} \, \langle Q_V\rangle^{k}) \,(V - VNa) \\
& \,\,- a_{ie} \, Z \, Q_Z + a_{ne} \, I \\
& \\
\dot{W}_k & = \phi \, \dfrac{m_K - W}{\tau_{K}} \\
& \nonumber\\
\dot{Z}_k &= b (a_{ni}\, I + a_{ei}\,V\,Q_V) \\
Q_{V} &= Q_{V_{max}} \, (1 + \tanh\left(\dfrac{V_{k} - VT}{\delta_{V}}\right)) \\
Q_{Z} &= Q_{Z_{max}} \, (1 + \tanh\left(\dfrac{Z_{k} - ZT}{\delta_{Z}}\right))
See Equations (7), (3), (6) and (2) respectively in [Breaksetal_2003_a]_.
Pag: 705-706
"""
# Define traited attributes for this model, these represent possible kwargs.
gCa = NArray(
label=":math:`g_{Ca}`",
default=numpy.array([1.1]),
domain=Range(lo=0.9, hi=1.5, step=0.1),
doc="""Conductance of population of Ca++ channels.""")
gK = NArray(
label=":math:`g_{K}`",
default=numpy.array([2.0]),
domain=Range(lo=1.95, hi= 2.05, step=0.025),
doc="""Conductance of population of K channels.""")
gL = NArray(
label=":math:`g_{L}`",
default=numpy.array([0.5]),
domain=Range(lo=0.45 , hi=0.55, step=0.05),
doc="""Conductance of population of leak channels.""")
phi = NArray(
label=r":math:`\phi`",
default=numpy.array([0.7]),
domain=Range(lo=0.3, hi=0.9, step=0.1),
doc="""Temperature scaling factor.""")
gNa = NArray(
label=":math:`g_{Na}`",
default=numpy.array([6.7]),
domain=Range(lo=0.0, hi=10.0, step=0.1),
doc="""Conductance of population of Na channels.""")
TK = NArray(
label=":math:`T_{K}`",
default=numpy.array([0.0]),
domain=Range(lo=0.0, hi=0.0001, step=0.00001),
doc="""Threshold value for K channels.""")
TCa = NArray(
label=":math:`T_{Ca}`",
default=numpy.array([-0.01]),
domain=Range(lo=-0.02, hi=-0.01, step=0.0025),
doc="Threshold value for Ca channels.")
TNa = NArray(
label=":math:`T_{Na}`",
default=numpy.array([0.3]),
domain=Range(lo=0.25, hi= 0.3, step=0.025),
doc="Threshold value for Na channels.")
VCa = NArray(
label=":math:`V_{Ca}`",
default=numpy.array([1.0]),
domain=Range(lo=0.9, hi=1.1, step=0.05),
doc="""Ca Nernst potential.""")
VK = NArray(
label=":math:`V_{K}`",
default=numpy.array([-0.7]),
domain=Range(lo=-0.8, hi=1., step=0.1),
doc="""K Nernst potential.""")
VL = NArray(
label=":math:`V_{L}`",
default=numpy.array([-0.5]),
domain=Range(lo=-0.7, hi=-0.4, step=0.1),
doc="""Nernst potential leak channels.""")
VNa = NArray(
label=":math:`V_{Na}`",
default=numpy.array([0.53]),
domain=Range(lo=0.51, hi=0.55, step=0.01),
doc="""Na Nernst potential.""")
d_K = NArray(
label=r":math:`\delta_{K}`",
default=numpy.array([0.3]),
domain=Range(lo=0.1, hi=0.4, step=0.1),
doc="""Variance of K channel threshold.""")
tau_K = NArray(
label=r":math:`\tau_{K}`",
default=numpy.array([1.0]),
domain=Range(lo=1.0, hi=10.0, step=1.0),
doc="""Time constant for K relaxation time (ms)""")
d_Na = NArray(
label=r":math:`\delta_{Na}`",
default=numpy.array([0.15]),
domain=Range(lo=0.1, hi=0.2, step=0.05),
doc="Variance of Na channel threshold.")
d_Ca = NArray(
label=r":math:`\delta_{Ca}`",
default=numpy.array([0.15]),
domain=Range(lo=0.1, hi=0.2, step=0.05),
doc="Variance of Ca channel threshold.")
aei = NArray(
label=":math:`a_{ei}`",
default=numpy.array([2.0]),
domain=Range(lo=0.1, hi=2.0, step=0.1),
doc="""Excitatory-to-inhibitory synaptic strength.""")
aie = NArray(
label=":math:`a_{ie}`",
default=numpy.array([2.0]),
domain=Range(lo=0.5, hi=2.0, step=0.1),
doc="""Inhibitory-to-excitatory synaptic strength.""")
b = NArray(
label=":math:`b`",
default=numpy.array([0.1]),
domain=Range(lo=0.0001, hi=1.0, step=0.0001),
doc="""Time constant scaling factor. The original value is 0.1""")
C = NArray(
label=":math:`C`",
default=numpy.array([0.1]),
domain=Range(lo=0.0, hi=1.0, step=0.01),
doc="""Strength of excitatory coupling. Balance between internal and
local (and global) coupling strength. C > 0 introduces interdependences between
consecutive columns/nodes. C=1 corresponds to maximum coupling between node and no self-coupling.
This strenght should be set to sensible values when a whole network is connected. """)
ane = NArray(
label=":math:`a_{ne}`",
default=numpy.array([1.0]),
domain=Range(lo=0.4, hi=1.0, step=0.05),
doc="""Non-specific-to-excitatory synaptic strength.""")
ani = NArray(
label=":math:`a_{ni}`",
default=numpy.array([0.4]),
domain=Range(lo=0.3, hi=0.5, step=0.05),
doc="""Non-specific-to-inhibitory synaptic strength.""")
aee = NArray(
label=":math:`a_{ee}`",
default=numpy.array([0.4]),
domain=Range(lo=0.0, hi=0.6, step=0.05),
doc="""Excitatory-to-excitatory synaptic strength.""")
Iext = NArray(
label=":math:`I_{ext}`",
default=numpy.array([0.3]),
domain=Range(lo=0.165, hi=0.3, step=0.005),
doc="""Subcortical input strength. It represents a non-specific
excitation or thalamic inputs.""")
rNMDA = NArray(
label=":math:`r_{NMDA}`",
default=numpy.array([0.25]),
domain=Range(lo=0.2, hi=0.3, step=0.05),
doc="""Ratio of NMDA to AMPA receptors.""")
VT = NArray(
label=":math:`V_{T}`",
default=numpy.array([0.0]),
domain=Range(lo=0.0, hi=0.7, step=0.01),
doc="""Threshold potential (mean) for excitatory neurons.
In [Breaksetal_2003_b]_ this value is 0.""")
d_V = NArray(
label=r":math:`\delta_{V}`",
default=numpy.array([0.65]),
domain=Range(lo=0.49, hi=0.7, step=0.01),
doc="""Variance of the excitatory threshold. It is one of the main
parameters explored in [Breaksetal_2003_b]_.""")
ZT = NArray(
label=":math:`Z_{T}`",
default=numpy.array([0.0]),
domain=Range(lo=0.0, hi=0.1, step=0.005),
doc="""Threshold potential (mean) for inihibtory neurons.""")
d_Z = NArray(
label=r":math:`\delta_{Z}`",
default=numpy.array([0.7]),
domain=Range(lo=0.001, hi=0.75, step=0.05),
doc="""Variance of the inhibitory threshold.""")
# NOTE: the values were not in the article.
QV_max = NArray(
label=":math:`QV_{max}`",
default=numpy.array([1.0]),
domain=Range(lo=0.1, hi=1., step=0.001),
doc="""Maximal firing rate for excitatory populations (kHz)""")
QZ_max = NArray(
label=":math:`QZ_{max}`",
default=numpy.array([1.0]),
domain=Range(lo=0.1, hi=1., step=0.001),
doc="""Maximal firing rate for excitatory populations (kHz)""")
t_scale = NArray(
label=":math:`t_{scale}`",
default=numpy.array([1.0]),
domain=Range(lo=0.1, hi=1., step=0.001),
doc="""Time scale factor""")
variables_of_interest = List(
of=str,
label="Variables watched by Monitors",
choices=("V", "W", "Z"),
default=("V",),
doc="""This represents the default state-variables of this Model to be
monitored. It can be overridden for each Monitor if desired.""")
# Informational attribute, used for phase-plane and initial()
state_variable_range = Final(
label="State Variable ranges [lo, hi]",
default={
"V": numpy.array([-1.5, 1.5]),
"W": numpy.array([-1.5, 1.5]),
"Z": numpy.array([-1.5, 1.5])},
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 inital
conditions when the simulation isn't started from an explicit
history, it is also provides the default range of phase-plane plots.""")
state_variables = tuple('V W Z'.split())
_state_variables = ("V", "W", "Z")
_nvar = 3
cvar = numpy.array([0], dtype=numpy.int32)
[docs] def dfun(self, state_variables, coupling, local_coupling=0.0):
r"""
Dynamic equations:
.. math::
\dot{V}_k & = - (g_{Ca} + (1 - C) \, r_{NMDA} \, a_{ee} \, Q_V + C \, r_{NMDA} \, a_{ee} \, \langle Q_V\rangle^{k}) \, m_{Ca} \, (V - VCa) \\
& \,\,- g_K \, W \, (V - VK) - g_L \, (V - VL) \\
& \,\,- (g_{Na} \, m_{Na} + (1 - C) \, a_{ee} \, Q_V + C \, a_{ee} \, \langle Q_V\rangle^{k}) \,(V - VNa) \\
& \,\,- a_{ie} \, Z \, Q_Z + a_{ne} \, I \\
& \\
\dot{W}_k & = \phi \, \dfrac{m_K - W}{\tau_{K}} \\
& \nonumber\\
\dot{Z}_k &= b (a_{ni}\, I + a_{ei}\,V\,Q_V) \\
Q_{V} &= Q_{V_{max}} \, (1 + \tanh\left(\dfrac{V_{k} - VT}{\delta_{V}}\right)) \\
Q_{Z} &= Q_{Z_{max}} \, (1 + \tanh\left(\dfrac{Z_{k} - ZT}{\delta_{Z}}\right))
"""
V, W, Z = state_variables
derivative = numpy.empty_like(state_variables)
c_0 = coupling[0, :]
# relationship between membrane voltage and channel conductance
m_Ca = 0.5 * (1 + numpy.tanh((V - self.TCa) / self.d_Ca))
m_Na = 0.5 * (1 + numpy.tanh((V - self.TNa) / self.d_Na))
m_K = 0.5 * (1 + numpy.tanh((V - self.TK ) / self.d_K))
# voltage to firing rate
QV = 0.5 * self.QV_max * (1 + numpy.tanh((V - self.VT) / self.d_V))
QZ = 0.5 * self.QZ_max * (1 + numpy.tanh((Z - self.ZT) / self.d_Z))
lc_0 = local_coupling * QV
derivative[0] = self.t_scale * (- (self.gCa + (1.0 - self.C) * (self.rNMDA * self.aee) * (QV + lc_0)+ self.C * self.rNMDA * self.aee * c_0) * m_Ca * (V - self.VCa)
- self.gK * W * (V - self.VK)
- self.gL * (V - self.VL)
- (self.gNa * m_Na + (1.0 - self.C) * self.aee * (QV + lc_0) + self.C * self.aee * c_0) * (V - self.VNa)
- self.aie * Z * QZ
+ self.ane * self.Iext)
derivative[1] = self.t_scale * self.phi * (m_K - W) / self.tau_K
derivative[2] = self.t_scale * self.b * (self.ani * self.Iext + self.aei * V * QV)
return derivative
