Jansen-Rit epileptiform regime

def _stable_expit(x):
    """Numerically stable sigmoid 1/(1+exp(-x)), scalar or array."""
    x = np.asarray(x, dtype=float)
    return np.where(
        x >= 0,
        1.0 / (1.0 + np.exp(-x)),
        np.exp(x) / (1.0 + np.exp(x)),
    )

def _jr_sigmoid(v, v_max=5.0, r=0.56, v0=6.0):
    """Jansen-Rit population sigmoid: S(v) = v_max / (1 + exp(r*(v0 - v)))."""
    return v_max * _stable_expit(r * (v - v0))

def rhs(t: float, y: np.ndarray) -> np.ndarray:
    x0, x1, x2, x3, x4, x5 = y
    S_pyr = _jr_sigmoid(x1 - x2)
    S_exc = _jr_sigmoid(C1 * x0)
    S_inh = _jr_sigmoid(C3 * x0)
    return np.array([
        x3,
        x4,
        x5,
        A * a * S_pyr - 2.0 * a * x3 - a2 * x0,
        A * a * (p + C2 * S_exc) - 2.0 * a * x4 - a2 * x1,
        B * b * C4 * S_inh - 2.0 * b * x5 - b2 * x2,
    ])