Neurophysical Modeling of Brain Dynamics

Abstract

A recent neurophysical model of brain electrical activity is outlined and applied to EEG phenomena. It incorporates single-neuron physiology and the large-scale anatomy of corticocortical and corticothalamic pathways, including synaptic strengths, dendritic propagation, nonlinear firing responses, and axonal conduction. Small perturbations from steady states account for observed EEGs as functions of arousal. Evoked response potentials (ERPs), correlation, and coherence functions are also reproduced. Feedback via thalamic nuclei is critical in determining the forms of these quantities, the transition between sleep and waking, and stability against seizures. Many disorders correspond to significant changes in EEGs, which can potentially be quantified in terms of the underlying physiology using this theory. In the nonlinear regime, limit cycles are often seen, including a regime in which they have the characteristic petit mal 3 Hz spike-and-wave form.

INTRODUCTION

Correlations of EEGs with brain function are widely used diagnostically, and close connections to dynamics, cognition, and mental disorders are inferred (Niedermeyer and Lopes da Silva, 1999). Yet the detailed link between EEGs and the underlying physiology is not well understood, despite over 125 years' work (Niedermeyer and Lopes da Silva, 1999). Similar remarks apply to evoked response potentials (ERPs). Still more cryptic are EEGs seen in disorders, including epilepsy, whose relation to normal EEGs is not understood. As a result EEG studies are not integrated within any overall framework, nor with other branches of neuroscience.

EEGs result from cortical electrical activity aggregated over scales much larger than individual neurons or that can be modeled using neural networks. Hence, in one class of models averages are taken over microscopic neural structure to obtain continuum descriptions on scales of millimeters to the whole brain, incorporating realistic anatomy of separate excitatory and inhibitory neural populations (pyramidal cells and interneurons), nonlinear neural responses, multiscale interconnections, dendritic, cell-body and axonal dynamics, and corticothalamic feedback (Wilson and Cowan, 1973; Lopes da Silva et al, 1974; Nunez, 1974, 1995; Freeman, 1975; Steriade et al, 1990; Wright and Liley, 1996; Jirsa and Haken, 1996; Robinson et al, 1997, 2001, 2002; Rennie et al, 2002).

We have developed a physiologically based continuum model of corticothalamic dynamics that reproduces and unifies many features of EEGs, including the discrete spectral peaks in the slow wave, ‘delta’, ‘theta’, ‘alpha’, and ‘beta’ bands, seen in waking and sleeping states (Robinson et al, 2001, 2002), ERPs (Rennie et al, 2002), measures of coherence (Robinson, 2003), generalized epilepsies, EEG entrainment and seizure activation by stimuli (Robinson et al, 2002), and low-dimensional seizure dynamics (Robinson et al, 2002). Many behaviors were derived from moderate parameter changes of a few mechanisms in a single model, thereby enabling classification of different states using these parameters. Our approach averages over microstructure to yield mean-field equations in a way that complements cellular-level and neural-network analyses: these other approaches can be employed to elucidate the connections between microstructure and mean-field quantities, while the large-scale fields provide the background against which microscopic neural activity takes place. In the following sections we outline our model and its main results to date.

METHOD

Corticothalamic Model

The details of the model have been published elsewhere (Robinson et al, 2002); owing to space limitations, we restrict ourselves here to a brief outline.

The first feature incorporated is the neural response to the cell-body potential. Mean firing rates Q_a of excitatory (a=e) and inhibitory (a=i) neurons are nonlinearly related to mean potentials V_a by Q_a(r,t)=Σ[V_a(r,t)], where Σ is a sigmoidal function that increases from 0 to Q as V_a increases from −∞ to +∞. We use

where θ is the mean neural firing threshold and is its standard deviation.

The potential V_a results after dendritic inputs have been filtered and smeared out in time while passing through the dendritic tree, then summed. It obeys a differential equation (Robinson et al, 1997, 2001)

where β and α are the inverse rise and decay times of the cell-body potential produced by an impulse at a dendritic synapse and β≈4α. The right-hand side of (2) involves contributions from φ_ei, from other cortical neurons, and inputs φ_s from thalamic relay nuclei, delayed by a time t₀/2 required for signals to propagate from thalamus to cortex. In (2) ν_ab=N_abs_b, where N_ab is the mean number of synapses from neurons of type b=e, i, s to type a=e, i and s_b is the strength of the response to a unit signal from neurons of type b.

Each part of the corticothalamic system produces a field φ_a of pulses that travels at ν=10 m s⁻¹ and obeys a damped wave equation (Robinson et al, 1997, 2001):

where γ_a=ν/r_a and r_a is the mean range of axons a. If intracortical connectivities are proportional to the numbers of synapses involved, V_i=V_e and Q_i=Q_e (Wright and Liley, 1996; Robinson et al, 2001), which let us concentrate on excitatory quantities. The smallness of r_i also lets us set γ_i≈∞ (Robinson et al, 1997).

The model incorporates corticothalamic connectivities and thalamic nonlinearities. Figure 1 shows the connectivities considered, including the thalamic reticular nucleus that inhibits relay nuclei. The latter relay external stimuli φ_n to the cortex as well as corticothalamic feedback. The cell-body potentials then satisfy

where there is a delay t₀/2 for signals to travel from cortex to thalamus, c=r, s, ν_cc=ν_rn=0, and φ_c(r,t)=Σ[V_c(r,t)] (Robinson et al, 2001) applies because the small size of the thalamic nuclei enables us to assume γ_c≈∞ and r_c≈0 in (4).

Figure 1

Schematic of corticothalamic interactions, showing the locations ab at which ν_ab and G_ab act.

Our model has 15 parameters: Q, θ, σ′ α, β, γ_e, t₀, ν_ee, ν_ei, ν_es, ν_se, ν_sr, ν_snφ_n, ν_re, and ν_rs, enough to allow realistic representation of the anatomy and physiology, but few enough to yield useful interpretations. The parameters are approximately known from the experiment (Robinson et al, 2001, 2002), leading to the nominal values in Table 1, which are indicative only—some vary severalfold between individuals, arousal states, and disorders. We use only values compatible with physiology.

Table 1 Nominal Parameters from Robinson et al (2002)

RESULTS

Steady States, Linear Waves, and Stability

Setting all derivatives to zero in (3) and (4) yields steady states when the system is driven by a constant, uniform mean stimulus level φ_n. The equations are easily solved numerically, yielding a single low-φ_e solution, which corresponds to a normal state.

Small perturbations of steady states allow use of linear analysis. A stimulus φ_n(k,ω) of angular frequency ω and wave vector k has the transfer function to φ_e(k,ω):

where L = (1 − iω/α)⁻¹ (1 − iω/β)⁻¹ and $\overline{φ}$_a is the steady-state value of φ_a. This function is the cortical excitatory response per unit external stimulus, and encapsulates the relative phase via its complex value (Robinson et al, 2001; Rennie et al, 2002); it is the key to linear properties of the system. The gain G_ab is the differential output produced by neurons a per unit input from neurons b, and the static gains for loops in Figure 1 are S_d=G_esG_se for feedback via relay nuclei only, S_i=G_esG_srG_re for the loop through reticular and relay nuclei, and S_r=G_srG_rs for the intrathalamic loop.

Waves obey the dispersion relation q²(ω)+k²=0, with instability boundaries where this equation is satisfied for real ω (Robinson et al, 1997, 2001). In most circumstances, waves with k=0 (ie spatially uniform) are the most unstable (Robinson et al, 1997), and it is found that only the first few (ie lowest frequency) spectral resonances can become unstable. Analysis for realistic parameter ranges finds just four k=0 instabilities, leading to global nonlinear dynamics (Robinson et al, 2002): (a) Slow-wave instability (f≈0) that leads to a low frequency spike-wave limit cycle. (b) Theta instability, which saturates in a nonlinear limit cycle near 3 Hz (see Figure 2(a), where t₀=0.2 s for clarity, giving a frequency nearer 2 Hz (see below)), with a spike-wave form unless its parameters are close to the instability boundary. (c) Spindle instability at ω≈(αβ)^1/2 (see Figure 2, in the alpha band for physiological α and β, leading to a limit cycle near 10 Hz. (d) Alpha instability giving a limit cycle near 10 Hz, with a waveform like in Figure 2(b).

Figure 2

Sample time series from the model in regimes corresponding to (a) theta instability and (b) spindle instability.

The occurrence of only a few instabilities, at low frequencies, enables the state and physical stability of the brain to be parameterized in a 3D space with axes

which parameterize cortical, corticothalamic, and thalamic stability, respectively (Robinson et al, 2002). In terms of these quantities, the brain occupies a stability zone illustrated in Figure 3. The back is at x=0 and the base at z=0. A pure spindle instability occurs at z=1, which couples to the alpha instability on the upper boundaries, with spindle dominating at top and left, and alpha at right. At small z the left surface is defined by a theta instability. The front right surface corresponds to slow-wave instability and follows the plane x+y=1 to y=y_c≈−0.2. The boundaries correspond to onsets of generalized seizures (Robinson et al, 2002).

Figure 3

Stability zone for nominal parameters in Table 1, except α=60 s⁻¹. The surface is shaded according to instability, as labeled (dark gray=spindle, light gray at right=alpha, light gray at left=theta), with the front right face left transparent as it corresponds to a slow-wave instability. Approximate locations are shown of EO, EC, S1, S2, S4, REM (R), anesthesia (A), and alpha coma (C) states, petit mal onset (P), and the parameters in Table 1 (T), with each state located at the top of its bar, whose x−y coordinates can be read from the grid.

Spectra, Evoked Potentials, and Coherence

The EEG frequency spectrum is obtained by squaring the modulus of φ_e and integrating over

Figure 4 shows excellent agreement with observed spectra if is white noise in space and time, including occurence of alpha and beta rhythms at frequencies f≈1/t₀, 2/t₀, and the asymptotic low- and high-frequency behaviors; key differences between waking and sleep spectra can also be reproduced (see below and Robinson et al, 2001). The low-frequency 1/f behavior is a signature of edge-of-stability dynamics, which allow complex behavior (Robinson et al, 1997, 2001). To test our model and estimate some of its parameters, we fitted its linear spectrum to 103 normal adults' eyes-closed (EC) and eyes-open (EO) spectra (Robinson et al, 2002). This yielded mean parameters near those in Table 1.

Figure 4

Example spectra (solid) and fits (dashed) from a typical subject in EO and EC states.

A 1D wave-number spectrum results if one integrates |φ_e(k, ω)|² over one component of k:

We find that this spectrum is flat at small k_x, then approximates a power law, with P(k_x,ω)∼k_x^−g(ω). Figure 5 compares the exponent with EC data from Shaw (1991), showing excellent agreement for physiologically realistic parameters, except at f∼5 Hz, where the data are affected by phase distortions (Shaw 1991). In particular, spectral steepenings at resonances are reproduced.

Figure 5

Wave number spectral index g vs frequency. Model results are shown solid and data are dashed.

The inverse Fourier transform of (6) gives the ERP that results from an impulse stimulus. Initial work shows that the result agrees well with the experiment (Rennie et al, 2002), as illustrated by the example of a background (as opposed to target, in an oddball paradigm) ERP in Figure 6. Significantly, the model parameters used are essentially the same as those that reproduce the same subject's prestimulus EEG spectrum.

Figure 6

Experimental (solid) and model (dashed) ERPs in response to a stimulus.

The cross spectrum P(r, r′, ω) is the phase average of φ_e (r, ω)φ_e(r′, ω). The coherence function is then

Figure 7 shows that this gives good agreement with observations for model parameters close to those used in obtaining the other plots in this work (Robinson, 2002). A rise seen in γ² at large R=|r−r′| (Nunez et al, 1999) is also reproduced (Robinson, 2002).

Figure 7

Model alpha-band coherence vs distance (solid) for EC conditions, compared with experimental curves.

States of Arousal

Nonseizure states lie within the stability zone in Figure 3. Detailed arguments regarding the arousal sequence, from alert to deep sleep, and including REM sleep and sleep stages 1–4 (S1–S4), constrain the relevant regions of parameter space and place this sequence as shown in Figure 3 (Robinson et al, 2002). States such as anesthesia can also be represented.

Figure 8 shows model time series for parameters illustrating EO, EC, S2, and S4 states, holding Q, γ_e, t₀, β/α, ν_ei, and ν_sn at their nominal values, and varying α and the other ν_ab only moderately. The features seen strongly resemble those of corresponding experimental data (Niedermeyer and Lopes da Silva, 1999). Known differences between EEG spectra for subjects with differing disorders enable classification of these conditions into different parts of the stability zone, while seizures correspond to departures from this zone, as discussed next.

Figure 8

Model time series at labeled points in Figure 3: (a) EO, (b) EC, (c) S2, and (d) S4.

Petit Mal Seizures

Petit mal is a common generalized epilepsy. Seizures are mostly seen at ages 4−20, last 5−20 s, cause loss of consciousness, and show a spike-wave cycle which starts and stops abruptly across the whole scalp (Niedermeyer and Lopes da Silva, 1999). The frequency falls from around 4 Hz to under 3 Hz in most cases, and non-REM sleep is a powerful seizure activator. Experiments show that the loops in Figure 1 are essential for petit mal, with the cortex synchronizing thalamic activity (Niedermeyer and Lopes da Silva, 1999; Steriade et al, 1990). GABA antagonists such as penicillin can start spike-wave oscillations, in some cases converting spindles to spike-wave complexes, similar to those also seen in some partial seizures (Niedermeyer and Lopes da Silva, 1999).

Our model gives an ≈3 Hz spike-wave cycle as the nonlinear stage of theta instability (Figure 2) and we conclude that this corresponds to petit mal (Robinson et al, 2002). Analysis shows that this cycle consists of a flip-flop (it alternates between two states) in the limit γ_e, α→∞, a residue of which is seen in Figure 2(a). The high-φ_e part corresponds to large φ_s incident on the cortex as a result of low φ_r a time t₀/2 earlier and low φ_e a time t₀ earlier; the low-φ_e part corresponds to the converse, with near silence in relay nuclei, as seen experimentally. Signals make two circuits of the system before it returns to its original state, giving a period 2t₀. At finite α, signals traveling via the reticular nucleus are delayed by ≈1/α more than those that only pass through relay nuclei. Hence, when φ_e flips to its upper state, there is a short period t₀/2 later when S_d>|S_i|, resulting another t₀/2 later in a spike of duration ≈α⁻¹. Finite α and γ_e also round off the other side of each square wave and finite α leads to damped spindle oscillations at ω=(αβ)^1/2. These mechanisms accord with the experimental inferences above, and observed EEGs often show a residual flip-flop plateau in each cycle (Niedermeyer and Lopes da Silva, 1999; Robinson et al, 2002).

Estimation of the petit mal period gives T≈2t₀+6/α+6/β+4/γ_e (Robinson et al, 2002), consistent with observations and the insensitivity of T to most parameters. The main features of the waveform, apart from spindles, are found in cases with β, γ_e=∞ but finite α, which implies that the 3D system resulting in that limit contains the essential dynamics. This accords with findings that dimensions of time series of petit mal and related seizures are low (Babloyantz and Destexhe, 1986).

A typical onset point for a petit mal seizure is shown in Figure 3. Transformation of spindles into petit mal is inferred to occur by moving from the vicinity of S2 to the theta instability zone, with a rapid switch of activity from roughly 10 Hz to 3 Hz. Large values of γ_e favor instability, which may explain the onset of petit mal at around age 4, since γ_e rises in children because of myelination.

DISCUSSION

We have developed a model that incorporates the main features of corticothalamic physiology and anatomy using only 15 parameters. Its predictions provide a unified description of a wide range of phenomena, with six parameters fixed across all states, and the others only varying moderately. Of key importance is the xyz parameter space in which the stability zone of the brain is easily visualized, and in which disorders, states of arousal, etc, can be classified. Within this zone, linear analysis gives good approximations to EEG spectra, ERPs, coherence and correlation functions, and related measures. Zone boundaries are identified with onsets of generalized seizures, consistent with known features of their time series and patterns of occurrence. The model also explains low-dimensional dynamics in petit mal seizures and other nonlinear behaviors such as alpha entrainment and seizure activation are also reproduced (Robinson et al, 2002). A key feature of our approach is that we extract a broad range of behavior from modest changes in the parameters of a single model, without postulating extra mechanisms.

Fitting the model's predictions to data provides a noninvasive probe of large-scale physiology that yields parameter values consistent with independent measures. This enables states of arousal, seizure onsets, and pathologies to be assigned to distinct regions of parameter space. We have found that the normal arousal sequence has a simple form in the xyz space, that clinically observed waking states lie at inferred locations, and that seizure onsets lie close to the most commonly seen precursor states (see also Robinson et al, 2002). This space thus provides a physiologically based organizing framework for a wide variety of phenomena. Its topography may indicate new connections among phenomena in neighboring regions, and enable the significance of the parameters that distinguish the various cases, and cause transitions between them, to be studied systematically.

The approach discussed here provides a powerful framework for further studies: It remains to investigate what factors control progression along inferred arousal sequences, or onset of seizures, for example. This will require inclusion of neuromodulator dynamics and brainstem reticular activation, possibly using additional feedback loops. Spatial variations are also under investigation, with successful initial application to understanding split alpha peaks seen in a few percent of the population (Robinson et al, 2003).