Note
Click here to download the full example code
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)
- 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):
where represents the average firing rate of the population, :math:` au` is a lumped population time constant,
is a coupling term,
is a refractory term, and
is the synaptic input.
Typically,
is chosen as a sigmoidal transfer function that transforms the pre-synaptic membrane potential
into an average firing rate:
Considering both the E and the I population, the Wilson-Cowan model can be described via two coupled ODEs:
where additional extrinsic inputs and
have been added.
Below, we will demonstrate how to load this a model into pyrates and perform numerical simulations with it.
References¶
[1] | H.R. Wilson and J.D. Cowan (1972) *Excitatory and Inhibitory Interactions in Localized Populations of Model |
Neurons.* Biophysical Journal 12.
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
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")
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')
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 and
of the model.
import matplotlib.pyplot as plt
plt.plot(results)
plt.show()
Total running time of the script: ( 0 minutes 0.000 seconds)