Repository: neurolib-dev/neurolib

Issue:
Bug: default input values for WC model
The default parameters of the WC model have this weird difference in E vs I. I'm not sure why this is in there, maybe @ChristophMetzner can tell, but I think this needs to change to be zero for both.

```
    # values of the external inputs
    params.exc_ext = 1.0 * np.ones((params.N,))
    params.inh_ext = np.zeros((params.N,))
```

Discussion:

Patch:
diff --git a/examples/example-0.4-wc-minimal.ipynb b/examples/example-0.4-wc-minimal.ipynb
index 76b1ddb..62eb74e 100644
--- a/examples/example-0.4-wc-minimal.ipynb
+++ b/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,
diff --git a/neurolib/models/multimodel/builder/wilson_cowan.py b/neurolib/models/multimodel/builder/wilson_cowan.py
index 1de22ba..9c7e898 100644
--- a/neurolib/models/multimodel/builder/wilson_cowan.py
+++ b/neurolib/models/multimodel/builder/wilson_cowan.py
@@ -7,7 +7,7 @@ from ..builder.base.constants import EXC, INH
 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]])
diff --git a/neurolib/models/wc/loadDefaultParams.py b/neurolib/models/wc/loadDefaultParams.py
index 34726ec..788dfe0 100644
--- a/neurolib/models/wc/loadDefaultParams.py
+++ b/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
 
 
diff --git a/neurolib/models/wc/timeIntegration.py b/neurolib/models/wc/timeIntegration.py
index 67af8c2..67dfee9 100644
--- a/neurolib/models/wc/timeIntegration.py
+++ b/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 timeIntegration_njit_elementwise(
 
             # 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 timeIntegration_njit_elementwise(
                         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 timeIntegration_njit_elementwise(
                     * 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
                 )
diff --git a/neurolib/models/ww/__init__.py b/neurolib/models/ww/__init__.py
new file mode 100644
index 0000000..576f30f
--- /dev/null
+++ b/neurolib/models/ww/__init__.py
@@ -0,0 +1,1 @@
+from .model import WWModel
diff --git a/neurolib/models/ww/loadDefaultParams.py b/neurolib/models/ww/loadDefaultParams.py
new file mode 100644
index 0000000..05d6810
--- /dev/null
+++ b/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
diff --git a/neurolib/models/ww/model.py b/neurolib/models/ww/model.py
new file mode 100644
index 0000000..12d72eb
--- /dev/null
+++ b/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)
diff --git a/neurolib/models/ww/timeIntegration.py b/neurolib/models/ww/timeIntegration.py
new file mode 100644
index 0000000..fcaeabf
--- /dev/null
+++ b/neurolib/models/ww/timeIntegration.py
@@ -0,0 +1,282 @@
+import numpy as np
+import numba
+
+from . import loadDefaultParams as dp
+
+
+def timeIntegration(params):
+    """Sets up the parameters for time integration
+
+    :param params: Parameter dictionary of the model
+    :type params: dict
+    :return: Integrated activity variables of the model
+    :rtype: (numpy.ndarray,)
+    """
+
+    dt = params["dt"]  # Time step for the Euler intergration (ms)
+    duration = params["duration"]  # imulation duration (ms)
+    RNGseed = params["seed"]  # seed for RNG
+
+    # ------------------------------------------------------------------------
+    # local parameters
+    a_exc = params["a_exc"]
+    b_exc = params["b_exc"]
+    d_exc = params["d_exc"]
+    tau_exc = params["tau_exc"]
+    gamma_exc = params["gamma_exc"]
+    w_exc = params["w_exc"]
+    exc_current = params["exc_current"]
+
+    a_inh = params["a_inh"]
+    b_inh = params["b_inh"]
+    d_inh = params["d_inh"]
+    tau_inh = params["tau_exc"]
+    w_inh = params["w_inh"]
+    inh_current = params["inh_current"]
+
+    J_NMDA = params["J_NMDA"]
+    J_I = params["J_I"]
+    w_ee = params["w_ee"]
+
+    # external input parameters:
+    # Parameter of the Ornstein-Uhlenbeck process for the external input(ms)
+    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)
+    exc_ou_mean = params["exc_ou_mean"]
+    # Mean external inhibitory input (OU process)
+    inh_ou_mean = params["inh_ou_mean"]
+    # ------------------------------------------------------------------------
+    # global coupling parameters
+
+    # Connectivity matrix
+    # Interareal relative coupling strengths (values between 0 and 1), Cmat(i,j) connnection from jth to ith
+    Cmat = params["Cmat"]
+    N = len(Cmat)  # Number of nodes
+    K_gl = params["K_gl"]  # global coupling strength
+    # Interareal connection delay
+    lengthMat = params["lengthMat"]
+    signalV = params["signalV"]
+
+    if N == 1:
+        Dmat = np.zeros((N, N))
+    else:
+        # Interareal connection delays, Dmat(i,j) Connnection from jth node to ith (ms)
+        Dmat = dp.computeDelayMatrix(lengthMat, signalV)
+        # no self-feedback delay
+        Dmat[np.eye(len(Dmat)) == 1] = np.zeros(len(Dmat))
+    Dmat_ndt = np.around(Dmat / dt).astype(int)  # delay matrix in multiples of dt
+    params["Dmat_ndt"] = Dmat_ndt
+
+    # # Additive or diffusive coupling scheme
+    # version = params["version"]
+    # # convert to integer for faster integration later
+    # if version == "original":
+    #     version = 0
+    # elif version == "reduced":
+    #     version = 1
+    # else:
+    #     raise ValueError('Paramter "version" must be either "original" or "reduced"')
+
+    # ------------------------------------------------------------------------
+
+    # Initialization
+    # Floating point issue in np.arange() workaraound: use integers in np.arange()
+    t = np.arange(1, round(duration, 6) / dt + 1) * dt  # Time variable (ms)
+
+    sqrt_dt = np.sqrt(dt)
+
+    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"]
+
+    # state variable arrays, have length of t + startind
+    # they store initial conditions AND simulated data
+    ses = np.zeros((N, startind + len(t)))
+    sis = np.zeros((N, startind + len(t)))
+
+    # holds firing rates
+    r_exc = np.zeros((N, startind + len(t)))
+    r_inh = np.zeros((N, startind + len(t)))
+
+    # ------------------------------------------------------------------------
+    # Set initial values
+    # if initial values are just a Nx1 array
+    if np.shape(params["ses_init"])[1] == 1:
+        ses_init = np.dot(params["ses_init"], np.ones((1, startind)))
+        sis_init = np.dot(params["sis_init"], np.ones((1, startind)))
+    # if initial values are a Nxt array
+    else:
+        ses_init = params["ses_init"][:, -startind:]
+        sis_init = params["sis_init"][:, -startind:]
+
+    # xsd = np.zeros((N,N))  # delayed activity
+    ses_input_d = np.zeros(N)  # delayed input to x
+
+    np.random.seed(RNGseed)
+
+    # Save the noise in the activity array to save memory
+    ses[:, startind:] = np.random.standard_normal((N, len(t)))
+    sis[:, startind:] = np.random.standard_normal((N, len(t)))
+
+    ses[:, :startind] = ses_init
+    sis[:, :startind] = sis_init
+
+    noise_se = np.zeros((N,))
+    noise_si = np.zeros((N,))
+
+    # ------------------------------------------------------------------------
+
+    return timeIntegration_njit_elementwise(
+        startind,
+        t,
+        dt,
+        sqrt_dt,
+        duration,
+        N,
+        Cmat,
+        Dmat,
+        K_gl,
+        signalV,
+        Dmat_ndt,
+        ses,
+        sis,
+        ses_input_d,
+        a_exc,
+        b_exc,
+        d_exc,
+        tau_exc,
+        gamma_exc,
+        w_exc,
+        exc_current,
+        a_inh,
+        b_inh,
+        d_inh,
+        tau_inh,
+        w_inh,
+        inh_current,
+        J_NMDA,
+        J_I,
+        w_ee,
+        r_exc,
+        r_inh,
+        noise_se,
+        noise_si,
+        exc_ou,
+        inh_ou,
+        exc_ou_mean,
+        inh_ou_mean,
+        tau_ou,
+        sigma_ou,
+    )
+
+
+@numba.njit
+def timeIntegration_njit_elementwise(
+    startind,
+    t,
+    dt,
+    sqrt_dt,
+    duration,
+    N,
+    Cmat,
+    Dmat,
+    K_gl,
+    signalV,
+    Dmat_ndt,
+    ses,
+    sis,
+    ses_input_d,
+    a_exc,
+    b_exc,
+    d_exc,
+    tau_exc,
+    gamma_exc,
+    w_exc,
+    exc_current,
+    a_inh,
+    b_inh,
+    d_inh,
+    tau_inh,
+    w_inh,
+    inh_current,
+    J_NMDA,
+    J_I,
+    w_ee,
+    r_exc,
+    r_inh,
+    noise_se,
+    noise_si,
+    exc_ou,
+    inh_ou,
+    exc_ou_mean,
+    inh_ou_mean,
+    tau_ou,
+    sigma_ou,
+):
+    """
+    Wong-Wang model equations (Deco2014, no long-rage feedforward inhibition):
+
+    currents
+    I_e = w_e * I_0 + w_ee * J_NMDA * s_e - J_I * s_i + K * J_NMDA * \sum Gij * s_e_j + I_ext
+    I_i = w_i * I_0 + J_NMDA * s_e - s_i
+
+    synaptic activity
+    d_se/dt = (-s / tau) + (1.0 - s) * gamma * r + noise
+    d_si/dt = (-s / tau) + r + noise
+
+    firing rate transfer function
+    r = (a * I - b) / (1.0 - exp(-d * (a * I - b)))
+
+    """
+    #  firing rate transfer function
+    def r(I, a, b, d):
+        return (a * I - b) / (1.0 - np.exp(-d * (a * I - b)))
+
+    ### integrate ODE system:
+    for i in range(startind, startind + len(t)):
+
+        # loop through all the nodes
+        for no in range(N):
+
+            # To save memory, noise is saved in the activity array
+            noise_se[no] = ses[no, i]
+            noise_si[no] = sis[no, i]
+
+            # delayed input to each node
+            ses_input_d[no] = 0
+
+            # input from other nodes
+            for l in range(N):
+                ses_input_d[no] += K_gl * Cmat[no, l] * (ses[l, i - Dmat_ndt[no, l] - 1])
+
+            # Wong-Wang
+            se = ses[no, i - 1]
+            si = sis[no, i - 1]
+
+            I_exc = w_exc * exc_current + w_ee * J_NMDA * se - J_I * si + J_NMDA * ses_input_d[no]
+            I_inh = w_inh * inh_current + J_NMDA * se - si
+
+            r_exc[no, i] = r(I_exc, a_exc, b_exc, d_exc)
+            r_inh[no, i] = r(I_inh, a_inh, b_inh, d_inh)
+
+            se_rhs = -(se / tau_exc) + (1 - se) * gamma_exc * r_exc[no, i] + exc_ou[no]  # exc_ou = ou noise
+            si_rhs = -(si / tau_inh) + r_inh[no, i] + inh_ou[no]
+
+            # Euler integration
+            ses[no, i] = ses[no, i - 1] + dt * se_rhs
+            sis[no, i] = sis[no, i - 1] + dt * si_rhs
+
+            # Ornstein-Uhlenberg process
+            exc_ou[no] = (
+                exc_ou[no] + (exc_ou_mean - exc_ou[no]) * dt / tau_ou + sigma_ou * sqrt_dt * noise_se[no]
+            )  # mV/ms
+            inh_ou[no] = (
+                inh_ou[no] + (inh_ou_mean - inh_ou[no]) * dt / tau_ou + sigma_ou * sqrt_dt * noise_si[no]
+            )  # mV/ms
+
+    return t, r_exc, r_inh, ses, sis, exc_ou, inh_ou
diff --git a/setup.py b/setup.py
index bdd7a67..1bd4e08 100644
--- a/setup.py
+++ b/setup.py
@@ -11,7 +11,7 @@ with open("requirements.txt") as f:
 
 setuptools.setup(
     name="neurolib",
-    version="0.5.11",
+    version="0.5.12",
     description="Easy whole-brain neural mass modeling",
     long_description=long_description,
     long_description_content_type="text/markdown",


Test patch:
diff --git a/tests/test_autochunk.py b/tests/test_autochunk.py
index 99a6909..87a2478 100644
--- a/tests/test_autochunk.py
+++ b/tests/test_autochunk.py
@@ -1,5 +1,6 @@
 import copy
 import unittest
+import logging
 
 import numpy as np
 import pytest
@@ -8,6 +9,7 @@ from neurolib.models.fhn import FHNModel
 from neurolib.models.hopf import HopfModel
 from neurolib.models.thalamus import ThalamicMassModel
 from neurolib.models.wc import WCModel
+from neurolib.models.ww import WWModel
 from neurolib.utils.loadData import Dataset
 
 
@@ -23,21 +25,22 @@ class AutochunkTests(unittest.TestCase):
     def single_node_test(self, model):
         for duration in self.durations:
             for chunksize in self.chunksizes:
-                for signalV in self.signalVs:
-                    # full run
-                    m1 = model()
-                    m1.params.signalV = signalV
-                    m1.params["duration"] = duration
-                    pars_bak = copy.deepcopy(m1.params)
-                    m1.run()
-                    # chunkwise run
-                    m2 = model()
-                    m2.params = pars_bak.copy()
-                    m2.run(chunkwise=True, chunksize=chunksize, append=True)
-                    # check
-                    self.assertTupleEqual(m1.output.shape, m2.output.shape)
-                    difference = np.sum(np.abs(m1.output - m2.output))
-                    self.assertAlmostEqual(difference, 0.0)
+                # full run
+                m1 = model()
+                m1.params["duration"] = duration
+                pars_bak = copy.deepcopy(m1.params)
+                m1.run()
+                # chunkwise run
+                m2 = model()
+                m2.params = pars_bak.copy()
+                m2.run(chunkwise=True, chunksize=chunksize, append=True)
+                # check
+                self.assertTupleEqual(m1.output.shape, m2.output.shape)
+                difference = np.sum(np.abs(m1.output - m2.output))
+                print(
+                    f"single node model: {model.name}, duration: {duration}, chunksize: {chunksize}, difference: {difference}"
+                )
+                self.assertAlmostEqual(difference, 0.0)
 
     def network_test(self, model):
         ds = Dataset("hcp")
@@ -57,6 +60,9 @@ class AutochunkTests(unittest.TestCase):
                     # check
                     self.assertTupleEqual(m1.output.shape, m2.output.shape)
                     difference = np.sum(np.abs(m1.output - m2.output))
+                    print(
+                        f"network model: {model.name}, duration: {duration}, chunksize: {chunksize}, difference: {difference}"
+                    )
                     self.assertAlmostEqual(difference, 0.0)
 
 
@@ -111,8 +117,16 @@ class TestWCAutochunk(AutochunkTests):
         self.network_test(WCModel)
 
 
+class TestWWAutochunk(AutochunkTests):
+    def test_single(self):
+        self.single_node_test(WWModel)
+
+    def test_network(self):
+        self.network_test(WWModel)
+
+
 class TestThalamusAutochunk(AutochunkTests):
-    @pytest.mark.xfail
+    # @pytest.mark.xfail
     def test_single(self):
         self.single_node_test(ThalamicMassModel)
 
diff --git a/tests/test_models.py b/tests/test_models.py
index 47c908f..e60762c 100644
--- a/tests/test_models.py
+++ b/tests/test_models.py
@@ -13,6 +13,7 @@ from neurolib.models.multimodel import MultiModel
 from neurolib.models.multimodel.builder.fitzhugh_nagumo import FitzHughNagumoNetwork, FitzHughNagumoNode
 from neurolib.models.thalamus import ThalamicMassModel
 from neurolib.models.wc import WCModel
+from neurolib.models.ww import WWModel
 from neurolib.utils.collections import star_dotdict
 from neurolib.utils.loadData import Dataset
 
@@ -134,11 +135,11 @@ class TestWC(unittest.TestCase):
     def test_single_node(self):
         logging.info("\t > WC: Testing single node ...")
         start = time.time()
-        wc = WCModel()
-        wc.params["duration"] = 2.0 * 1000
-        wc.params["sigma_ou"] = 0.03
+        model = WCModel()
+        model.params["duration"] = 2.0 * 1000
+        model.params["sigma_ou"] = 0.03
 
-        wc.run()
+        model.run()
 
         end = time.time()
         logging.info("\t > Done in {:.2f} s".format(end - start))
@@ -147,14 +148,48 @@ class TestWC(unittest.TestCase):
         logging.info("\t > WC: Testing brain network (chunkwise integration and BOLD simulation) ...")
         start = time.time()
         ds = Dataset("gw")
-        wc = WCModel(Cmat=ds.Cmat, Dmat=ds.Dmat)
-        wc.params["signalV"] = 4.0
-        wc.params["duration"] = 10 * 1000
-        wc.params["sigma_ou"] = 0.1
-        wc.params["K_gl"] = 0.6
-        wc.params["x_ext_mean"] = 0.72
-
-        wc.run(chunkwise=True, bold=True, append_outputs=True)
+        model = WCModel(Cmat=ds.Cmat, Dmat=ds.Dmat)
+        model.params["signalV"] = 4.0
+        model.params["duration"] = 10 * 1000
+        model.params["sigma_ou"] = 0.1
+        model.params["K_gl"] = 0.6
+
+        # local node input parameter
+        model.params["exc_ext"] = 0.72
+
+        model.run(chunkwise=True, bold=True, append_outputs=True)
+        end = time.time()
+        logging.info("\t > Done in {:.2f} s".format(end - start))
+
+
+class TestWW(unittest.TestCase):
+    """
+    Basic test for WW model.
+    """
+
+    def test_single_node(self):
+        logging.info("\t > WW: Testing single node ...")
+        start = time.time()
+        model = WWModel()
+        model.params["duration"] = 2.0 * 1000
+        model.params["sigma_ou"] = 0.03
+
+        model.run()
+
+        end = time.time()
+        logging.info("\t > Done in {:.2f} s".format(end - start))
+
+    def test_network(self):
+        logging.info("\t > WW: Testing brain network (chunkwise integration and BOLD simulation) ...")
+        start = time.time()
+        ds = Dataset("gw")
+        model = WWModel(Cmat=ds.Cmat, Dmat=ds.Dmat)
+        model.params["signalV"] = 4.0
+        model.params["duration"] = 10 * 1000
+        model.params["sigma_ou"] = 0.1
+        model.params["K_gl"] = 0.6
+
+        model.run(chunkwise=True, bold=True, append_outputs=True)
         end = time.time()
         logging.info("\t > Done in {:.2f} s".format(end - start))