Merge pull request #130 from neurolib-dev/model/wong_wang

New model: Wong Wang

examples/example-0.4-wc-minimal.ipynb

@@ -107,7 +107,7 @@
     "for exc_ext in exc_inputs:\n",
     "    # Note: this has to be a vector since it is input for all nodes\n",
     "    # (but we have only one node in this example)\n",
-    "    wc.params['exc_ext'] = [exc_ext]\n",
+    "    wc.params['exc_ext'] = exc_ext\n",
     "    wc.run()\n",
     "    # we add the maximum and the minimum of the last second of the \n",
     "    # simulation to a list\n",
@@ -234,7 +234,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "wc.params['exc_ext'] = [0.65] * wc.params['N']\n",
+    "wc.params['exc_ext'] = 0.65\n",
     "\n",
     "wc.params['signalV'] = 0\n",
     "wc.params['duration'] = 20 * 1000 \n",
@@ -297,13 +297,6 @@
     "print(\"Correlation per subject:\", [f\"{s:.2}\" for s in scores])\n",
     "print(\"Mean FC/FC correlation: {:.2f}\".format(np.mean(scores)))"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
@@ -322,7 +315,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.9"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,

neurolib/models/multimodel/builder/wilson_cowan.py

@@ -7,7 +7,7 @@
 from ..builder.base.network import Network, SingleCouplingExcitatoryInhibitoryNode
 from ..builder.base.neural_mass import NeuralMass
 
-WC_EXC_DEFAULT_PARAMS = {"a": 1.5, "mu": 3.0, "tau": 2.5, "ext_input": 1.0}
+WC_EXC_DEFAULT_PARAMS = {"a": 1.5, "mu": 3.0, "tau": 2.5, "ext_input": 0.0}
 WC_INH_DEFAULT_PARAMS = {"a": 1.5, "mu": 3.0, "tau": 3.75, "ext_input": 0.0}
 # matrix as [to, from], masses as (EXC, INH)
 WC_NODE_DEFAULT_CONNECTIVITY = np.array([[16.0, 12.0], [15.0, 3.0]])

neurolib/models/wc/loadDefaultParams.py

@@ -50,9 +50,9 @@ def loadDefaultParams(Cmat=None, Dmat=None, seed=None):
 
     # external input parameters:
     params.tau_ou = 5.0  # ms Timescale of the Ornstein-Uhlenbeck noise process
-    params.sigma_ou = 0.0  # mV/ms/sqrt(ms) noise intensity
-    params.exc_ou_mean = 0.0  # mV/ms (OU process) [0-5]
-    params.inh_ou_mean = 0.0  # mV/ms (OU process) [0-5]
+    params.sigma_ou = 0.0  # noise intensity
+    params.exc_ou_mean = 0.0  # OU process mean
+    params.inh_ou_mean = 0.0  # OU process mean
 
     # neural mass model parameters
     params.tau_exc = 2.5  # excitatory time constant
@@ -66,6 +66,10 @@ def loadDefaultParams(Cmat=None, Dmat=None, seed=None):
     params.mu_exc = 3.0  # excitatory firing threshold
     params.mu_inh = 3.0  # inhibitory firing threshold
 
+    # values of the external inputs
+    params.exc_ext = 0  # baseline external input to E
+    params.inh_ext = 0  # baseline external input to I
+
     # ------------------------------------------------------------------------
 
     params.exc_init = 0.05 * np.random.uniform(0, 1, (params.N, 1))
@@ -75,10 +79,6 @@ def loadDefaultParams(Cmat=None, Dmat=None, seed=None):
     params.exc_ou = np.zeros((params.N,))
     params.inh_ou = np.zeros((params.N,))
 
-    # values of the external inputs
-    params.exc_ext = 1.0 * np.ones((params.N,))
-    params.inh_ext = np.zeros((params.N,))
-
     return params
 
 

neurolib/models/wc/timeIntegration.py

@@ -37,9 +37,9 @@ def timeIntegration(params):
     tau_ou = params["tau_ou"]
     # Parameter of the Ornstein-Uhlenbeck (OU) process for the external input ( mV/ms/sqrt(ms) )
     sigma_ou = params["sigma_ou"]
-    # Mean external excitatory input (OU process) (mV/ms)
+    # Mean external excitatory input (OU process)
     exc_ou_mean = params["exc_ou_mean"]
-    # Mean external inhibitory input (OU process) (mV/ms)
+    # Mean external inhibitory input (OU process)
     inh_ou_mean = params["inh_ou_mean"]
 
     # ------------------------------------------------------------------------
@@ -72,6 +72,7 @@ def timeIntegration(params):
     max_global_delay = np.max(Dmat_ndt)
     startind = int(max_global_delay + 1)  # timestep to start integration at
 
+    # noise variable
     exc_ou = params["exc_ou"]
     inh_ou = params["inh_ou"]
 
@@ -95,8 +96,8 @@ def timeIntegration(params):
         inh_init = params["inh_init"][:, -startind:]
 
     # xsd = np.zeros((N,N))  # delayed activity
-    exc_input_d = np.zeros(N)  # delayed input to x
-    inh_input_d = np.zeros(N)  # delayed input to y
+    exc_input_d = np.zeros(N)  # delayed input to exc
+    inh_input_d = np.zeros(N)  # delayed input to inh (note used)
 
     np.random.seed(RNGseed)
 
@@ -202,7 +203,6 @@ def S_I(x):
 
             # delayed input to each node
             exc_input_d[no] = 0
-            inh_input_d[no] = 0
 
             for l in range(N):
                 exc_input_d[no] += K_gl * Cmat[no, l] * (excs[l, i - Dmat_ndt[no, l] - 1])
@@ -218,7 +218,7 @@ def S_I(x):
                         c_excexc * excs[no, i - 1]  # input from within the excitatory population
                         - c_inhexc * inhs[no, i - 1]  # input from the inhibitory population
                         + exc_input_d[no]  # input from other nodes
-                        + exc_ext[no]
+                        + exc_ext
                     )  # external input
                     + exc_ou[no]  # ou noise
                 )
@@ -232,8 +232,7 @@ def S_I(x):
                     * S_I(
                         c_excinh * excs[no, i - 1]  # input from the excitatory population
                         - c_inhinh * inhs[no, i - 1]  # input from within the inhibitory population
-                        + exc_input_d[no]  # input from other nodes
-                        + inh_ext[no]
+                        + inh_ext
                     )  # external input
                     + inh_ou[no]  # ou noise
                 )

neurolib/models/ww/__init__.py

@@ -0,0 +1 @@
+from .model import WWModel

neurolib/models/ww/loadDefaultParams.py

@@ -0,0 +1,109 @@
+import numpy as np
+
+from ...utils.collections import dotdict
+
+
+def loadDefaultParams(Cmat=None, Dmat=None, seed=None):
+    """Load default parameters for the Wong-Wang model
+
+    :param Cmat: Structural connectivity matrix (adjacency matrix) of coupling strengths, will be normalized to 1. If not given, then a single node simulation will be assumed, defaults to None
+    :type Cmat: numpy.ndarray, optional
+    :param Dmat: Fiber length matrix, will be used for computing the delay matrix together with the signal transmission speed parameter `signalV`, defaults to None
+    :type Dmat: numpy.ndarray, optional
+    :param seed: Seed for the random number generator, defaults to None
+    :type seed: int, optional
+
+    :return: A dictionary with the default parameters of the model
+    :rtype: dict
+    """
+
+    params = dotdict({})
+
+    ### runtime parameters
+    params.dt = 0.1  # ms 0.1ms is reasonable
+    params.duration = 2000  # Simulation duration (ms)
+    np.random.seed(seed)  # seed for RNG of noise and ICs
+    params.seed = seed
+
+    # ------------------------------------------------------------------------
+    # global whole-brain network parameters
+    # ------------------------------------------------------------------------
+
+    # signal transmission speec between areas
+    params.signalV = 20.0
+    params.K_gl = 0.6  # global coupling strength
+
+    if Cmat is None:
+        params.N = 1
+        params.Cmat = np.zeros((1, 1))
+        params.lengthMat = np.zeros((1, 1))
+
+    else:
+        params.Cmat = Cmat.copy()  # coupling matrix
+        np.fill_diagonal(params.Cmat, 0)  # no self connections
+        params.N = len(params.Cmat)  # number of nodes
+        params.lengthMat = Dmat
+
+    # ------------------------------------------------------------------------
+    # local node parameters
+    # ------------------------------------------------------------------------
+
+    # # the coupling parameter determines how nodes are coupled.
+    # # "original" for original wong-wang model, "reduced" for reduced wong-wang model
+    # params.version = "original"
+
+    # external noise parameters:
+    params.tau_ou = 5.0  # ms Timescale of the Ornstein-Uhlenbeck noise process
+    params.sigma_ou = 0.0  #  noise intensity
+    params.exc_ou_mean = 0.0  # OU process mean
+    params.inh_ou_mean = 0.0  # OU process mean
+
+    # neural mass model parameters
+    params.a_exc = 0.31  # nC^-1
+    params.b_exc = 0.125  # kHz
+    params.d_exc = 160.0  # ms
+    params.tau_exc = 100.0  # ms
+    params.gamma_exc = 0.641
+    params.w_exc = 1.0
+    params.exc_current = 0.382  # nA
+
+    params.a_inh = 0.615  # nC^-1
+    params.b_inh = 0.177  # kHz
+    params.d_inh = 87.0  # ms
+    params.tau_inh = 10.0  # ms
+    params.w_inh = 0.7
+    params.inh_current = 0.382  # nA
+
+    params.J_NMDA = 0.15  # nA, excitatory synaptic coupling
+    params.J_I = 1.0  # nA, inhibitory synaptic coupling
+    params.w_ee = 1.4  # excitatory feedback coupling strength
+
+    # ------------------------------------------------------------------------
+
+    params.ses_init = 0.05 * np.random.uniform(0, 1, (params.N, 1))
+    params.sis_init = 0.05 * np.random.uniform(0, 1, (params.N, 1))
+
+    # Ornstein-Uhlenbeck noise state variables
+    params.exc_ou = np.zeros((params.N,))
+    params.inh_ou = np.zeros((params.N,))
+
+    return params
+
+
+def computeDelayMatrix(lengthMat, signalV, segmentLength=1):
+    """Compute the delay matrix from the fiber length matrix and the signal velocity
+
+    :param lengthMat:       A matrix containing the connection length in segment
+    :param signalV:         Signal velocity in m/s
+    :param segmentLength:   Length of a single segment in mm
+
+    :returns:    A matrix of connexion delay in ms
+    """
+
+    normalizedLenMat = lengthMat * segmentLength
+    # Interareal connection delays, Dmat(i,j) in ms
+    if signalV > 0:
+        Dmat = normalizedLenMat / signalV
+    else:
+        Dmat = lengthMat * 0.0
+    return Dmat

neurolib/models/ww/model.py

@@ -0,0 +1,48 @@
+from . import loadDefaultParams as dp
+from . import timeIntegration as ti
+from ..model import Model
+
+
+class WWModel(Model):
+    """
+    Wong-Wang model. Original version and reduced version.
+
+    Main reference:
+        [original] Wong, K. F., & Wang, X. J. (2006). A recurrent network mechanism
+        of time integration in perceptual decisions. Journal of Neuroscience, 26(4),
+        1314-1328.
+
+    Additional references:
+        [reduced] Deco, G., Ponce-Alvarez, A., Mantini, D., Romani, G. L., Hagmann,
+        P., & Corbetta, M. (2013). Resting-state functional connectivity emerges
+        from structurally and dynamically shaped slow linear fluctuations. Journal
+        of Neuroscience, 33(27), 11239-11252.
+
+        [original] Deco, G., Ponce-Alvarez, A., Hagmann, P., Romani, G. L., Mantini,
+        D., & Corbetta, M. (2014). How local excitation–inhibition ratio impacts the
+        whole brain dynamics. Journal of Neuroscience, 34(23), 7886-7898.
+    """
+
+    name = "wongwang"
+    description = "Wong-Wang neural mass model"
+
+    init_vars = ["r_exc", "r_inh", "ses_init", "sis_init", "exc_ou", "inh_ou"]
+    state_vars = ["r_exc", "r_inh", "se", "si", "exc_ou", "inh_ou"]
+    output_vars = ["r_exc", "r_inh", "se", "si"]
+    default_output = "r_exc"
+
+    def __init__(self, params=None, Cmat=None, Dmat=None, seed=None):
+
+        self.Cmat = Cmat
+        self.Dmat = Dmat
+        self.seed = seed
+
+        # the integration function must be passed
+        integration = ti.timeIntegration
+
+        # load default parameters if none were given
+        if params is None:
+            params = dp.loadDefaultParams(Cmat=self.Cmat, Dmat=self.Dmat, seed=self.seed)
+
+        # Initialize base class Model
+        super().__init__(integration=integration, params=params)