diff --git a/examples/example-0.4-wc-minimal.ipynb b/examples/example-0.4-wc-minimal.ipynb index 76b1ddbc..62eb74e9 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 1de22bac..9c7e8985 100644 --- a/neurolib/models/multimodel/builder/wilson_cowan.py +++ b/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]]) diff --git a/neurolib/models/wc/loadDefaultParams.py b/neurolib/models/wc/loadDefaultParams.py index 34726ecd..788dfe0d 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 67af8c2b..67dfee94 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 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 ) diff --git a/neurolib/models/ww/__init__.py b/neurolib/models/ww/__init__.py new file mode 100644 index 00000000..576f30f5 --- /dev/null +++ b/neurolib/models/ww/__init__.py @@ -0,0 +1 @@ +from .model import WWModel diff --git a/neurolib/models/ww/loadDefaultParams.py b/neurolib/models/ww/loadDefaultParams.py new file mode 100644 index 00000000..05d68108 --- /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 00000000..12d72ebc --- /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 00000000..fcaeabf9 --- /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 bdd7a678..1bd4e08e 100644 --- a/setup.py +++ b/setup.py @@ -11,7 +11,7 @@ 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", diff --git a/tests/test_autochunk.py b/tests/test_autochunk.py index 99a69098..87a24784 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.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 @@ def network_test(self, model): # 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 @@ def test_network(self): 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 47c908f4..e60762c7 100644 --- a/tests/test_models.py +++ b/tests/test_models.py @@ -13,6 +13,7 @@ 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 @@ def test_network(self): 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))