Wilson–Cowan model

Proportion of cells in refractory period (absolute refractory period ${\displaystyle r}$) ${\displaystyle \int _{t-r}^{t}E(t')dt'}$

Proportion of sensitive cells (complement of refractory cells) ${\displaystyle 1-\int _{t-r}^{t}E(t')dt'}$

If ${\displaystyle \theta }$ denotes a cell's threshold potential and ${\displaystyle D(\theta )}$ is the distribution of thresholds in the tissue, then the expected proportion of neurons receiving an excitation at or above threshold level per unit time is

${\displaystyle S({\bar {N}})=\int _{0}^{{\bar {N}}(t)}D(\theta )d\theta }$,

where ${\displaystyle {\bar {N}}(t)}$ is the mean integrated excitation at time t. The term "integrated" in this case means that each neuron sums up (i.e. integrates) all incoming excitations in a linear fashion to receive its total excitation. If this integrated excitation is at or above the neuron's excitation threshold, it will in turn create an action potential. Note that the above equation relies heavily on the homogeneous distribution of neurons, as does the Wilson-Cowan model in general.

Subpopulation response function based on the distribution of afferent synapses per cell (all cells have the same threshold)

${\displaystyle S({\bar {N}})=\int _{\frac {\theta }{{\bar {N}}(t)}}^{\infty }C(w)dw}$

Average excitation level of an excitatory cell at time ${\displaystyle t}$

${\displaystyle \int _{-\infty }^{t}\alpha (t-t')[c_{1}E(t')-c_{2}I(t')+P(t')]dt'}$

where ${\displaystyle \alpha (t)}$ is the stimulus decay function, ${\displaystyle c_{1}}$ and ${\displaystyle c_{2}}$ are respectively the connectivity coefficient giving the average number of excitatory and inhibitory synapses per cell, P(t) is the external input to the excitatory population.

Excitatory subpopulation expression

${\displaystyle E(t)=[1-\int _{t-r}^{t}E(t')dt']S({\bar {N}})dt}$

Complete Wilson–Cowan model

${\displaystyle E(t+\tau )=[1-\int _{t-r}^{t}E(t')dt']S_{e}\left\{\int _{-\infty }^{t}\alpha (t-t')[c_{1}E(t')-c_{2}I(t')+P(t')]dt'\right\}}$

${\displaystyle I(t+\tau )=[1-\int _{t-r}^{t}I(t')dt']S_{i}\left\{\int _{-\infty }^{t}\alpha (t-t')[c_{3}E(t')-c_{4}I(t')+Q(t')]dt'\right\}}$

Time Coarse Graining ${\displaystyle \tau {\frac {d{\bar {E}}}{dt}}=-{\bar {E}}+(1-r{\bar {E}})S_{e}[kc_{1}{\bar {E}}(t)+kP(t)]}$

Isocline Equation ${\displaystyle c_{2}I=c_{1}E-S_{e}^{-1}\left({\frac {E}{k_{e}-r_{e}E}}\right)+P}$

Sigmoid Function ${\displaystyle S(x)={\frac {1}{1+\exp[-a({\bar {N}}-\theta )]}}-{\frac {1}{1+\exp(a\theta )}}}$

The transition from normal state of brain activity to epileptic seizures was not formulated theoretically until 2008, when a theoretical path from a baseline state to large-scale self-sustained oscillations, which spread out uniformly from the point of stimulus, has been mapped for the first time.[7]

A realistic state of baseline physiological activity has been defined, using the following two-component definition:[7]

(1) A time-independent component represented by subthreshold excitatory activity E and superthreshold inhibitory activity I.

(2) A time-varying component which may include singlepulse waves, multipulse waves, or periodic waves caused by spontaneous neuronal activity.

This baseline state represents activity of the brain in the state of relaxation, in which neurons receive some level of spontaneous, weak stimulation by small, naturally present concentrations of neurohormonal substances. In waking adults this state is commonly associated with alpha rhythm, whereas slower (theta and delta) rhythms are usually observed during deeper relaxation and sleep. To describe this general setting, a 3-variable ${\displaystyle (u,I,v)}$ spatially dependent extension of the classical Wilson–Cowan model can be utilized.[10] Under appropriate initial conditions,[7] the excitatory component, u, dominates over the inhibitory component, I, and the three-variable system reduces to the two-variable Pinto-Ermentrout type model[11]

${\displaystyle {\partial u \over \partial t}=u-v+\int _{R^{2}}\omega (x-x',y-y')f(u-\theta )\,dxdy+\zeta (x,y,t),}$

${\displaystyle {\partial v \over \partial t}=\epsilon (\beta u-v).}$

The variable v governs the recovery of excitation u; ${\displaystyle \epsilon >0}$ and ${\displaystyle \beta >0}$ determine the rate of change of recovery. The connection function ${\displaystyle \omega (x,y)}$ is positive, continuous, symmetric, and has the typical form ${\displaystyle \omega =Ae^{-\lambda {\sqrt {-(x^{2}+y^{2})}}}}$.[11] In Ref.[7] ${\displaystyle (A,\lambda )=(2.1,1).}$ The firing rate function, which is generally accepted to have a sharply increasing sigmoidal shape, is approximated by ${\displaystyle f(u-\theta )=H(u-\theta )}$, where H denotes the Heaviside function; ${\displaystyle \zeta (x,y,t)}$ is a short-time stimulus. This ${\displaystyle (u,v)}$ system has been successfully used in a wide variety of neuroscience research studies.[11][12][13][14][15] In particular, it predicted the existence of spiral waves, which can occur during seizures; this theoretical prediction was subsequently confirmed experimentally using optical imaging of slices from the rat cortex.[16]

The expansion of hypersynchronized regions exhibiting large-amplitude stable bulk oscillations occurs when the oscillations coexist with the stable rest state ${\displaystyle (u,v)=(0,0)}$. To understand the mechanism responsible for the expansion, it is necessary to linearize the ${\displaystyle (u,v)}$ system around ${\displaystyle (0,0)}$ when ${\displaystyle \epsilon >0}$ is held fixed. The linearized system exhibits subthreshold decaying oscillations whose frequency increases as ${\displaystyle \beta }$ increases. At a critical value ${\displaystyle \beta ^{*}}$ where the oscillation frequency is high enough, bistability occurs in the ${\displaystyle (u,v)}$ system: a stable, spatially independent, periodic solution (bulk oscillation) and a stable rest state coexist over a continuous range of parameters. When ${\displaystyle \beta \geq \beta ^{*}}$ where bulk oscillations occur,[7] "The rate of expansion of the hypersynchronization region is determined by an interplay between two key features: (i) the speed c of waves that form and propagate outward from the edge of the region, and (ii) the concave shape of the graph of the activation variable u as it rises, during each bulk oscillation cycle, from the rest state u=0 to the activation threshold. Numerical experiments show that during the rise of u towards threshold, as the rate of vertical increase slows down, over time interval ${\displaystyle \Delta t,}$ due to the concave component, the stable solitary wave emanating from the region causes the region to expand spatially at a Rate proportional to the wave speed. From this initial observation it is natural to expect that the proportionality constant should be the fraction of the time that the solution is concave during one cycle." Therefore, when ${\displaystyle \beta \geq \beta ^{*}}$, the rate of expansion of the region is estimated by[7]

${\displaystyle Rate=(\Delta t/T)*c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(1)}$

where ${\displaystyle \Delta t}$ is the length of subthreshold time interval, T is period of the periodic solution; c is the speed of waves emanating from the hypersynchronization region. A realistic value of c, derived by Wilson et al.,[17] is c=22.4 mm/s.

How to evaluate the ratio ${\displaystyle \Delta t/T?}$ To determine values for ${\displaystyle \Delta t/T}$ it is necessary to analyze the underlying bulk oscillation which satisfies the spatially independent system

${\displaystyle {{du} \over {dt}}=u-v+H(u-\theta ),}$

${\displaystyle {{dv} \over {dt}}=\epsilon (\beta u-v).}$

This system is derived using standard functions and parameter values ${\displaystyle \omega =2.1e^{-\lambda {\sqrt {-(x^{2}+y^{2})}}}}$, ${\displaystyle \epsilon =0.1}$ and ${\displaystyle \theta =0.1}$[7][11][12][13] Bulk oscillations occur when ${\displaystyle \beta \geq \beta ^{*}=12.61}$. When ${\displaystyle 12.61\leq \beta \leq 17}$, Shusterman and Troy analyzed the bulk oscillations and found ${\displaystyle 0.136\leq \Delta t/T\leq 0.238}$. This gives the range

${\displaystyle 3.046mm/s\leq Rate\leq 5.331mm/s~~~~~~~~~~~~(2)}$

Since ${\displaystyle 0.136\leq \Delta t/T\leq 0.238}$, Eq. (1) shows that the migration Rate is a fraction of the traveling wave speed, which is consistent with experimental and clinical observations regarding the slow spread of epileptic activity.[18] This migration mechanism also provides a plausible explanation for spread and sustenance of epileptiform activity without a driving source that, despite a number of experimental studies, has never been observed.[18]

The rate of migration of hypersynchronous activity that was experimentally recorded during seizures in a human subject, using chronically implanted subdural electrodes on the surface of the left temporal lobe,[8] has been estimated as[7]

${\displaystyle Rate\approx 4mm/s}$,

which is consistent with the theoretically predicted range given above in (2). The ratio ${\displaystyle Rate/c}$ in formula (1) shows that the leading edge of the region of synchronous seizure activity migrates approximately 4–7 times more slowly than normal brain wave activity, which is in agreement with the experimental data described above.[8]

To summarize, mathematical modeling and theoretical analysis of large-scale electrophysiological activity provide tools for predicting the spread and migration of hypersynchronous brain activity, which can be useful for diagnostic evaluation and management of patients with epilepsy. It might be also useful for predicting migration and spread of electrical activity over large regions of the brain that occur during deep sleep (Delta wave), cognitive activity and in other functional settings.