Walter J. Freeman Journal Article e–Reprint

Methods of Analysis of Brain Electrical & Magnetic Signals:

Analytic Techniques Used in the Search for the Physiological Basis of the EEG

WALTER J. FREEMAN

Department of Molecular & Cell Biology, University of California, Berkeley, CA 94720, USA

Reprinted from Chapt 18 of EEG Handbook (revised series, Vol. 1), A.S. Gevins and A. Rémond (Eds.) 9 1987 Elsevier Science Publishers B.V. (Biomedical Division)

In the course of generating behavior, the brain also generates time-varying electric and magnetic fields. The aim of electroencephalography is to record and measure samples of these fields during certain states and sequences of behavior, in order to explain some of the mechanisms by which behavior is generated.

In this chapter the brain is treated as a hierarchically organized network of neurons whose dynamics determine behavior and vice versa. The dynamics is largely determined by the strengths of interaction among the neurons. The tasks of physiologists are:

(a) to measure accurately the time courses and spatial patterns of potential fields,

(b) to simulate the wave shapes of these events with the solutions to differential equations or their equivalent that represent the dynamics of neural networks,

(c) to estimate strengths of action ('gains') from rates of change of potential, and

(d) to interpret changes in strengths of neural action in relation to changes in behavioral states.

Procedures are outlined for acquiring and measuring suitable data bases, and for constructing and evaluating explanatory equations. Although the brain is a distributed and highly non-linear system, the treatment begins with the development of lumped, linear systems. This is because the construction of a neural dynamic model requires specification of the network of elements and their connections, definition of state variables and observables, types and locations of non-linear operators, measurement of open loop time and space constants, and evaluation of state-dependent coefficients. Analysis of non-linear properties presupposes an understanding in detail of these and other features of a system. Fortunately, linear domains of complex systems are ubiquitous in the brain when they are sought; they are the best preparatory path of access to the distributed, non-linear, multiply looped feedback systems that generate behavior and the EEG. Emphasis is given to techniques for fine-grain spatial analysis of the EEG, for use of behavioral correlation to validate and optimize procedures of measurement, and for empirical detection of chaos, limit cycles, and random components of EEG recordings. Examples are drawn from the mammalian olfactory bulb, in which the neural non-linear dynamics manifested in the EEG sustain the basic operations in olfactory perception.

18.1. The need for models

The neuron releases electromotive forces (EMFs) at synapses and trigger zones that drive electric current along dendrites and axons, so that activity is transmitted from one part of the neuron to another. The current must flow in a closed path; it crosses the membrane in one direction at the site of its EMF, crosses in the other direction elsewhere including its sites of action, and passes through the extraneuronal tissue space in its return path. As the current flows across the resistance of the tissue, it causes a low amplitude field of potential. That field constitutes a sign of the activity, but this sign is not itself the activity. The activity is the space-time pattern of depolarization and hyperpolarization of the neuron that is established by the flow of current across the membrane resistance, it is not the potential difference established by the flow across the extraneuronal tissue resistance. The manifestation is like the sound of a single part of an engine; it mixes with the sounds of all other parts to create a roar. Our task of analysis is to localize and measure the pings and rumbles of particular parts, so as to infer what the parts contribute to the function of the whole engine. While it is conceivable that extracellular EEG currents might contribute directly to brain function, just as sound vibrations might partly determine the performance of an engine, that possibility is disallowed here.

The key to understanding brain function lies in determining the strength of action (represented by a coefficient called the 'gain') of one group of neurons on another, and the changes in gain with altered behavior. When we can isolate two groups and show that the action is unidirectional, then the problem is mainly to demonstrate for each group the relation between its activity and its manifestation of activity (its sign to the observer, not to the other group). The gain is given by the ratio of the amplitudes of the two activities.

Groups in the intact brain are not isolated. They function in networks with parallel forward and feedback loops. The problem of determining activities from the signs of activities is the same as in the isolated case, but the procedures for estimating the forward and feedback gains are more complex than taking ratios of amplitudes.

The basic procedures in dynamics are to measure the rates of change and the frequencies in signals, and to construct a model so as to infer strengths of action as gains from the measurements. The model typically is a differential equation that allows us to infer strengths of synaptic interaction among groups of neurons from recordings of evoked potentials and the EEG.

As an assembly of neurons the brain is structured by their interconnections into a hierarchy of parts. Locally dense interconnections, synaptic and non-synaptic, form local modules such as columns in cortex and clusters in nuclei of the basal ganglia. Axonal tracts, usually in bidirectional pairs, serve to interconnect these local parts into subsystems comprising areas of cortex and their related parts of the basal ganglia. The subsystems are integrated to form major entities such as the motor, visual and limbic systems; the interactive total comprises the functioning brain.

In this review the description of these parts for computational purposes is by a hierarchy of sets (Freeman 1975, pp. 25-34). A KO set represents a spatially distributed collection (KO mass) of neurons having a common input (e.g., a set of sensory receptors, or the first order neurons to which they project) and a common sign of output (excitatory, +, or inhibitory, -) but no interaction between them. A KI set models a collection (KI mass) of neurons that has the first two properties but that is densely interactive by interconnections, either mutually excitatory (KI_e) or mutually inhibitory (KI_i). A KII set describes the interaction of a collection (KII mass) of neurons consisting of a Kl_e mass and a KII mass. KI and KII sets represent the neuronal assemblies commonly comprised by areas of cortex of nuclei; KI and KII subsets suffice for local collections such as cortical columns and nuclear modules. A network comprising a subsystem is modeled by the conjoining of KI and KII sets into a KIII set. A network of KIII sets may serve to represent a neural system at the KIV level, although this step remains to be demonstrated.

The dynamic action of neurons in the brain is regulated by the input from sensory receptors (interoceptor and exteroceptor) and by the strengths of the interconnections determined by genes, by past experience incorporated by learning, and by the chemical environment in the brain. The play of the active states is observed in goal-oriented behavior and in the extraneuronal electric and magnetic fields generated by the neurons. The aim of electroencephalography is to describe regular relationships between these two observables, such that underlying neural active states can be discerned and the dynamics of the generators can be understood.

To achieve this aim 3 levels of inference using models are required. First, a volume conductor model is used to establish relations between neural active states and the fields of potential that signify them. Given an electric or magnetic event and a behavioral event, the neurons generating the electric currents must be identified. Generally the neural elements for behavioral correlation are KO and KI masses rather than single neurons. Estimates of their activity are derived by appropriate ensemble averaging over the activities of neurons comprising them. By use of neuroanatomy, neurophysiology, neuropharmacology, and volume conductor theory an electrostatic or electric current model is devised for each postulated generator of a component in the event, and the spatial distribution of that component is tested for compatibility with the prediction of the model. Upon completion of the analysis we should be able to infer that a pattern of potential difference of a certain kind over points near or in the brain manifests excitation or inhibition of a certain KO or KI mass or a portion thereof.

Second, a model is required to relate given neural active states to neural operations. When a neural mass is given a certain space-time pattern of input, it generates another pattern of output. In effect it operates on the input to change it, add to it, replace it or delete it. The nature of the operation is inferred by comparing the input and output active states. The preferred language of description is by an integrodifferential equation, in which an operation is defined by the ratio of output to input. By constructing and solving a set of equations for K sets representing the dynamics of a neural mass, and by further transforming the solution in accordance with the volume conductor model, we can expect to derive curves that can be closely fitted to observed neuropotential patterns.

Third, a model is required to relate neural operations to behavioral operations. This is the turn of a coin by which a neural event is identified with a psychological event. This, the most difficult of the 3 steps, is embodied in the experimental paradigm in which the electric fields are measured. Unless behavior is properly structured and controlled, the most precise measurements of potential are likely to be meaningless.

This chapter is mainly about the second stage, in which data are collected with properly positioned electrodes in an appropriate behavioral context, and in which the task, is to construct a set of differential or difference equations to describe the space-time dynamics of neural systems. The choice of equations is a matter of convenience, but the physiologist must choose on the basis of prior knowledge of the cellular anatomy and physiology of the system of interest, including data from unit and evoked potential studies that are independent of EEG recordings. In order to lay out the underlying principles of analysis, this discussion begins with lumped systems and ordinary differential equations. The extension to distributed non-linear systems is straightforward for those who have or will acquire the necessary mathematical tools. For the same reason the examples begin with ERPs and poststimulus time histograms (PSTHs) from electrical stimulation of brain pathways; the approach is the same for treatment of events related to sensory stimulation and the EEG.

Readers who are interested in mathematical and systems analysis of large scale neural activity and the EEG will find related and alternative approaches in the references at the end of this chapter. Particularly noteworthy are books by Griffith (1971), Simon (1972), Lin and Segal (1974), Basar (1976), MacGregor and Lewis (1977), Thatcher and John (1977), Childers (1978), and Nuñez (1981), and reports by Wilson and Cowan (1972, 1973), Lopes da Silva et al. (1976), and Ingber (1982).

18.2. Establishing the data base

The study of neural dynamics is based on input-output pairs: a stimulus is measured and given, a response is recorded and digitized, and the numbers representing input, output and context are stored. The procedure is repeated 10²-10⁵ times until there are no further surprises. The pairs are grouped in accordance with input and context parameters; they might also be grouped tentatively in relation to output, if there are salient features such as the apparent presence or absence of major components or distinctive wave shapes with no apparent intermediates to suggest gradation.

Four aspects of the nervous system must be kept constantly in mind during these explorations: background activity (often mislabeled 'noise'), instability, time variance, and non-linearity. First, all normal nervous systems are continually active whether or not stimuli are given. The data base should include a null class with no known stimuli as a basis for characterizing intrinsic activity. The use of randomized stimuli with respect to time windows of observation is unsatisfactory, because the distributions of amplitudes may be warped by evoked potentials. These records should be processed in the same way as all others, including ensemble averaging, spectral analysis, amplitude histograms, calculation of moments, covariance matrices, etc.

Second, in collecting pairs the assumptions are made that the neural system is in a certain state prior to stimulation, that it is perturbed by the stimulus away from that state, and that it returns to the same state. More specifically it is essential to include in each output member a segment of prestimulus baseline, in order to verify that the presumed state does hold, and that the preceding response has terminated. More generally we should remain aware that the stimulus may elicit an instability, such that the neural system may change from a basal state to another state and perhaps back again some time later. Some well-known examples include seizures, spreading depression, kindling, long-term potentiation, and learning, but subtler shifts may hold the key to analysis. The point is important because the solutions to the neural dynamics in each state can be described by a set of time-dependent equations, and a new set may be required for each new state. This is most likely to be the case if each state is described by a linear approximation.

Third, it is a truism that the nervous system is never twice the same, and that even in a certain state, defined by both behavioral and neural variables, there is time variance such as by fluctuation, drift or learning. These variations are entered into a set of differential equations by use of coefficients that vary randomly or with time or with contextual parameters such as satiety or adaptation, on the assumption that the time scale of this kind of variance is very long or very short in comparison to the duration of an output event. The coefficients are fixed for the duration of the simulation of each event, and they are changed to new values with each simulation, if necessary.

Fourth, linearity denotes the property of a limited domain of dynamics in which superposition holds. In this kind of state outputs from multiple inputs are additive, and their magnitudes are proportional to the input magnitudes. Experimentally this domain is demonstrated by paired-shock testing, in which two stimuli are given first separately and then together. The algebraic sum of the separate digitized responses equals the compound digitized response if the neural system is perturbed within one of its linear domains. If the sum of single responses differs from the compound response, the system is in a non-linear domain.

Neurons and neural systems have multiple linear and non-linear domains. The demarcation requires extensive parametric testing using systematic variation of stimulus intensity, duration, rate and interval. A non-linear domain may be linearized by re-definition of variables, or it may be divided into smaller domains that are described by piece-wise linear approximations. For example, the action potential manifests a non-linear dynamic process, but trains of action potentials can be described by pulse frequency or density, which in itself may be a linear variable with respect to transmembrane potential.

However tedious parametric analysis might be, it is essential as the basis for constructing a first dynamic model of a neural system. We will do well to recall that 20 years of subthreshold testing of axons in their linear domain was a necessary basis for the construction of the explicitly non-linear Hodgkin-Huxley equations for the nerve action potential. The dynamics of a system in one of its linear domains is described by a set of linear differential equations. Once we establish linear domains we can gently push the system outside and begin to ask questions of the location and nature of non-linearities, whether static (as in much evoked potential research) or dynamic (as with the action potential), and of the nature of the thresholds and stability properties that characterize one or more non-linear domains. Indeed, it is in these non-linear domains where answers to questions on the neural bases of behavior are to be found. The doorway to access is provided by the linear differential equation or its equivalent.

18.3. The use of basis functions

The search for linear domains requires use of low intensity stimuli in order to keep the neural system in a low amplitude range. This does not guarantee linearity, but conversely it is nearly certain that if the evoked potential amplitude exceeds 2 or 3 standard deviations of background activity, the neural mechanism will be driven into a non-linear domain. Owing to the background activity, the low-amplitude evoked potential can be observed and measured only with the aid of averaging over ensembles to give an event related potential (ERP). The variability about the average decreases with the square root of the number of events in the average (see Ch. 5). However, to the extent that the evoked potential varies over the time required to collect the ensemble, the information in the ERP undergoes destructive averaging; the true shape may be lost, together with information about the manner of variation.

Description of the particular wave shapes of ERPs is often considered of secondary importance, the principal aim being to determine the variation of shape within and between sets of ERPs in relation to associated patterns of behavior. However, the quantitative description of that variability by use of the variance in sets of numbers requires that the shapes be clearly characterized and precisely measured. Further, the information cannot be used to construct explanations of neural dynamics until the wave shapes are fitted with curves.

Most ERPs have irregular shapes that are difficult to describe and comprehend, and their shapes vary erratically and often hugely from trial to trial. The two main sources of the variation are the background EEG activity and the changes from trial to trial in the neural mechanisms generating the ERPs. When the experimental conditions relating to behavioral stimulus and response are fixed to the extent possible, the variation in wave-form is reflected in the numbers by what we will call 'uncorrelated variance.' When an experimental condition is deliberately changed across trials, the ERPs can be expected to change as the neural mechanisms change, and this introduces 'correlated variance.' A good system of measurement should produce matrices of numbers that preserve the correlated variance and minimize the uncorrelated variance. The aim is not to suppress the variance but to partition it most effectively.

The primary step in measuring the form and variance of an ERP is to decompose it into components. Each has 3 aspects. First, each component is an elementary curve. These elementary curves are called 'basis functions.' When a selected set of such curves is added together, the sum should closely conform to the shape of the ERP. Second, each component has a set of numbers or coefficients that denote its height, width, rate of change, and other features. Third, each component defines one or more coordinates, and the set of numbers denotes distances along the axes. The number of coefficients in the set of basis functions specifies the dimensions of the measurement space. In this point of view the measurement of each ERP gives a vector in that space.

The basis functions must be chosen so that each one represents unique aspects of the ERP. This implies that the basis functions are independent of each other, and the coordinate axis of each of the coefficients is not parallel to any other axis, although it need not be orthogonal to any other. If they are not orthogonal, then pairs of coefficients have non-zero correlation across sets of ERPs. For example, the phase and frequency of a damped cosine function fitted to ERPs will vary inversely, if the frequency varies from one ERP to the next, but the duration of the first peak of the ERP will be fixed. In general, the shapes of ERPs tend to vary in certain ways and not others owing to constraints in the underlying neural mechanisms. These constraints become evident on correlation and factor analysis of sets of coefficients (Freeman 1975, pp. 408-414). They are exceedingly valuable in model-building.

Each behavioral measure is also represented by an axis. The variation in wave-forms of the ERPs and in behavior is represented by the variance of the numbers on each axis. If there is covariance between the ERPs and behavior, the axes of the ERP measurement space should be positioned (by the proper choice of the basis functions and by varimax rotation in factor analysis - see Ch. 5 and 16) to project strongly onto the behavioral axes in order to display that covariance.

From an alternative perspective the ERP is conceived as a signal embedded in noise. Here the term 'noise' includes background EEG activity that is uncorrelated with the imposed stimulus train as well as artifacts and electrode thermal noise. The signal is asserted to have certain measurable attributes such as amplitude, latency, frequency, etc. and the aim is to measure these in the presence of noise. With the use of a digital computer we build a template or filter that is expected to match the signal in order to identify it in the noise. The form of the signal is expected to vary with changes in the behavioral conditions. That is, the filter must be adaptive to the varying form of the signal. Furthermore, the adaptive filter consists of the sum of a set of elementary curves (basis functions) that is fitted to the ERP. Hence in both views the measurement process consists in curve-fitting, and the requisite step is to choose the basis functions, which may represent the parts of the ERP.

18.4. Decomposition of ERPs with analytic basis functions from the AERP

An ERP is commonly interpreted as consisting of a sequence of peaks labeled N1, P1, N2, etc., each of which is measured by height and latency. This basis function is an empirical upright or inverted 'U.' Although useful in preliminary work this approach has several drawbacks. Much of the information in the ERP is lost to the extent that the amplitude and latency numbers extracted fail to allow reconstruction of the ERP. The numbers give an inadequate basis on which to evaluate the rates of change in the ERP, which are more valuable than the amplitudes of the peaks. This is because interactive brain systems are more likely to reveal changes in their sensitivities by changes in the rates and durations of their activity patterns than by fluctuations in the amplitudes of activity. Uncertainties arise when adjacent peaks merge or overlap. There is no way to partition the uncorrelated variance from the variance that is related to stimulus and response variables.

The deficiencies of this approach can be remedied by using sets of elementary curves that are called 'analytic basis functions.' Each is a continuous curve that extends over the entire duration of each ERP and is defined by an algebraic term or expression containing one or more coefficients.

The procedure for using them has two stages. In the first stage an ensemble average ERP (AERP) is made over a delimited set of ERPs. The individual ERP contains both signal and noise; grand averaging minimizes the noise, although it does not entirely eliminate it. Then a curve consisting of the sum of a set of analytic basis functions is fitted to the AERP by adjusting the values of all the coefficients. This curve is called a 'matched filter.' The fitted curve is treated as an approximation for the signal, and the difference between the AERP and the curve is treated as residual noise.

We assume that the noise differs on successive trials, and that the shape of the signal also differs due to the correlated variance. The shape of the average for the signal must be unlike any of the signals; it is only an estimate for the shape of each one. For example, if we take the ensemble average of a set of exponential curves with a rate coefficient that varies randomly on successive trials, the mean curve is not exponential but can be fitted with one to an approximation. In addition to the distortion of the signal, the AERP of course contains none of the variance of the ERPs.

In the second stage the matched filter is used to measure the signal embedded in noise in each of the ERPs. A curve from the same equation is fitted to each of the ERPs in the set by readjusting the coefficients. That is, the filter is adapted to each ERP. The variance in the set of ERPs is incorporated into a matrix of coefficients (usually 4-10) of the basis functions, and the distortions from the exponential form introduced into the AERP may be largely corrected.

The standard deviation (S.D.) of the AERP may serve to evaluate the proportions of correlated and uncorrelated variance in a preliminary assay. If the variance is determined by the noise and not by changes in the shape of the ERPs, then the S.D. is relatively constant over the time of the ensemble average. If there are changes in the ERPs, then the S.D. is not uniform. An example is shown in Fig. 1 of the AERP over ERPs from 20 trials at each of 3 stages in the training of a cat to respond to an electrical shock to the lateral olfactory tract (LOT) as a conditioned stimulus, while the evoked potential is recorded from the olfactory cortex. Early in training with no significant behavioral responses the S.D. is uniform. During acquisition the wave shape changes, and the S.D. increases. In the later stage of consolidation the wave shape reverts near to its original form, but the S.D. remains high in the neighborhood of the high amplitude of the AERP, reflecting the variation in shape of the component ERPs. Higher moments (the cubic and quartic) serve mainly to detect outliers and artifacts, which can be deleted to improve the quality of the AERP.

Fig. 1. The relations are shown between the standard deviation (solid curves) and mean (dotted curves) of the AERPs of 3 sets of 9 event related potentials (ERPs) from implanted electrodes in a cat being trained to press a bar to the electrical stimulus (olfactory tract) used to evoke the potentials (prepyriform cortex). Above: response on <10% of trials; middle: >50% below: >90% (Freeman 1975).

This procedure overcomes the deficiencies in raw digitized data. The number of axes in the measurement space is reduced by an order of magnitude. Each coefficient is evaluated on the basis of what is happening over the entire duration of an ERP, so that information on both amplitudes and their rates of change is preserved. Most importantly, curve-fitting introduces a structured hypothesis into the measurement process. It asserts on the basis of the grand average that particular parts or components such as exponentials or damped sine waves may be used to describe the ERPs. When bursts of background EEG waves are superimposed on the potential related to a behavioral event, and the signal in the ERP is obscured by noise, the basis functions impose a form that is adaptive but cannot be twisted arbitrarily.

More explicitly, if a neural mechanism has certain stable properties, such that the generic character of its neural dynamics manifested in a set of ERPs is invariant, then those ERPs from single trials can be approximated by a set of curves all from the same equation over sets of trials. The shape of ERPs varies from condition to condition due to variations in the mechanism, and these variations are expressed in a matrix of the coefficients in the equations of the fitted curves. Also each ERP contains noise contributed by the background EEG. The noise can be expected to distort the values of the coefficients from their putative true values (in the absence of noise) on each trial.

However, a difference will persist between the fitted curve and the ERP. The residuals may be rejected after inspection, because the set of basis functions is not established in respect to individual ERPs but to the AERP. This rejection is the main basis for partitioning the variance, with retention only of that part that is related to the coordinate space for measurement, i.e., the set of basis functions, because the AERP serves to impose a structured hypothesis onto the data. The matrix of coefficients will contain less information than the digitized raw data, but that content should be more richly related to the neural events underlying the ERPs and therefore to behavior. However, this hypothesis should not go untested; we should always inspect the residuals and test them for correlation with behavior before discarding them.

These procedures are applied to the EEG in Section 18.10, where an ensemble average of a spatially distributed set of EEG segments simultaneously recorded from an electrode array (an AEEG) plays the role assigned here to the AERP.

18.5. Spectral analysis of ERPs

At this stage we may use the Fourier transform and spectral analysis to decompose the wave shape of the ERPs into sets of sines and cosines over a range of frequencies (see Ch. 4). There are several reasons for this. First, the ERPs may contain oscillatory components that will require use of oscillatory basis functions. Second, the results may suggest improvements in recording procedures to enhance or diminish one or more parts of the spectrum by pre-filtering the ERPs; or, if the main frequencies are separated by at least a factor of two, the ERPs may be decomposed into two or more parts by pre-filtering. Third, spectral analysis may be applied to the background EEG, which is a prime source of the variance of the ERP. If there is a strong periodicity in the EEG, the stimulus repetition rate should be set at an odd multiple of the half-length of the period. For example, if a set of cortical evoked potentials is recorded in the presence of strong gamma activity at 40 Hz, the stimulus rate should not be 5/sec (which might lock on to a burst of gamma at any phase) but might be 4.5 or 5.5/sec; or the interstimulus interval might be randomized.

Fig. 2. The characteristic effect on ERPs of increasing the stimulus intensity is a spectral shift to the left. The lowest curve is the amplitude spectrum of the prepyriform EEG in a cat. The remaining curves in ascending order are the amplitude spectra of the Fourier transforms of ERPs (normalized to equal maximal amplitude) over a range of increasing stimulus intensity. The asterisk shows the level of the threshold for behavioral responses to the electrical stimulation, at which the single evoked potentials are close in amplitude to the background EEG. At this level the ERP spectrum most closely conforms to the EEG spectrum. (From Freeman 1975.)

Yet another use of the EEG spectrum is to compare it with the spectrum of the ERP (Fig. 2). Do the spectral peaks occur at the same frequencies in the EEG and the ERP, and is it desirable or not that they do so? From the standpoint of observability of the single ERP it is not desirable, because the EEG and ERP cannot be separated by bandpass filters. Moreover, such correspondence of spectral peaks happens most often when the amplitude of the single-trial ERP is the same or less than the amplitude of the EEG. Because the only distinguishing characteristic of the ERP is the time of its occurrence, averaging in this case is a necessity.

From an alternative standpoint, the stimulus is seen as superimposed on background input to a functional group of neurons, and the characteristics of the response to the stimulus should closely resemble those of the response to the background input. With a strong stimulus the neural system may respond abnormally, usually involving saturation or partial blockage of the neural mechanism, so that the normal range of ERP variance is diminished. Then strong stimulation, which is commonly used to improve signal to, noise ratios, should be avoided. After all, the aim of analysis is not the suppression of variance but decomposition and selective removal. Parametric testing should be used, in which the stimulus intensity or some other experimental variable is changed in small steps over a suitable range. For example, a series of ERPs and standard deviations can be obtained over a range of stimulus intensities or output amplitudes (as appropriate), the dependence of wave shape on intensity noted, and a range with optimal variance chosen for further study such that there is an optimal covariance with behavior (Fig. 3). It may also prove useful to adjust the stimulus intensity to maximize the non-uniformity of the second moment with respect to ERP time (Fig. 1).

Fig. 3. Two sets of ERPs (dotted curves) are shown from the olfactory bulb on stimulation of (right) the lateral olfactory tract (LOT) and (left) the primary olfactory nerve (PON). From above to below the stimulus intensity is raised in steps, which changes the amplitude and shape of the ERPs. Solid curves: sum of fitted basis functions. Dashed curves: dominant and subsidiary components (respectively ‘signal' and 'clutter' - see Section 18.6). (From Freeman 1975.)

Typically we have some control over the correspondence between the spectral peaks of the ERP and EEG because we can change the stimulus rate or intensity which will change the spectrum of the ERP. The most likely result of increasing the stimulus intensity is an increase in amplitude of the ERP above the EEG and a decrease in the frequencies of the ERP (Fig. 3), which is reflected in a prominent spectral shift (Fig. 2). This makes an ERP much easier to retrieve from the background activity, but it also decreases its variance. This may be an instance of throwing out the baby with the bath water, because in attempting to observe the signal, we may excessively reduce the correlated variance as well.

The Fourier transform of an ERP or AERP also provides a basis for estimating the amount of information that can be extracted. The phase spectrum should be displayed along with the gain spectrum preferably in the form of a polar plot (Fig. 4). Each peak or trough in the gain spectrum may represent a characteristic frequency (called a 'pole') or a null (called a 'zero') in the system (Appendix B). We can deduce the order of the system from the goodness of fit of curves to the data in both the time and frequency domains, i.e., the number of simultaneous first order differential equations to be evaluated. Each 90° of phase lead or lag will require a first order differential equation in the set needed to describe the dynamics, or equivalently yet another term in an n-th order equation. In the example in Fig. 4 the polar plot can be viewed as the track of the tip of a vector that rotates clockwise around its base at the origin from low to high frequency. It sweeps through 6 quadrants comprising 540°; the neural system is minimally sixth order. The goal of curve-fitting at this stage is not parsimony (the smallest number) but evaluation of as many rate coefficients as are indicated by the system. If the neural system is operating in a linear domain, each term or equation will represent a real dynamic transformation. If it is not, then some of the terms or equations will reflect the dynamics spuriously. For example, a saturation non-linearity may give rise to harmonics.

Fig. 4. Above left: ERP from a cat with implanted electrodes (dotted curve) fitted with sum of basis functions (solid curve). Above right: amplitude and phase of the two curves at left after Fourier transform by numerical integration. Below left: polar plot of the Fourier transform. Lower right: dominant and subsidiary components, which are subsets of basis functions comprising the curve fitted to the ERP. (From Freeman 1975.)

Linearization is principally desired by the fact that the curves comprising the solutions to linear differential equations include pulses, steps, ramps, exponentials, sines, cosines, damped sine waves, and several other elementary wave-forms. This may be called the 'family of linear basis functions.' There are numerous electrophysiological systems at both the single neuron level and the level of masses of neurons for which the dynamics in certain ranges can be described with linear differential equations. On driving by electrical impulses each output must then consist of the sum of a set of linear basis functions. The Fourier transform is an example of linear decomposition of a signal into sinusoidal basis functions; it is most effective when the system that generates the signal is operating in a linear range.

At this stage we return to the experiment and again carry out careful parametric testing (which presumably has already been done), i.e., systematic variation of stimulus interval and intensity (as appropriate) to determine where linear ranges of function exist (Freeman 1975, Ch. 2) in which additivity and proportionality hold. An example is shown in Fig. 5 of the ERPs from single-shock responses recorded in the hippocampus on electrical stimulation of the fornix (Horowitz 1969). Over a range of low stimulus intensities (a-c) the ERP increases in amplitude without change in its form (proportionality). Over a middle range (d-h) the generic character of the ERP is invariant (approximating a damped sine wave), but the frequency and decay rate of the oscillation change as the initial amplitude increases. At the highest stimulus level shown (i) the ERP changes to a new generic character. The 3 ranges are linear, piece-wise linear, and non-linear.

Fig. 5. ERPs are shown from the ventral hippocampus at increasing amplitudes of single shocks given to the dorsal fornix. (a) 2.5 V, N = 6600; (b) 3.0 V, N = 3600; (c) 4.0 V, N = 1800; (d) 5.2 V, N = 850; (e) 7.0 V, N = 650; (f) 9.0 V, N = 510; (g) 11.0 V, N = 417; (h) 15.0 V, N = 357; (i) 20.0 V, N = 290. Duration of each ERP is 125 msec. (From Horowitz 1972.)

Alternatively we may systematically explore the background conditions for recording. An illustration of this procedure is shown in Fig. 6. The set of ERPs is taken from a cat with stimulating electrodes in the optic tract and recording electrodes in the superior colliculus. The cat is placed in a light-tight box and the uppermost ERP is taken. Thereafter the ERP is taken at the designated times during presumptive dark adaptation. Control sequences are taken in the same box with illumination to control for the effects on the ERP of habituation, changes in motivation, etc. (Pickering and Freeman 1967). The sets of ERPs from the hippocampus and superior colliculus resemble damped sine waves. The degree of conformance is tested by fitting to each one a curve generated from the equation for a damped sine wave. The fit is rather close over the linear (Fig. 5, a-c) and near-linear (Fig. 5, d-h) ranges for the hippocampal ERPs (for example, lower frame in Fig. 7) but not for the non-linear range (Fig. 5, i). The collicular responses conform only when the initial positive peak (Fig. 6) is excluded from the curve-fitting procedure (Fig. 7, upper frame). In both sets the frequency, decay rate and phase of the damped sine wave change in relation to the antecedent variable. Readers may profitably compare these illustrations of ERPs from the olfactory bulb and cortex, the hippocampus, and the superior colliculus with examples commonly seen in the literature and with the common varieties of human ERPs (Ch. 5) in order to see how much can be learned by the systematic search for linear domains.

Fig. 6. A set of ERPs is shown from a cat with electrodes chronically implanted to stimulate the optic tract and record from the superior colliculus. At the time of the upper ERP the cat is placed in a dark box, and thereafter the ERP is taken at the times indicated. The frequency of the ERP increases slightly and the decay rate decreases during dark adaptation (Pickering and Freeman 1967).

18.6. Fitting curves and testing the results

Each basis function is specified by an algebraic term containing one or more coefficients, and the sum of basis functions is generated by the sum of the terms. The next step is to fit the sum of basis functions to an AERP. This entails finding the proper values of the coefficients by curve-fitting. A number of standard computer procedures are available for this. The most useful procedure is non-linear regression (Golub and Pereyra 1973; Kaufman 1975; Freeman 1975, pp. 103-106), in which all of the coefficients (usually 4-10) are simultaneously adjusted. There are 5 tests that must be passed for provisional acceptance of a fitted curve.

Fig. 7. Damped sine wave curves are fitted to the ERPs from the superior colliculus (above) and hippocampus (below) to measure frequency, decay rate, phase and amplitude. (From Freeman 1980.)

The first and most basic criterion for determining the optimal values for the coefficients is minimization of the sum of squares of the differences between the AERP and the sum of basis functions. The ratio of mean square differences to mean square data should be well under 5%. Typically a first attempt yields an unsatisfactory match, and repeated trials must be made while adding new terms or changing existing terms. A successful outcome of the first test is when the residue on subtracting the fitted curve from the AERP resembles white noise or, better still, the appearance of records of averages over trials on which no event occurred.

A second and essential test of such a fitted curve is to fit the same curve to each of the ERPs that comprise the AERP. Here the basis functions are adapted, in that the values of their coefficients are optimized to fit each subaverage. If this cannot be done, the most likely reasons are that the basis functions were not properly chosen to reflect the variation, or that averaging was sufficiently destructive to remove certain components, or that one or more components occurred on some trials and not on others. Either a new matching set of basis functions must be sought, or the ensemble must be divided into subsets on the basis of the results on hand.

The third test of whether the decomposition is effective lies in study of the relations between sets of coefficients and concomitant measurements of behavioral variables. When the experimental conditions have been fixed, we should observe ‘spontaneous' variation of the coefficients from trial to trial. When a behavioral factor is varied, some further variation of coefficients should appear. In both cases the values for each coefficient should be 'well-behaved' in the sense that they fall into a unimodal distribution (not necessarily normal), and that we can see patterns of covariation within the matrix of coefficients.

The patterns of covariation between coefficients are best explored with a graph of each pair to determine the ranges and distributions, and the form of the relation.

Three examples are shown in Fig. 8. The upper frames show the AERPs from 3 cats (triangles), the fitted curves (solid), and the dominant basis function and clutter (dashed curves). The lower frames show the values for the frequency and decay rate of the ERPs plotted each as a triangle. The regression line for each set is flanked by the 95% confidence intervals (dashed curves). Also, negative correlations hold between amplitude and both frequency and decay rate; the arrows indicate the direction of increasing ERP amplitude. Whenever they are needed, suitable transformations should be sought to normalize the distributions and to linearize as many of the relations among the coefficients as is feasible (see Ch. 16). Then factor analysis with Varimax rotation is used to characterize the dimensionality of the patterns (Emery and Freeman 1969) (see Ch. 5 and 16).

Fig. 8. Three examples are shown in the upper frames of AERPs (triangles), the sum of basis functions (solid curves), and the dominant component and clutter (dashed curves). The recording sites are monopolar on the surface of the olfactory bulb; the stimulus sites are the primary olfactory nerve (PON). The frequencies in radians/sec are plotted against the decay rates in I/sec of the component ERPs in the 3 lower frames (see Appendix B). The ERPs in the left frame are collected concomitantly from 64 electrodes in an 8 X 8 array; the decay rate decreases in proportion to the increase in response amplitude (arrow) at fixed stimulus intensity. The stimulus is increased stepwise for the other two sets of ERPs (increasing amplitude in the direction of the arrow). Each solid line shows the regression line; the dashed curves show the 95% confidence limits. The variation in the slope of the relation between frequency and decay rate is analyzed elsewhere (Freeman 1975, Ch. 6); the product moment correlation between frequency and decay rate does not by itself adequately describe the relationship. The small circle on the abscissa shows the average decay rate of the bulbar ERP under deep anesthesia, when the EEG is suppressed and the bulbar evoked potential becomes non-oscillatory; this 'open loop' response is discussed in Section 18.7 and is illustrated in Appendix A.

The fourth test is to determine whether the same set of basis functions can be applied to the ERPs from several subjects, and whether the same factors emerge from the matrix of coefficients. Here we must come to terms with between-subject variability. It should now be possible to fit the ERPs of considerable variety of shape with the same sum of basis functions, though the mean values of the coefficients differ significantly. Alternatively it may be necessary to add, delete, or change terms within the sum. It is reasonable to postulate that the same neural mechanisms operate in all subjects under study, but that their manifestations in the form of ERPs may undergo distortions in differing ways or be overlapped by subsidiary events in unpredicted ways. If we can identify one or more components that are always present (though variable in shape), we may characterize this subset of basis functions as a signal and treat the other sporadic basis functions as background 'clutter.'

There is another way of looking at the problem of extracting behaviorally relevant information. The procedure of curve-fitting is essentially the design of a matched filter to detect a desired signal against a white noise background. Any one basis function or any subset of them within a set comprising a fitted curve can be specified as a filter that is matched to a signal. Then the remaining basis functions serve as a filter that is matched to events in the ERP other than the desired signal.

For example, the ERP in Fig. 4 has 2 damped sine wave components. The larger one is a manifestation of the response of an intracortical negative feedback loop in the olfactory cortex to single-shock stimulation, that is manifested by a reverberation between sets of excitatory and inhibitory neurons. This is the signal. The smaller component is the manifestation of a number of complicated feedback loops that include or project onto the intracortical loop, including neurons in the olfactory bulb, nucleus and tubercle, and possibly recurrent projections from the cortical efferent systems. The contributions to the dynamics of those longer pathways are not easily distinguishable from each other, but the sum of their effects can be approximated in form by a damped sine wave. This is the 'clutter.' When the cortical dynamics change, both the signal and the clutter vary, and the two filters are adapted (by curve-fitting) simultaneously to reflect the change. To a certain extent the variations in the two filters are correlated, as they should be, because they both involve the intracortical feedback loop.

We can expect that for differing stimulus sites or different subjects the category of the signal will be invariant. That is, the equation for the dominant component is always a damped sine wave, though the values for the coefficient are variable. However, the form of the clutter may vary, meaning that the equation for the sum of the other subsidiary basis functions may vary as well as the values of the coefficients (Fig. 3, left and right frames).

The use of subsidiary basis functions can be regarded as a pre-whitening transformation of the ERP, before the AERP filter is adapted to ERPs in order to recover the desired signal. In the interest of economy, if the ranges of the frequencies of the filters for the signal and the clutter are sufficiently far apart, the coefficients of the filter matched to the clutter in the AERP may be fixed, so that only the coefficients of the filter for the signal are varied during adaptation to the signal by non-linear regression. The variance of the dominant component in relation to behavior may still be optimally preserved.

The fifth and ultimately most important test is correlation with behavior. The underlying hypothesis is that the electrical signs in each part of the brain depend on neural activity of that part, and that some albeit undefined aspects of on-going behavior also depend on that part, because the intact brain functions as 'an organic whole.' Our own choice for a 'behavioral touchstone' is to measure the rate of work done by hungry cats for food in an ergometer over a period of 10 sec that is required to derive each ERP (Freeman 1962b). The ERPs are fitted with curves consisting of 2 damped cosines, yielding values for 8 coefficients. From a total of 30 trials, 30 sets are collected from each cat. The sets of 8 X 30 values are introduced into a set of 30 simultaneous linear equations with 8 ERP coefficients as the independent variables and the rate of work as the dependent variable. Then we calculate the multiple correlation coefficient R² to evaluate the fraction of the total variance of rate of work that can be predicted from the ERPs. The criterion for the optimization of the process of measurement is the maximization of R² (Freeman 1964). Upon completion of the development the fraction ranges from 0.19 to 0.55 among 7 cats; on the average 36% of the variance in rate of work for food is predictable from the ERPs of one prepyriform cortex. We conclude that the ERP data are optimally cast in forms suitable for neural modeling.

From this example, the desired properties of the behavioral touchstone are that it be an event which extends through the time course of an evoked potential (or a set of them for an ERP), which is continuously variable in magnitude, and for which the magnitude is determined by brain function, whatever might be the antecedent conditions (e.g., food deprivation). Further, the real numbers in which events are expressed should be normally distributed on an ordinate scale. The development of such a criterion may tax the ingenuity, but the reward is reliable data.

18.7. Models of the dynamics of neural masses from ERPs

In this section we consider procedures used to explain the shapes of ERPs in terms of the anatomical and physiological properties of neurons. The properties we rely on are well known to hold for almost all types of neurons; no special or unusual characteristics are needed. In the spirit of an introductory review we restrict consideration to ERPs that accrue from electrical stimulation, that conform to superposition under paired-shock testing, and that are fitted with sums of linear basis functions. In this circumstance each basis function is the solution to a differential equation that can be found in a table of mathematical functions.

The differential equation is not in itself a model; it must be structured so as to represent some aspect of the dynamics of the neural system generating the ERPs. For example, let us suppose that one of the basis functions for potential v is a damped cosine as a function of time t (Fig. 8)

where v₀ is the amplitude in microvolts, w is the frequency in radians/sec (2π times the frequency in Hz), f is phase in radians, and a is decay rate in 1/sec (the time constant is 1/a msec). The generic differential equation for the damped cosine is

where v (with dot above) = dv/dt, I is the magnitude of an input impulse ∂(t) including a scaling parameter, and b₁ and b₂ are coefficients that are interpreted in a worked example below.

It is not sufficient to state the obvious that an 'oscillator' exists in the neural system; we would like to determine whether it is within or between neurons. In order to learn this we must set up two hypotheses in the form of models, one to describe an oscillation of non-interactive neurons, the other to describe the oscillation among sets of interactive neurons. We must use each model to make predictions about ERPs and PSTHs in order to validate or reject the hypotheses. We know on the one hand that under restricted, controlled conditions an isolated nerve can generate a damped cosine fluctuation of membrane potential in response to a brief shock (as distinct from a periodic pulse train). The Hodgkin-Huxley equations admit to corresponding solutions, but they are far from robust. On the other hand we know that cortical excitatory neurons that are excited by an afferent volley, in turn excite inhibitory interneurons, and thereafter receive synaptic inhibition. A typical diphasic cortical evoked potential manifests an extracellular compound EPSP followed by a compound IPSP, both generated by the excitatory neurons, the IPSP mediated by the interneurons. This phenomenon is quite robust. Our aim is to devise a model to simulate a hypothetical process in which the diphasic ERP manifests the impulse response of a neural negative feedback loop. That is, the feedback interaction within masses of excitatory and inhibitory neurons is described by a linear differential equation, and the solution of the equation with impulse input must give us a curve that closely fits the ERP.

The best way to analyze the neural mass is to open the loop and measure separately the passive membrane time constant a_e of the excitatory cell mass and that of the inhibitory cell mass a_i. The loop is opened by inducing very deep anesthesia to suppress the EEG and render the ERP non-oscillatory. What we observe is the compound EPSP of the excitatory cells without the following IPSP. With further search we then find a means to deliver an electrical stimulus to excitatory axons ending synaptically on the inhibitory interneurons to elicit the interneuronal compound EPSP (not the IPSP). We fit curves to these ERPs in the same manner as for other ERPs (Appendix A).

For illustrative purposes we begin with a first order approximation; we ignore the rise time and consider only the decay. The potential V_e as a function of time for the EPSP of the excitatory mass is

and that for the inhibitory mass v_i is

where as before I is the magnitude of an input pulse, a_e and a_i are rate coefficients measured by regression, and v_e and v_i represent the time courses of the active states of the two neural masses. These are known as the state variables of the two neural masses making up the hypothetical neural system being considered (Freeman 1975, pp. 34-36).

Experimentally we find that the rate constants of the two masses of cells are equal, so we set a = a_e = a_i, and the differential equations in the open loop state become

For the waking 'closed loop' state the inhibitory mass receives its input not from a stimulus directly but from the output of the excitatory mass.

where k_e is a forward gain coefficient introduced to express the strength of action of the excitatory mass onto the inhibitory mass. The excitatory mass receives both the initiating pulse input and the output of the inhibitory mass.

where k_i represents the strength of action of the inhibitory mass onto the excitatory mass, and the sign (-) denotes the fact that the action is inhibitory. We combine Eqs. 7 and 8 to find an expression for each of the two state variables. We do this by taking the derivative with respect to time of both sides of Eq. 7, solving Eq. 8 for v_e, and using Eq. 7 to define V_e

and similarly for v, where K_n = k_ek_I is defined as the negative feedback gain. Eq. 9 should be compared with Eq. 2 in order to see the way in which the coefficients b₁ and b₂ have been evaluated.

The two roots of Eq. 9 by the quadratic equation are

so the solution for the inhibitory state variable is

By using Euler's formula,