2.2. The Wilson-Cowan Neural Mass Model

Here, we will introduce the Wilson-Cowan model, a popular neural mass model of the dynamic interactions between 2 populations: excitatory projection neurons (E), and inhibitory interneurons (I).

The model is of phenomenological nature, but captures well the steady-states in average firing rate of populations of leaky integrate-and-fire neurons [1]. It is probably the most-used existing neural mass model and has been applied to a wide range of brain systems and neuroscientific questions. The dynamics of each population are described by a 1st-order ordinary differential equation (ODE):

\[\tau \dot r = -r + (k - q r) S(m),\]

where \(r\) represents the average firing rate of the population, :math:` au` is a lumped population time constant, \(k\) is a coupling term, \(q\) is a refractory term, and \(S(m)\) is the synaptic input. Typically, \(S\) is chosen as a sigmoidal transfer function that transforms the pre-synaptic membrane potential into an average firing rate:

\[S(m) = \frac{1}{1 + e^{-m}}.\]

Considering both the E and the I population, the Wilson-Cowan model can be described via two coupled ODEs:

\[ \begin{align}\begin{aligned}\tau_e \dot r_e &= -r_e + (k_e - q_e r_e) S(c_{ee} r_e + c_{ei} r_i + I_e(t)),\\\tau_i \dot r_i &= -r_i + (k_i - q_i r_i) S(c_{ie} r_e + c_{ii} r_i + I_i(t))\end{aligned}\end{align} \]

where additional extrinsic inputs \(I_e\) and \(I_i\) have been added. Below, we will demonstrate how to load this a model into pyrates and perform numerical simulations with it.

References

2.2.1. Step 1: Importing the frontend class for defining models

As a first step, we import the pyrates.frontend.CircuitTemplate class, which allows us to set up a model definition in PyRates.

from pyrates.frontend import CircuitTemplate

2.2.2. Step 2: Loading a model template from the model_templates library

In the second step, we load a model template for the Wilson-Cowan model that comes with PyRates via the from_yaml() method of the CircuitTemplate. This method returns a CircuitTemplate instance. Have a look at the yaml definition of the model that can be found via the path used for the from_yaml() method. You will see that all variables and parameters are already defined there. As an alternative, you can also use the circuit_from_yaml() function that you can simply import via from pyrates import circuit_from_yaml. These are the basic steps you perform, if you want to load a model that is defined inside a yaml file. To check out the different model templates provided by PyRates, have a look at the PyRates.model_templates module.

wc = CircuitTemplate.from_yaml("model_templates.neural_mass_models.wilsoncowan.WC")

2.2.3. Step 3: Numerical simulation of a the model behavior in time

After loading the model template, numerical simulations can be performed via the run() method. Calling this function will solve the initial value problem of the above defined differential equations for a time interval from 0 to the given simulation time. This solution will be calculated numerically by a differential equation solver in the backend, starting with a defined step-size. As part of this process, a pyrates.backend.BaseBackend instance is created that contains a graph representation of the network equations. How exactly these equations will be solved, can be customized via the arguments of the run() method. You can choose between different backends, numerical solvers and integration step-sizes, for instance. Here, we choose the scipy backend, which uses a 4(5)th order Runge-Kutta method as default numerical solver. The default of the run() method is the forward Euler method that is implemented in PyRates itself.

import numpy as np

# simulation parameters
T = 125.0
dt = 5e-4
dts = 1e-2
steps = int(np.round(T/dt))
inp = np.zeros((steps,))
inp[int(25/dt):int(100/dt)] = 1.0

# perform simulation
results = wc.run(simulation_time=T,
                 step_size=dt,
                 sampling_step_size=dts,
                 outputs={'E': 'e/rate_op/r',
                          'I': 'i/rate_op/r'},
                 inputs={'e/se_op/r_ext': inp},
                 backend='default',
                 solver='euler')

2.2.4. Step 4: Visualization of the solution

The output of the run() method is a pandas.Dataframe, which comes with a plot() method for plotting the timeseries it contains. This timeseries represents the numerical solution of the initial value problem solved in step 4 with respect to the state variables \(r_e\) and \(r_i\) of the model.

import matplotlib.pyplot as plt
plt.plot(results)
plt.show()

Gallery generated by Sphinx-Gallery