"""Wilson‑Cowan neural population dynamics with empirical parameters."""
from typing import Optional, Tuple
import numpy as np
from mercurial.params.empirical import WilsonCowanParams
[docs]
class WilsonCowanPopulation:
"""
Wilson‑Cowan population with empirical parameters from literature.
"""
[docs]
def __init__(self, params: Optional[WilsonCowanParams] = None, **kwargs):
"""
Parameters
----------
params : WilsonCowanParams, optional
Empirical parameter dataclass. If None, uses default.
**kwargs : individual parameter overrides.
"""
if params is None:
params = WilsonCowanParams()
# Set parameters from dataclass, then override with kwargs
self.tau_e = kwargs.get("tau_e", params.tau_e)
self.tau_i = kwargs.get("tau_i", params.tau_i)
self.r_e = kwargs.get("r_e", 0.8)
self.r_i = kwargs.get("r_i", 0.8)
self.w_ee = kwargs.get("w_ee", params.w_ee)
self.w_ei = kwargs.get("w_ei", params.w_ei)
self.w_ie = kwargs.get("w_ie", params.w_ie)
self.w_ii = kwargs.get("w_ii", params.w_ii)
self.a_e = kwargs.get("a_e", params.a_e)
self.theta_e = kwargs.get("theta_e", params.mu_e)
self.a_i = kwargs.get("a_i", params.a_i)
self.theta_i = kwargs.get("theta_i", params.mu_i)
self.sigma = kwargs.get("sigma", 0.01)
self.rng = np.random.default_rng()
# State variables
self.E = 0.1
self.I = 0.1
[docs]
def sigmoid(self, x: float, a: float, theta: float) -> float:
arg = a * (x - theta)
# Clip to avoid overflow in exp
if arg > 500:
return 1.0
if arg < -500:
return 0.0
return 1.0 / (1.0 + np.exp(-arg))
[docs]
def derivative(
self, t: float, state: np.ndarray, P_ext: float = 0.0, Q_ext: float = 0.0
) -> np.ndarray:
E, I = state
input_e = self.w_ee * E - self.w_ei * I + P_ext
S_e = self.sigmoid(input_e, self.a_e, self.theta_e)
dE = (-E + (1 - self.r_e * E) * S_e) / self.tau_e
input_i = self.w_ie * E - self.w_ii * I + Q_ext
S_i = self.sigmoid(input_i, self.a_i, self.theta_i)
dI = (-I + (1 - self.r_i * I) * S_i) / self.tau_i
return np.array([dE, dI])
[docs]
def clamped_derivative(
self, t: float, state: np.ndarray, P_ext: float = 0.0, Q_ext: float = 0.0
) -> np.ndarray:
"""Derivative with clamping to prevent negative activity."""
deriv = self.derivative(t, state, P_ext, Q_ext)
E, I = state
dE, dI = deriv
if E + dE * 1e-3 < 0:
dE = -E / 1e-3
if I + dI * 1e-3 < 0:
dI = -I / 1e-3
return np.array([dE, dI])
[docs]
def step(self, dt: float, P_ext: float = 0.0, Q_ext: float = 0.0) -> None:
"""Euler step with noise."""
dE, dI = self.clamped_derivative(0, np.array([self.E, self.I]), P_ext, Q_ext)
self.E += dE * dt + self.sigma * np.sqrt(dt) * self.rng.normal()
self.I += dI * dt + self.sigma * np.sqrt(dt) * self.rng.normal()
self.E = np.clip(self.E, 0.0, 1.0)
self.I = np.clip(self.I, 0.0, 1.0)
[docs]
def evolve(
self,
dt: float,
n_steps: int,
P_ext: float = 0.0,
Q_ext: float = 0.0,
initial_state: Optional[np.ndarray] = None,
) -> Tuple[np.ndarray, np.ndarray]:
"""
Evolve the population for n_steps using the step method.
Returns (E_history, I_history).
"""
if initial_state is not None:
self.E, self.I = initial_state
E_hist = np.zeros(n_steps + 1)
I_hist = np.zeros(n_steps + 1)
E_hist[0] = self.E
I_hist[0] = self.I
for i in range(n_steps):
self.step(dt, P_ext, Q_ext)
E_hist[i + 1] = self.E
I_hist[i + 1] = self.I
return E_hist, I_hist
[docs]
def steady_state(
self, P_ext: float = 0.0, Q_ext: float = 0.0, tol: float = 1e-6, max_iter: int = 1000
) -> Tuple[float, float]:
"""
Find fixed point by iteration.
Returns (E_ss, I_ss).
"""
E, I = 0.1, 0.1
for _ in range(max_iter):
input_e = self.w_ee * E - self.w_ei * I + P_ext
S_e = self.sigmoid(input_e, self.a_e, self.theta_e)
E_new = (1 - self.r_e * E) * S_e
input_i = self.w_ie * E - self.w_ii * I + Q_ext
S_i = self.sigmoid(input_i, self.a_i, self.theta_i)
I_new = (1 - self.r_i * I) * S_i
if abs(E_new - E) < tol and abs(I_new - I) < tol:
return E_new, I_new
E, I = E_new, I_new
return E, I