Wilson-Cowan
The Wilson-Cowan model of interacting excitatory and inhibitory populations of neurons.
- Name:
wilson_cowan - Field:
biology/neuroscience - Formalism: ordinary differential equations
- References: Wilson and Cowan 1972; arXiv:2112.12147
Factory
@register_model(
field='biology/neuroscience',
tags=('neural mass', 'oscillations', 'excitation-inhibition'),
references=('Wilson and Cowan 1972', 'arXiv:2112.12147'),
defaults={'until': 100.0, 'dt': 0.2},
)
def wilson_cowan(
E: float = 0.1,
I: float = 0.05,
*,
P: float = 1.25,
Q: float = 0.0,
c1: float = 16.0,
c2: float = 12.0,
c3: float = 15.0,
c4: float = 3.0,
) -> System:
"""Return the Wilson-Cowan model of interacting excitatory and inhibitory populations of neurons.
With the published couplings, a moderate excitatory input `P = 1.25` gives a limit cycle of population activity,
and no input lets the activity decay to rest.
Args:
E:
The initial excitatory activity.
I:
The initial inhibitory activity.
P:
The external input to the excitatory population.
Q:
The external input to the inhibitory population.
c1:
The excitatory to excitatory coupling.
c2:
The inhibitory to excitatory coupling.
c3:
The excitatory to inhibitory coupling.
c4:
The inhibitory to inhibitory coupling.
Returns:
system:
The model, with the activities `cortex.E` and `cortex.I`.
Examples:
>>> from skcomplex.models.biology.neuroscience import wilson_cowan
>>> wilson_cowan(P=0.0).elements[0].P
0.0
"""
cortex = WilsonCowanElement('cortex', E=E, I=I, P=P, Q=Q, c1=c1, c2=c2, c3=c3, c4=c4)
relaxations = [
RelaxationInteraction(
'excitatory_relaxation',
target='cortex',
state='E',
goal=_excitation,
time=lambda s: s.target.tau_e,
),
RelaxationInteraction(
'inhibitory_relaxation',
target='cortex',
state='I',
goal=_inhibition,
time=lambda s: s.target.tau_i,
),
]
return System(elements=[cortex], interactions=relaxations)
Components
class WilsonCowanElement(Element):
"""Excitatory and inhibitory populations of Wilson and Cowan, as fractions of active cells.
The excitatory activity follows `tau_e E' = -E + (k_e - r_e E) S_e(c1 E - c2 I + P)` and the inhibitory activity
`tau_i I' = -I + (k_i - r_i I) S_i(c3 E - c4 I + Q)`, where `S(x) = 1 / (1 + exp(-a (x - theta))) - 1 / (1 +
exp(a theta))` is a sigmoid shifted to vanish at zero input and `k` is its largest value. The element holds the
states and parameters and has no laws of its own: the relaxations of `wilson_cowan` read its parameters.
Args:
label:
The label of the populations.
visual:
How the element is drawn, or `None` for `default_visual_`.
E:
The initial excitatory activity.
I:
The initial inhibitory activity.
c1:
The excitatory to excitatory coupling.
c2:
The inhibitory to excitatory coupling.
c3:
The excitatory to inhibitory coupling.
c4:
The inhibitory to inhibitory coupling.
a_e:
The gain of the excitatory response.
theta_e:
The threshold of the excitatory response.
a_i:
The gain of the inhibitory response.
theta_i:
The threshold of the inhibitory response.
r_e:
The refractory period of excitatory cells.
r_i:
The refractory period of inhibitory cells.
P:
The external input to the excitatory population.
Q:
The external input to the inhibitory population.
tau_e:
The time constant of the excitatory population.
tau_i:
The time constant of the inhibitory population.
Examples:
>>> from skcomplex.models.biology.neuroscience import WilsonCowanElement
>>> WilsonCowanElement('cortex', P=1.25).c1
16.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
E: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
I: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
c1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 16.0,
c2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 12.0,
c3: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 15.0,
c4: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 3.0,
a_e: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.3,
theta_e: Annotated[float, ParameterSpec()] = 4.0,
a_i: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
theta_i: Annotated[float, ParameterSpec()] = 3.7,
r_e: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
r_i: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
P: Annotated[float, ParameterSpec()] = 1.25,
Q: Annotated[float, ParameterSpec()] = 0.0,
tau_e: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
tau_i: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.E = E
self.I = I
self.c1 = c1
self.c2 = c2
self.c3 = c3
self.c4 = c4
self.a_e = a_e
self.theta_e = theta_e
self.a_i = a_i
self.theta_i = theta_i
self.r_e = r_e
self.r_i = r_i
self.P = P
self.Q = Q
self.tau_e = tau_e
self.tau_i = tau_i
class RelaxationInteraction(RateInteraction):
"""Relaxation `(goal(ns) - x) / time(ns)` of a state `x` of a target towards a goal, with optional noise.
With a constant goal and time and an additive noise of amplitude `sigma`, the state follows the Ornstein-Uhlenbeck
process `dx = (goal - x) / time dt + sigma dW`, whose stationary variance is `sigma**2 time / 2`. Without noise it
is first-order smoothing, or adaptive expectation, towards the goal. A goal or a time function takes the namespace
of the interaction, as for `RateInteraction`, and a number stands for a constant. A goal that is a number applies
to every component of a vector state.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state relaxes.
state:
The name of the state of the target that relaxes.
goal:
A function of the namespace that returns the goal, or a number.
time:
A function of the namespace that returns the relaxation time, or a positive number.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
Raises:
ValidationError:
The state is not a name, or a function is neither callable nor a number.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import RelaxationInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> smoothing = RelaxationInteraction('smoothing', target='belief', state='n', goal=10.0, time=0.5)
>>> system = System(elements=[CompartmentElement('belief', n=0.0)], interactions=[smoothing])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=5.0).trajectory_
>>> round(float(trajectory.values('belief.n')[-1]), 3)
10.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'x',
goal: Callable[[SimpleNamespace], object] | float = 0.0,
time: Callable[[SimpleNamespace], object] | float = 1.0,
noise: Callable[[SimpleNamespace], object] | float | None = None,
) -> None:
super().__init__(label, target=target, state=state, noise=noise)
self.goal = goal
self.time = time
def _check_functions(self) -> None:
"""Raise an error for a state that is not a name, or a function that is neither callable nor a number."""
super()._check_functions()
for name in ('goal', 'time'):
value = getattr(self, name)
if not (callable(value) or isinstance(value, int | float)):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace or a number.'
raise ValidationError(error_msg)
def _get_goal(self, namespace: SimpleNamespace) -> object:
"""Return the goal, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.goal, namespace)
def _get_time(self, namespace: SimpleNamespace) -> object:
"""Return the relaxation time, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.time, namespace)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the relaxation `(goal - x) / time` of the state.
Args:
namespace:
The namespace of the symbols the rate can use.
Returns:
rate:
The rate of the state.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
goal = self._broadcast(self._get_goal(namespace), state)
return (goal - state) / self._get_time(namespace)
Equations
cortex.E: diff(E_cortex(t), t) == -((r_e_cortex*E_cortex(t) + 1/(e^(a_e_cortex*theta_e_cortex) + 1) - 1)*(1/(e^(-(c1_cortex*E_cortex(t) - c2_cortex*I_cortex(t) + P_cortex - theta_e_cortex)*a_e_cortex) + 1) - 1/(e^(a_e_cortex*theta_e_cortex) + 1)) + E_cortex(t))/tau_e_cortex
cortex.I: diff(I_cortex(t), t) == -((r_i_cortex*I_cortex(t) + 1/(e^(a_i_cortex*theta_i_cortex) + 1) - 1)*(1/(e^(-(c3_cortex*E_cortex(t) - c4_cortex*I_cortex(t) + Q_cortex - theta_i_cortex)*a_i_cortex) + 1) - 1/(e^(a_i_cortex*theta_i_cortex) + 1)) + I_cortex(t))/tau_i_cortex