Simulation module

Simulations utilities

Simulate MEEG-like signals with different connectivity patterns or methods.

Todo

• Follow up on these:

  • RNN of rate models: RNN of rate models

  • ctRNN revision, see ctRNN revision

• Simulation based on connectivity matrix: implement a network of N nodes simply using the connectivity matrix and node-specific implementations

• Neural mass models to be added:

  • Wilson-Cowan

class pyeeg.simulate.CTRNN(N, W, input_dim=1, output_dim=1, dt=0.001, seed=42, nonlinearity=<function sigmoid>, theta=None)

Continuous Time Recurrent Neural Network (CTRNN) model.

\[\tau \dot{x} = -x + W o + I + \theta\]

simulate(x0, tmax=1, noise=0.0, I=<function CTRNN.<lambda>>)

Simulate the CTRNN model and monitor the output.

Parameters:

• x0 (array_like) – The initial state of the system. Shape (N,).

• tmax (float) – The maximum time to simulate.

• noise (float) – The standard deviation of the noise to add to the system.

Returns:

x – The simulated time series. Shape (n, N).

Return type:

array_like

step(I=None, noise=0.0)

Compute one step of the CTRNN model.

class pyeeg.simulate.JRNetwork(N=2, W=array([[0, 1], [0, 0]]), delay=0.01, w=0.8, node_dynamics=None, dt=0.001, seed=42)

Abstract class for neural mass models. Defines the function to be implemented for the simulation.

Notes

Two types of networks are modelled in the literature: either one when the mean input and variance of the input are controlled for each node, such that the input received from connected nodes is normalised and will relatively shut down contribution from external input, or one where the input is not normalised and the external input is simply summed over the input from connected nodes. The latter is seen in [2] and [3], while the former is seen in [4] & [5].

• [1] Jansen, B. H., & Rit, V. G. (1995). Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biological cybernetics, 73(4), 357-366.

• [2] Kazemi & Jamali, (2022). Phase synchronization and measure of criticality in a network of neural mass models. Nature.

• [3] Forrester et al. (2020). Network Neuroscience. The role of node dynamics in shaping emergent functional connectivity patterns in the brain. PMC.

• [4] David & Friston (2006). NeuroImage. A Neural mass model for MEG/EEG: coupling and neuronal dynamics. ScienceDirect.

• [5] David et al., (2004). Evaluation of different measures of functional connectivity using a neural mass model. ScienceDirect.

[1]

Jansen, B. H., & Rit, V. G. (1995). Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biological cybernetics, 73(4), 357-366.

[2] (1,2)

Kazemi & Jamali, (2022). Phase synchronization and measure of criticality in a network of neural mass models. Nature.

[3] (1,2)

Forrester et al. (2020). Network Neuroscience. The role of node dynamics in shaping emergent functional connectivity patterns in the brain. PMC.

[4] (1,2)

David & Friston (2006). NeuroImage. A Neural mass model for MEG/EEG: coupling and neuronal dynamics. ScienceDirect.

[5] (1,2)

David et al., (2004). Evaluation of different measures of functional connectivity using a neural mass model. ScienceDirect.

update_connectivity(x, sigma_p=1)

x represents the firing rates output of all nodes, needed to compute the standard deviation Shape of x: (N, ntimes)

class pyeeg.simulate.JansenRit(dt=0.0001, seed=42, nonlinearity=<function sigmoid>)

Jansen-Rit model.

3 populations: excitatory, inhibitory and pyramidal:

___________    ___________
|         |    |         |
!  Inhib  !    !  Excit  !
|         |    |         |
-----------    -----------
C2,C4 \             / C1, C3
       \___________/
        |         |
        ! Pyramid !
        |         |
        -----------

Parameters and typical values as in Grimbert & Faugeras, 2006:

• C1, C2, C3, C4: Average number of synapses between populations: 135 * [1, 0.8, 0.25, 0.25]

• tau_e: Time scale for excitatory population: 100 ms

• tau_i: Time scale for inhibitory population: 50 ms

• G_exc: Average excitatory synaptic gain: 3.25

• G_inh: Average inhibitory synaptic gain: 22

• rmax: Amplitude of sigmoid: 5 s^-1

• beta: Slope of sigmoid: 0.56 mV^-1

• theta: Threshold of sigmoid: 6 mV

• Conduction velocity: 10 m/s

• h: Integration time step: 0.0001 s (by default)

• P: External input to each of the neural masses: 150 Hz (constant input)

• Coupling: Coupling between the neural masses: [0.1:0.012:0.292]

This table is from the paper: Kulik et al, Network Neurosci. (2023)

simulate(x0, tmax=1, noise=0.0, P=None)

Simulate the Jansen-Rit model and monitor the output.

Parameters:

• x0 (array_like) – The initial state of the system. Shape (N,).

• tmax (float) – The maximum time to simulate.

• noise (float) – The standard deviation of the noise to add to the system.

Returns:

x – The simulated time series. Shape (n, N).

Return type:

array_like

step(I=0.0)

Compute one step of the Jansen-Rit model.

class pyeeg.simulate.JansenRitExtended(w=0.5, dt=0.0001, seed=42, nonlinearity=<function sigmoid>)

Jansen-Rit model - Exctended version: dual kinetic model.

See David & Friston, NeuroImage (2003)

We model two parallel subpopulations with different kinematics in order to capture multiband or broadband dynamics.

simulate(x0, tmax=1, noise=0.0, P=None)

Simulate the Jansen-Rit model and monitor the output.

Parameters:

• x0 (array_like) – The initial state of the system. Shape (N,).

• tmax (float) – The maximum time to simulate.

• noise (float) – The standard deviation of the noise to add to the system.

Returns:

x – The simulated time series. Shape (n, N).

Return type:

array_like

step(I=0.0)

Compute one step of the Jansen-Rit model.

class pyeeg.simulate.NeuralMassNetwork(N, W, delay=0, node_dynamics=None, dt=0.001, seed=42)

Abstract class for neural mass models. Defines the function to be implemented for the simulation.

class pyeeg.simulate.NeuralMassNode(dt=0.001, seed=42)

Abstract class for neural mass models. Defines the function to be implemented for the simulation.

pyeeg.simulate.simulate_ar(order, coefs, n, sigma=1, seed=42)

Simulate an autoregressive process of order order.

Parameters:

• order (int) – The order of the autoregressive process.

• coefs (array_like) – The coefficients of the autoregressive process. The first element is the coefficient of the lag (t-1).

• n (int) – The number of samples to simulate.

• sigma (float) – The standard deviation of the additive noise process.

Returns:

x – The simulated time series. Shape (n,).

Return type:

array_like

pyeeg.simulate.simulate_var(order, coef, nobs=500, ndim=2, seed=42, verbose=False)

Simulate a VAR model of order order.

The VAR model is defined as:

\[x_t = A_1 x_{t-1} + A_2 x_{t-2} + ... + A_p x_{t-p} + \epsilon_t x_t = \sum_{i=1}^p A_i x_{t-i} + \epsilon_t\]

where \(x_t\) is a vector of shape (ndim, 1), \(A_i\) is a matrix of shape (ndim, ndim)

Note

The coefficients at a given lag are such as \(C_ij\) is i->j, so it will be the coefficients for dimension j! For example, each row of the first column are determining the contributions of each component onto the first component.

pyeeg.simulate.simulate_var_from_cov(cov, nobs=500, ndim=2, seed=42, verbose=False)

Simulate a VAR model of order order from a covariance matrix.

The VAR model is defined as:

\[x_t = A_1 x_{t-1} + A_2 x_{t-2} + ... + A_p x_{t-p} + \epsilon_t x_t = \sum_{i=1}^p A_i x_{t-i} + \epsilon_t\]

where \(x_t\) is a vector of shape (ndim, 1), \(A_i\) is a matrix of shape (ndim, ndim)

Note

The coefficients at a given lag are such as \(C_ij\) is i->j, so it will be the coefficients for dimension j! For example, each row of the first column are determining the contributions of each component onto the first component.

Some Simulation Models:

Classes