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,

and by defining a term for phase f,

we predict the output of the inhibitory mass for impulse input to the excitatory mass.

By the same procedure the output of the excitatory mass is

For both state variables the decay rate of the predicted outputs in Eq. 1 is the same as the open loop rate constant a, and the frequency w = K_n^0.5 is the square root of the negative feedback gain. Importantly, the oscillation of the excitatory mass (f_e = 0 radians) leads that of the inhibitory mass (f_i = π/2 radians) by one quarter cycle, if k_e = k_i, or with some further lead or lag if they are unequal (Appendix B).

This result means that a necessary (but not sufficient) condition for the interaction hypothesis (as distinct from the single neuron hypothesis) to hold in cortex is that the oscillatory activity of excitatory and inhibitory neurons must have the same frequency and decay rate in each time epoch, but that the one must lead the other on the average by a quarter of a cycle. This condition has been found to hold in the olfactory bulb, prepyriform cortex, and hippocampus (see Appendix A for an example). It should be noted that the frequency and decay rate serve to evaluate the closed loop rate constants of the system of excitatory and inhibitory neurons. This is distinct from the open loop rate constants defined in Eqs. 3 and 4, that were evaluated in the open loop state under deep anesthesia.

However, the solution of this elementary model further predicts that the decay rates a are equal for both open and closed loop responses. This prediction immediately fails; for example, the decay rate of ERPs in the bulb and cortex under deep anesthesia is near -220/sec, whereas the decay rates of ERPs from waking animals range from near zero to less than -150/sec (Fig. 8). This result means that the closed loop gains K_n cannot be estimated directly from the frequencies W₂ Of ERPs by use of this model.

The basis for the failure of this elementary model is 3-fold. First, the first order approximation for the open loop EPSP includes only the passive membrane decay rate and not the rise time. It is essential for the theoretical studies to use at least a second order differential equation in order to include the cable properties of dendrites (Appendix A), and for modeling ERPs it is desirable to use a third order approximation that includes synaptic and axonal delays as well (Freeman 1975, pp. 314-321). The resulting equations are cumbersome, but the procedures for solution are basically the same.

The second reason is that the use of a single gain coefficient k_e or k_i does not allow us to separate the contributions to overall gain of the operations at synapses and at trigger zones. Each mass of non-interacting neurons in parallel (a KO set see Section 18.1) performs 4 operations in series: conversion of pulse input to synaptic current, dendritic integration, conversion of the integrand to pulses at trigger zones, and axonal transmission. Of these only the conversion at trigger zones is non-linear. For masses of neurons this function is sigmoidal (Appendix C); for any specified low-level stimulus the function can be approximated by a straight line with slope g. This action of replacing the non-linear operation in an equation for the KO set is known as linearization; when the replacement is done over a range of inputs with a different value of g in each part of the range, the procedure is known as piece-wise linearization. We find experimentally that this non-linear function is the rate-limiting step in neural feedback loops that holds the other operations within their linear domains. For this reason we represent the synaptic gain from the i-th to the i-th neuron or subset of neurons by a coefficient k_ii, and within local spatial domains of short axons we combine the time delays of axons and dendrites into one linear differential equation. The overall forward gain is the product k_ijg_j.

The utility of this formulation lies in 3 experimental findings, namely that the open loop rate coefficients do not vary in relation to behavioral changes or to the induction of anesthesia, that the coefficients g_j vary in relation to level of arousal or motivation and also to stimulus intensity, and that the coefficients k_ij change in relation to conditioning, including associational learning and habituation. In brief, the closed loop rate constants from ERPs change markedly in relation to behavioral variables, but this is the result of selective changes in synaptic and trigger zone sensitivities and not of changes in the open loop rate constants.

The third reason is that negative feedback is not the only type of feedback in cortex. To be sure it is a major feature of the excitatory and inhibitory neurons in the cortex, but there are also varying degrees of mutual excitation and mutual inhibition that constitute two forms of positive feedback. An effective approach to this problem is to describe the mutually excitatory mass of cells as a KI, set, the mutually inhibitory mass as a k_ij set, and their interaction by negative feedback in a KII set as the basic dynamic element of cerebral cortex (see Section 18. 1). The topologies of connection of the KO, KI, and KII sets as they are found in the olfactory bulb are presented in Fig. 9. The key step in the analysis of the KI set is to treat it as the interaction of two KO sets having the same sign of action. This allows representation of the feedback gain by a coefficient K_e = k_ee² or K_i = k_ii² and it avoids the inclusion of self-excitation and self-inhibition, which are inadmissible on anatomical and physiological grounds. The negative feedback gain is K_n k_eik_ie as before.

Fig. 9. At left is shown a schematic flow diagram for the olfactory bulb. The receptors (R) form a KO_e set that is excitatory to the bulb (+). A KI_e mutually excitatory set represents the periglomerular (P) cells in the outer bulb. A KII_MG set represents the mitral (M) and granule (G) cells in the inner bulb, with output from the bulb to the cortex by mitral axons. In the center is shown the equivalent circuit when the KI and KII sets are replaced by KO sets with the appropriate 3 types of feedback (+ +, - -, and + -). At the right are shown the prototypical time courses for the KO, KI and KII masses. e, excitatory, i, inhibitory; t, transmitting; r, receiving. (From Freeman 1975.)

The minimal representation of the dynamics of the lumped KII set is by 4 coupled second-order differential equations. Two state variables are required for each of the 4 component KO sets; these are to represent its activity as a function of time in the dendritic wave mode v and in the pulse mode p. The relations between these two forms of the active state are experimentally analyzed from extended time series (10-20 min) of digitized (1 msec) EEG activity from a defined KII neural mass and concomitantly recorded pulse activity from a single neuron in the same KII mass. The amplitude histogram of the entire set of digitized EEG values is formed. The conditional pulse probability is defined by the ratio of the number of pulses that accompany the EEG amplitudes in a given bin of the histogram to the number of amplitude entries into that bin. It is determined over a range of EEG values from -3 to +3 S.D. It is normalized by dividing each value by the mean pulse probability. The procedure is repeated with successive time lags in both directions between the EEG trace and the pulse train, in order to allow for intrinsic delays within the generating mechanism (Freeman 1975). The sigmoid curve (Appendix C) is fitted to the normalized conditional pulse probability in order to evaluate the non-linear gain. In the linearized model the state variables in the pulse p and wave v modes are related by p = g_jv.

When this function has been evaluated, the state variables may be cast either in the wave mode or in the pulse mode. In the olfactory bulb and the KII_MG mass the activity of the excitatory mitral (M) cells is normally observable only in the pulse mode in the form of pulse trains and poststimulus time histograms (PSTHs), and the activity of the inhibitory granule (G) cells is observable only in the wave mode in the form of EEG waves and ERPs. These are the forms in which events are obtained for testing and evaluating the gains within the KII set (Figs. 3 and 8 and Appendix A).

The evaluating procedure consists of setting a reference gain called K₀ at which K_n = K_e = K_i in the linearized KII equations and solving them for a median value of K₀ = 2.0 (small rectangle in Fig. 10). Then the 3 linearized gains are changed in small steps according to certain rules expressed in subsidiary equations that we call describing functions (Freeman 1975). The frequency and decay rate of the solution change in a consistent way that allows us to plot a curve of decay rate vs. frequency (Fig. 10). This is repeated for different values of K₀ giving families of curves such that, if we have the frequency and decay rate from an ERP or PSTH, we can estimate the linearized gains from a graph. The describing functions are devised in accordance with hypotheses about the effects on the KII gains of changing the stimulus parameters, of administering anesthetics or other drugs, or of training animals with implanted electrodes so as to manipulate attention, motivation, or other behavioral variables (Freeman 1975). Two examples are taken up in the next section. To reiterate, the most important result here is that frequencies and decay rates of ERPs and PSTHs serve to evaluate neural sensitivities, provided that a valid model is used.

Fig. 10. The triangles show the frequencies and decay rates from sets of ERPs that are derived by adapting to them the matched filters from their AERPs. The curves represent extrapolations from evaluations of the frequencies and decay rates of repeated solutions of the equations of the lumped, linearized KII model for differing values of the gain coefficients that are specified by two sets of describing functions (see Ch. 6 in Freeman 1975). Right frame: the data show the effects of increasing the stimulus intensity in small steps. One factor suffices to describe the results for each subject and session, but two other factors are needed to account for differences between subjects and sessions. Left frame: the data show the values from one subject and session with stimulus intensity fixed near threshold ('spontaneous' variation). Two main factors are apparent from factor analysis of the coefficients from the adapted filters; their independent effects are shown by means of the family of curves. Each curve is generated for a value of a reference gain K₀ that is increased in steps of 0.5 (small rectangle at K₀ = 2) and that represents increasing arousal. The curve that runs upward to the right is the locus of values for K₀ = K_n = k_e = k_i with its origin at w = 0 and a = -220/sec. The position of points on each curve is determined by a mutually excitatory gain k_e that represents increasing selective attention with decreasing frequency. In both graphs the vertical line shows the border for stability. The diagonal line at the top shows a boundary that reflects onset of another form of instability; it occurs by the existence of a pole on the abscissa to the right of the origin for certain combinations of the gain values for K_n K_e and K_i. The generic properties of this family hold over subjects and sessions, but the degree of curvature varies in accordance with another parameter in the linearized KII equations (Freeman 1975, pp. 355-378).

18.8. Using ERPs to evaluate the stability of neural masses

Each part of the brain is, directly or indirectly, under ceaseless bombardment from sensory input, as well as from the activity generated by internal dynamic processes. Each part must have the property of maintaining its integrity and status quo under most kinds of external and internal activation, but it must also be capable of undergoing changes to different states when the appropriate circumstances arise. In other words, it must be stable under most kinds of perturbation, but it cannot be absolutely stable, or it would be incapable of participating in the evolutionary processes that underlie behavior. In this section we take up the problems of how to make quantitative estimates of the degree of stability of a neural mass by measuring its ERPs, how to determine the conditions in which instabilities might occur relating to arousal and attention, and how to predict the forms they might take.

Any number from half a dozen to over 100 ERPs can be measured from their common AERP. Each ERP yields a value for each of the several coefficients representing the amplitudes, frequencies, phases and decay rates of the component basis functions. Of these coefficients the most crucial in respect to stability are the decay rates of exponential curves and damped oscillations, or more particularly for the latter, the ratio of frequency to decay rate. The most useful graphic display is the decay rate of each basis function plotted on the abscissa against its frequency on the ordinate as in Figs. 8 and 10, with zero frequency for exponential terms. Fig. 11 shows prototypes of linear basis functions with various frequencies and decay rates, each positioned in the approximate location of its point in these coordinates.

Fig. 11. The relationship is shown between wave shapes of ERP components and the values for frequency w in radians/sec and decay rate a in 1/sec. Each trace is an example of a component 100 msec in duration that is centered at the point for its w and a. For example, the upper sine wave has w 200 radians/sec (31.8 Hz) and a = 0/sec. (From Freeman 1975, pp. 278-284.)

The outstanding feature of this display is the fact that the sign of the decay rate shows whether the neural mass has responded to an input in a stable or an unstable way. If the decay rate is positive, the output in response to, for example, a single shock does not decay to zero but continually augments. If the output for a brief input is a steady-state voltage or a maintained oscillation, the neural mass is marginally stable. Only if the decay rate is negative and the response decays to the prestimulus baseline after the input has ended can the system be described as stable, whether or not the response is oscillatory.

The experimental conditions in which ERPs are taken preclude direct observation of an instability, because we require that the potential return to baseline after each stimulus, and before we give the next one. Therefore, all of the decay rates must be negative. However, the more strongly negative the decay rates are, the more stable the system is, and conversely the closer a point representing the frequency and decay rate of an ERP lies to the line for zero decay rate, the less stable is the system. The pattern of the collection of points from a set of ERPs may suggest the experimental conditions in which an instability may arise, and the frequency or rate of change that will then be manifested.

There are two kinds of instability revealed in such patterns that are particularly important. One predicts the onset of an epileptiform seizure and the other one the appearance of EEG waves.

A set of values from ERPs is shown in Fig. 10 (right frame) by the triangular plotting symbols. As in Fig. 11 the vertical line shows zero decay rate. The set is taken while increasing the current strength of the stimulus in small steps. There is successive increase in amplitude, decrease in frequency, and increase in the negativity of the decay rate of the damped cosine wave fitted to the ERPs as indicated by the downward arrow. This result shows that minimal stability occurs with minimal input intensity and minimal response amplitude. When the peak amplitude of the single shock evoked potential rises above the level of the ongoing EEG, the response of the neural mass indicates that the system is made more stable in proportion to the magnitude of the impulse input.

This is a common response configuration observable in all parts of the cerebrum. It manifests a simple and reliable neural mechanism for maintaining the stability of each part in the face of massive but brief overload. It occurs when the system is saturated by excessive input. The initial excitation of excitatory neurons by a large impulse input is virtually unlimited, and their output is given pari passu to inhibitory interneurons. However, the effectiveness of the inhibition back onto the excitatory neurons is limited to the amount of pre-existing background activity. Once the excitatory neurons have been shut off they cannot transmit their degree of further inhibition; they saturate at threshold, because pulse rate goes to zero and cannot be negative.

This mechanism can be modeled by a negative feedback system in which the feedback gain is reduced in proportion to the response amplitude and thereby the driving input. The curves in Fig. 10 (right frame) show the solutions for such a model (the KII set) for a set of values of the level of pre-existing background activity expressed in K_O. The rectangle shows the reference point for K_O = 2.0; each curve predicts the frequencies and decay rates of a set of ERPs with increasing stimulus intensity and with a persisting level of background activity between successive responses. The curves form a family in which the system is made more stable with an increase in the input intensity whenever the single shock evoked potential rises above the EEG activity.

The same equations that predict the form of the ERP also predict that at very high stimulus intensity one of the components of the ERP should be an exponential term, for which with increasing input intensity the decay rate goes to zero and then becomes positive. That is, the system is predicted to become unstable. This has been observed under massive repeated single shock stimulation of the hippocampus, prepyriform cortex, or septum, leading to the onset of a petit mal-type seizure in waking animals (Freeman 1962a). Further quantitative modeling of this phenomenon is needed, as well as study to determine whether it also occurs in neocortex (Freeman 1986, 1987).

It is now obvious that to increase the stimulus intensity in order to make an ERP more visible is to render abnormal the neural dynamics under study. Yet the tendency of the ERPs to vary in their shapes increases at near-threshold stimulus intensities. We take advantage of this variation by collecting sets of ERPs under fixed stimulus intensity and minimal behavioral variation. The sets of points show the 'spontaneous' or background variation, and thereby they reveal something about the properties of the neural mass. An example in Fig. 10 (left frame) shows a set of values for the frequency and decay rate from the olfactory bulb of a cat. As before, all of the decay rates are negative, but here we find that the pattern of variation is not the same as that with changing the stimulus intensity.

Two factors dominate the variations in the shapes of these ERPs. One is the level of arousal or motivation that can be increased by such measures as food deprivation. This results in increased background activity that is manifested by increased rates of background single unit firing, by augmented amplitude of the EEG, and by corresponding changes in ERPs and PSTHs from the bulb and cortex. These effects are simulated in the output of the linearized KII model by increasing the reference gain coefficient K₀ above 2.0. By the criterion of decay rate this decreases the stability of the system. Each curve of the family of curves in Fig. 10 (right frame) is generated for a particular value of K₀ increasing to the right (Freeman 1975, pp. 370-378).

The second factor is one that we identify with the degree of attentiveness to the evoking stimulus. We do this by implanting electrodes in the olfactory tract and prepyriform cortex of cats, training them to respond to electric stimuli used as conditioned stimuli under reinforcement, and applying factor analysis with varimax rotation to the sets of coefficients of ERPs taken before training, after training, and after extinction of a conditioned response (Emery and Freeman 1969). The factor pattern that emerges and that is related to attentiveness by analysis of variance consists of concerted changes in phase, frequency, amplitude and decay rate of the ERPs. We simulate this pattern change in the impulse response of the KII set by increasing the values for the mutually excitatory gain coefficients k_ee. other change in the model suffices (Freeman 1979b). Each curve in Fig. 10 (left frame) is generated by increasing k_e, for fixed K₀ causing increasing response amplitude and decreasing frequency (downward arrow). The decreasing negativity of decay rate implies decreasing stability of the model; correspondingly the less negative decay rates of the ERPs imply lessened stability of the KII mass in respect to the evoking electrical stimulus under reinforcement as learning takes place.

We also find that the curves in the left frame correspond to the factors traced out by sets of frequencies and decay rates from ERPs undergoing 'spontaneous' variation (Fig. 10, left frame, triangles). We infer that the expressions of independent properties of the neural KII mass relating to arousal and attention (whether or not the subjects are attending to the stimulus) are subject to random variations in respect to the timing of the stimulus trains (Freeman 1975, pp. 386-390).

These curves differ from those in the left frame in an aspect that has major importance. Whereas all the curves in the left frame cross the vertical line from right to left with increasing response amplitude (arrows), some of the curves in the right frame cross the vertical line from left to right with increasing amplitude. Those that do so have relatively high values of K₀ which are associated with increased arousal, and the parts of the curves that cross to the right have relatively high values of k_ee, which are associated with selective attention. According to the stability criterion these curves imply that the combination of arousal and attention may cause the KII mass to become unstable, when the stimulus that is attended actually arrives. That is, oscillation might be initiated that would augment rather than decay.

This cannot be observed with the evoked potential method; simulations with the non-linear KII model (to be discussed in the next section) show that single shocks are too brief to allow this instability to occur. Even if it did, the fluctuations in phase and frequency of its onset would lead to destructive averaging, causing the appearance of a faster decay rate for the ensemble average than those for single responses. However, in aroused, attentive animals we observe that with each inhalation the inputs to the bulb and prepyriform increase owing to augmented receptor input, and we know that then an oscillatory burst commonly occurs in the EEG of the olfactory bulb and cortex (Fig. 12). The modal frequency of the oscillation found in cats is 40 Hz (250 radians/sec divided by 2π), which corresponds to the convergence point of the curves in Fig. 10 (left frame) where they re-cross the vertical line from right to left with increasing amplitude. (These curves are derived from ERPs in cats; the simulation of the modal gamma activity in rats and rabbits that occurs more often in the range of 50-70 Hz is done by adjusting the gain values.) Intuitively this suggests that a small fluctuation that is initiated in this condition might grow in amplitude at the indicated starting frequency, but that as it grows its frequency might decrease in concert with the increase of amplitude, until it re-stabilizes at a steady frequency and amplitude of oscillation. During exhalation there is cessation of receptor input, allowing the oscillation to end. That is, the neural mass can be inferred to switch from one stable state to another and back again with each cycle of respiration. This pattern of function cannot be described quantitatively with linear approximations. It requires simulation with an explicitly non-linear model. We take this up in the next section, where we explore the physiological basis for gamma activity in the range of 35-90 Hz.

Fig. 12. A: the upper trace shows respiration of a waking rabbit measured with a pneumograph; the lower trace is an EEG from the olfactory bulb. Characteristically there is a surface-negative wave that peaks after each inhalation; riding on the negative crest is a brief burst of gamma activity. B: the upper trace shows the air flow through the nose of an anesthetized cat with a respirator connected to the distal end of the trachea; respiration occurred at an independent and higher frequency than that shown for air flow. The lower trace shows the bulbar EEG with the slow wave removed by filtering. The burst occurs with air inflow and not with inhalation, showing that the burst depends on sensory input and not on a centrifugal pathway from the brain-stem to the bulb. (From Freeman 1976.)

18.9. Non-linear dynamics of the EEG

We turn now to the EEG, noting that in the half century since the electronic amplifier first gave easy access to this phenomenon, empirical exploration has done little to resolve its baffling origins and complexities. We need theories expressed in models to give us explicit hypotheses on what to took for and how to took. We have used the modern computer to analyze our data with digital filters, integral transforms, spectral displays, non-linear regression, and multivariate statistics; we must also use it to express our hypotheses and predictions in mathematical forms.

We restrict consideration to the full KII set. We use the lumped KII set (Fig. 9) as the kernel of a new model in which the linear dynamics of each component KO set is approximated by a second order equation or equivalently by 2 first order equations. We make two advances. First, we retain the non-linear wave-pulse function to present the operation at the trigger zones of neurons in each component KO subset, instead of replacing it by a straight line approximation. Second, in recognition of the fact that each neural mass is distributed in some part of the brain, we form an array of kernels (non-linear KII subsets) in the form of a surface, and interconnect them by mutually excitatory and mutually inhibitory synaptic lines as shown in Fig. 13. The input and output pathways are likewise distributed but are locally determined in respect to KII subsets as already defined for the lumped linear approximation. The equations are summarized in Appendix C.

Fig. 13. The topology of connections for the bulbar KII set is shown for the distributed case. Each component KII subset has the same form as that in Fig. 10; it is replicated in one-dimensional circular arrays with each element connected to every other by mutually excitatory and mutually inhibitory lines (not shown in full) or as two-dimensional arrays in a toroid (Freeman 1979c). The subscripted parameters are forward gain coefficients having origin and termination on excitatory (e) and/or inhibitory (1) neurons. An input line is designated 1. The superscripts denote synaptic lines that are modifiable by learning.

There are no closed or analytic solutions known to exist for the non-linear distributed KII set, nor are any likely to be found. Solutions must be found by numerical integration. For this purpose the model is expressed in an 8 X 8 array of 64 kernels, each having 8 first order equations, so that the dynamics is described by 51.2 coupled first order ordinary equations instead of by partial or integrodifferential equations. Owing to the use of this small number of kernels in this discrete approximation, we must suppress instabilities due to edge effects by using a toroidal boundary condition. In effect every kernel is connected to every other in all directions over the surface. For simplicity, propagational delays are omitted; interactions have no dependence on distance over the surface of the array. Again for simplicity the initial conditions for amplitudes and rates of change are all set to zero, and the synaptic strengths (the k_ijs) are made spatially uniform. The input is either a set of impulses given to selected input channels to simulate low-level orthodromic or antidromic electrical stimulation, or (in the example from the olfactory system) a surge designed to simulate receptor input to the bulb during an inhalation.

In order that it be admissible, the model must meet at least 3 criteria. The first is that it be stable over a sufficiently wide range of input and of parameter settings to allow ERP and EEG simulations. This aim is achieved by including several features, such as: dynamic range compression at the input corresponding to an operation done by the periglomerular interneurons in the outer bulbar layer (Freeman 1975, pp. 291-298, 1979b), the proper boundary condition, a degree of symmetry in the strengths of mutual excitation and inhibition (k_ee and k_ii), and minimization of the spatial divergence of negative feedback k_ei and Kl_e (diagonal connections in Fig. 13 are unnecessary). The second criterion is that the model replicate all of the known properties of ERPs in relation to input parameters, behavioral manipulations, and drug effects. This is done by use of impulse input, adjusting the synaptic and non-linear gain parameters, and measuring the two kinds of output (waves and pulses) with the same procedures used to measure ERPs and PSTHs. The third is that the model serves to generate the main varieties of EEG wave-forms (Fig. 12). Examples of simulated bursts are shown in Fig. 14 (Freeman 1979b,c, 1986, 1987).

Fig. 14. Various shapes of burst are simulated with the non-linear KII set by using slightly different values for the gain parameters. The right lower frame shows a limit cycle induced by sustained simulated input (see also Gonzalez- Estrada and Freeman 1980). (From Freeman 1979b.)

The class of coupled non-linear oscillators to which the KII set belongs has certain properties that are crucial for understanding the EEG. Each system has multiple stable states. Each stable state is characterized by parameter values and a domain of input, such that when the system is placed within the domain it will stay within that domain unless and until it is forced out by a change of input or by a parameter change. If then left to itself the system will settle into one of 3 types of output: either it will go to rest and stay there, or it will go into regular oscillation, or it will go into a restless aperiodic fluctuation in some limited region of its output range. If further inputs are given within the designated domain the system will settle back to its designated end state, manifesting what is called an 'attractor.'

The 3 types are called 'equilibrium,' 'limit cycle,' and 'chaotic' or 'strange' attractor. The input domain over which the attractor exerts its influence effectively is called its 'basin.' Transition from one basin to another is called a 'bifurcation' (Hassard et al. 1981).

Examination of the KII set shows that it has at least 4 equilibrium points. One is a rest state in which the non-linear gain is set to zero. This corresponds to the open loop state imposed by deep anesthesia. Another is a rest state that may be inferred to exist in the conditions of steady state that are used for deriving ERPs and PSTHs. (More will be said about this later.) The other two stable states occur when the aforementioned stability conditions are violated, and the system goes into steady maximal excitation or inhibition. These two states do not exist in the bulb. Although they might represent departures into states leading to or comprising seizures, the KII model is not sufficient to simulate epilepsy. By using a different non-linear gain curve Babloyantz and Kaczmarek (1979) have made some progress in doing so (Freeman 1986, 1987).

A near-sinusoidal limit cycle attractor has been demonstrated by giving the model a step input (Fig. 14D). The amplitude, frequency and spatial patterns of phase and amplitude depend on the input and parameter values. The oscillation decays to rest when the step input is terminated. This type of limit cycle behavior has been induced in the olfactory bulb and cortex by drugs and by hyperthermia (Freeman 1975, pp. 381-390). Recurrent, abortive limit cycle behavior is simulated by phasic input corresponding to the level of receptor input to the bulb during normal respiration (Figs. 12 and 14).

Chaotic activity is generated by the model, but only when it is structurally changed. This is done by forming a feedback connection in each kernel (diagonally between M_t and G_t in Fig. 9 or e2 and i2 in Fig. 13). Then simulated bursts appear to have two or more frequency components with unpredictable amplitudes and phases. Indeed in long stretches the output often resembles the erratic fluctuations of the EEG at rest. Each sequence is reproducible if re-started with the same input and parameters, but small changes in the input or in a parameter lead to increasingly different and divergent wave shapes, though always within certain bounds of frequency and amplitude. Analysis of the KII model is incomplete since the conditions for transition from chaos to limit cycle and back again have not been identified (short of the arbitrary and physiologically meaningless insertion and deletion of connections). The KIII model is discussed elsewhere (Freeman 1987).

A type of output that is very similar to chaos has been generated by giving simulated low-level white noise (random number sequences) as input to the model, when it is in the basin of an equilibrium attractor. The amplitude is limited by the non-linear gain curve, the spectral range is ordained by the bandpass characteristics of the KII set, and the time series is unpredictable and unreproducible whether or not the initial conditions are changed. The similarity poses the questions now to be considered, whether the background EEG manifests an equilibrium attractor under random perturbation or a chaotic attractor, and what criteria might be used to tell the difference.

The distinction is very important for understanding brain function. Chaos is a deterministic and not a random or stochastic process (Helleman 1980; Garfinkel 1983). If the initial conditions and the input that place a system in the basin of a chaotic attractor are reproduced, the system will give the same output no matter how random it might appear. If the input or initial conditions are changed (however slightly), the output will differ increasingly over time and not converge to a periodic or reproducible wave-form. Unlike the random variable in which large deviations occur however rarely, the chaotic variable is bounded; its degrees of freedom are relatively low, its trajectory circumscribed. It can arise from a tendency to instability that is counteracted by a stabilizing influence. For example, if a mutually excitatory population (KI_i set) has a normalized gain greater than unity, it is intrinsically unstable. If its output is delivered to inhibitory interneurons with negative feedback in a KII set, the excitation is counteracted, and restless perturbation results. If we can devise criteria to identify the basal EEG of the bulb as chaotic, we can assign it to a deterministic mode of function of the KII mass; otherwise we must identify sources of random activity in or outside the bulb and demonstrate that they have the requisite properties to drive the KII set in its background state. This problem is considered further in Section 18.13 and Appendix D.

18.10. Spatial analysis of the EEG

The main conclusion to be drawn from non-linear analysis of the olfactory EEG is that the bulbar mechanism is highly non-stationary. The bulb changes its state as it enters a burst with each inhalation and again as it leaves the burst with each exhalation. At a typical respiratory rate for rabbits of 3-5/sec the duration of a given state may be on the order of 75-150 msec. No two successive bursts are identical, nor should they be if there are continual changes in the odor content of air inhaled by the animals. Herein lies a major difficulty in the use of averaging techniques to the extent that they assume stationarity over the time period of averaging. Yet, in order to make estimates of means, variances and confidence intervals, some form of averaging is essential to the measurement process. In this section we present a solution to this problem by introducing the techniques of simultaneous recording of the EEG from multiple electrodes that are placed on or in a single neural structure.

The problem can be viewed in another way. Up to this point we have postulated that a recording is made from a single electrode pair over a time period on the order of 10 sec, which is sufficient to allow data collection for calculation of an ERP or an autocovariance function and power spectrum of the EEG. These procedures give access only to the stationary properties manifested by the dynamics on the average over the whole period but not to the time-varying properties within the period. If we chop the recording into segments on the order of 0. 1 sec in duration, we may be hard-pressed to describe the EEG patterns in those segments with any confidence. We cannot form ensemble time averages in the manner of deriving ERPs, because the phases, frequencies and envelopes of the EEG burst activity vary continually, and we have no time marker by which to set the start of each segment. However, we can form an ensemble time average by recording an EEG segment simultaneously from 60-64 electrodes monopolarly with respect to a distant electrode. We call this ensemble average an AEEG; it serves the same purpose for the nucleus of measurement as does the AERP (Section 18.4) for measuring ERPs.

The technique requires that the electrodes be placed on or in a contiguous area of cortex or nucleus. We will first consider the data and then undertake their measurement. An example in Fig. 15 shows 60 unaveraged traces 100 msec in duration of an EEG burst from a 6 X 10 array of electrodes 0.8 mm apart (a 4 mm. x 7 mm 'window' onto the lateral surface of the olfactory cortex of a rabbit). Each trace is amplified 10K with high and low pass filters at 10 Hz and 300 Hz. The same wave-form appears on all channels with equal numbers of peaks and zero crossings, though with local variations in the envelope of the oscillation, in peak amplitude, and in phase. The commonality holds in spectral analysis, as shown in Fig. 16 by comparing the spectrum of the ensemble average with the average of the 64 spectra from an 8 X 8 array, and it holds also for background segments between bursts both in displays of time series and spectra. It is confirmed by examination of the distributions of the 2016 product moment correlation coefficients between pairs of the 64 traces, and of the 64 coefficients between each trace and the ensemble average of the burst. For example, the grand mean (by z-transform) of the correlation coefficient from 60 bursts by each of 6 rabbits is 0.82.

Fig. 15. Unaveraged EEG recordings are shown of a single burst (100 msec) recorded from a 6 X 10 array of electrodes placed over the olfactory cortex. The commonality of wave-form is apparent over the total window of 4 mm X 7 mm, though there are differences in amplitude and in phase with respect to the ensemble average (AEEG). The amplitude scale is in microvolts. (From Freeman 1978a.)

Fig. 16. The spectrum of the ensemble average (AEEG) of the set of 64 traces recorded from an 8 x 8 array is compared with the average spectrum of the 64 spectra by the FFT. The two curves are in close conformance both in a burst (above) and in an interburst period (below). This reflects the high degree of coherence of the bulbar EEG spatially, whether or not a burst is present. (From Freeman and Viana Di Prisco 1986b.)

We describe each burst as a 'carrier wave' at a certain frequency (in rabbits averaging 68.2 ±12.1 Hz) with spatial modulation of amplitude and phase. We determine the amplitude from the root mean square (rms) of each trace, and the phase by correlating each trace with the ensemble average at 5 time lags, fitting a parabola to the 5 values, calculating the time lag of the extremum, and dividing that lag by the cycle duration at the peak frequency from the Fast Fourier Transform (FFT). Amplitude and phase patterns are displayed in density or contour plots (Fig. 17) either individually or as the ensemble average and standard deviation of sets from 10-20 bursts. The latter show that each animal generates its own characteristic spatial patterns of amplitude and phase. Successive bursts vary unpredictably in shape. No two are identical, yet each is readily identified as characteristic of the pattern of a given animal. The ensemble averages of the spatial patterns of interburst segments correspond closely to those of bursts, but the standard deviations are larger.

Fig. 17. Contour plots are shown in perspective for the root mean square (rms) amplitude of a burst recorded over the anterior olfactory nucleus of a cat, and (below) for the phase pattern measured with respect to the phase of the AEEG. In this case the phase gradient (lead to lag from left to right) reflects the conduction velocity of the underlying lateral olfactory tract. Similar gradients are found in bulbar EEG foci but at lower phase velocities; they have not been found to vary consistently in relation to behavior, and they do not conform to PON velocity or direction (Freeman and Baird 1986). (From Freeman 1978.)

We have devised or adapted statistical procedures based on chi-square and t tests to measure the pattern variability and to assess the significance of pattern differences relating to training animals to respond to odors presented as conditional stimuli; these are described and discussed in a different context (Freeman and Schneider 1982; Viana Di Prisco and Freeman 1985; Grajski et al. 1987).

We undertake decomposition and precise measurement of the multiple EEG segments by curve-fitting in much the same way as for ERPs. There being no analytic solutions known for the full KII set, we choose the cosine as the basis function. One of several variants of the procedure is as follows. For each burst we calculate the ensemble average (an AEEG) and its power spectrum to estimate the peak frequency and its phase, fit a cosine by non-linear regression to optimize the phase and frequency values, and subtract the fitted curve from the average. We repeat the procedure 4 times, and then fit the sum of 5 cosines at the optimized frequencies and phases to each of the 64 traces in order to determine the 5 amplitudes at each channel by linear regression. The result is a 5 X 64 matrix of amplitude values at 5 frequencies, which incorporate over 80% of the burst energy on the average. The graph in Fig. 18 shows the proportion of bursts for which the 5 components cumulatively account for a specified fraction of the energy. There is striking variability in the extent to which energy is compressed into one narrow frequency band.

Fig. 18. The results are shown of measurement by decomposition using cosine basis functions at 5 frequencies on 900 bursts without odor stimulation (control) and on 450 bursts with odor stimulation under reinforcement (CS+). The abscissa shows the fraction of total energy incorporated cumulatively by the 5 components. The ordinate shows the percentage of bursts that have at least the designated fraction of energy. For example, roughly 50% of bursts have at least 0.50 of their energy in the dominant component, but in 10% it is less than 0.25, and in only 10% is it greater than 0.90. On the average about 0.20 of the energy is spatially uncorrelated and is removed by filtering. This fraction increases to about 0.30 in odor bursts (see Figs. 23 and 24). (From Freeman and Viana Di Prisco 1986a.)

To be sure, these components do not imply the existence in the bulb of independent oscillators or harmonics. The variations of their frequencies from burst to burst, the irrational ratios of their frequencies, and their frequency modulation within bursts all argue against this. The aim here, as before, is to incorporate as much as possible of the variance of the wave-forms in order to determine where the correlated variance lies and how its measurement is to be optimized. Here again there is a dominant component; the remaining 4 components are usually sufficiently small that they can be regarded as clutter. Again, the precise measurement of the desired component requires us to design an adaptive filter for the clutter. With this aim in mind we modify the cosine for the dominant component by introducing terms for its amplitude and frequency modulation into that equation

where AM and FM are coefficients providing for linear modulation over time respectively in the amplitude and frequency of the dominant component about its mean amplitude and frequency at the midpoint of the burst t_m.

Fig. 19. Examples are shown of AEEGs of 3 bursts recorded during a control period (C) and during presentation of an odor (TI and T2). At the left are the time series and at the right are the spectra by the FFT. The dots show the digital values. The curves at left show the dominant component that has both amplitude and frequency modulation (treated as linear over time). The faint curves in the right frames show the spectra of this component. The dark curves show the spectra of the residuals after subtraction of the dominant component. This is the clutter plus noise, which is fitted with 4 additional cosines. The reason for devising an adaptive filter for the clutter is to measure optimally the amplitude of the dominant component. (From Freeman and Viana Di Prisco 1986b.)

Fig. 19 shows examples for 3 bursts of the ensemble average (AEEG), their spectra, the fitted curves of the dominant components, and the clutter. The tendency for the energy of the clutter to form spectral clumps justifies the use of cosine basis functions. The 5 cosines routinely incorporate over 99% of the energy in the ensemble average. Occasionally the sum of energy of the components exceeds 100% and rarely 150%, showing that the procedure may recover energy that is otherwise lost to observation owing to interference, or else that the linear decomposition injects structure into the data that does not belong there. Finally, when the components are identified at specific frequencies, we can make corrections for the amplitude and phase distortions imposed by the filters that are used on the time series prior to analysis. An example in Fig. 20 shows the spatial patterns of the 5 components in the form of density plots for 3 subgroups of bursts from the olfactory bulb of a rabbit.

Fig. 20. Density plots are shown for the ensemble averages of components over 60 control, 30 CS+, and 30 CS- bursts from a rabbit. The mean frequency for each set of components is shown above the frames. Below each frame is shown the average fraction as percentage of total energy for the components. (From Freeman and Viana Di Prisco 1986b.)

Measurement of EEGs by decomposition opens the way to spatial pattern analysis, which should be distinguished from use of multi-electrode recording with widely spaced scalp or depth placements to sample from different brain structures simultaneously (Ch. 12). For example, if we have 30 electrodes we may place them all in one cortical area, such as a part of the olfactory bulb, in order to 'digitize' the spatial configuration of monopolarly recorded potential in that part. Or we may place 2 electrodes in each of 15 cortical areas (the bulb, anterior nucleus, prepyriform, amygdala, hippocampus, etc.) in order to correlate the bipolarly recorded time series from multiple spatially distinct structures. The latter is not spatial analysis, because it does not allow the determination of potential as a function of distance. The aim of spatial analysis is to test the hypothesis that within- structure EEG spatial patterns (as distinct from between en- structure correlation patterns) may correlate with behaviorally significant aspects of stimuli, responses, or brain states and operations. However, before we consider such correlations we must deal with the problem of the design and use of spatial filters.

18.11. Spatial spectral analysis, filtering, and deconvolution

The overwhelming advantage that EEG decomposition offers us lies in the solutions it opens up to one of the most vexing problems in EEG work. The electric currents generated by neurons singly or in sets spread outwardly in all directions, often (depending on neuron geometry) to distances greatly exceeding the size of the neurons. This is the biophysical basis for the scalp EEG. Even in proximity to the generating neurons the fine spatial structures of the fields of potential are blurred or smoothed. Worse still, the extracellular compartment in the brain is the common return path for currents from all neurons, and their contributions to potential fields overlap 'and sum in ways related to their geometries and not to their functions. This is known as the 'cocktail party' effect. Solutions to the problems of current spread and intermingling lie in the proper design of spatial filters. It is unnecessary to apply such filters to raw digitized EEGs, because it can easily be done to each of the five 8 X 8 matrices of amplitude values of the components or of the complex gains at peak frequencies from the FFT, which already incorporate most of the variance.

The first datum needed for this work is the interelectrode distance for a regularly spaced array. It is equivalent to the choice of the time interval used for digitizing a time series: if the spectrum of the signal is unknown at the outset, the interval must be made arbitrarily small in order to avoid aliasing from unsuspected high frequency spectral peaks. A useful device is a linear array of electrodes, say 1 X 64, that can be placed across an area of cortex. For small structures such as the olfactory bulb we use a circular array. Fig. 21 (upper frame) shows the 64 rms amplitudes (dots) of a burst recorded from the bulb of a rabbit with an array 4 mm in diameter and an interelectrode interval of 0.2 mm. This is a typical cross-section of the potential of an EEG focus on the bulb.

Fig. 21. Upper frame: the first step in spatial analysis is to determine the proper interelectrode distance that is equivalent to the digitizing interval in time series analysis. For this determination the interval must be arbitrarily small; a I X 64 linear or curvilinear array best serves this purpose. For this example the array is circular with an interval of 0.2 mm over a circumference of 12.6 mm. The dots show the rms amplitudes of a typical burst digitized at 2 msec intervals over a duration of 100 msec. The curve is a 1-D Gaussian density function fitted to the data by non-linear regression and scaled from zero to the point of maximal rms amplitude placed at the origin of the spatial dimension. Lower frame: the potential distribution in the upper frame is simulated by means of a core conductor model of the granule cells in the olfactory bulb (Rall and Shepherd 1968). The distribution of activity density in an active population is represented by a bivariate Gaussian function (inset solid curve) with its origin 0.58 mm beneath the point of origin at the surface on the circular electrode array; the dotted curve shows the fall-off in activity density under the circular segment of the array. The dipole moment density of the granule cell generator is shown to scale beneath the abscissa as if beneath the bulbar surface at a mean depth of 0.58 mm, with discrete approximation by dipoles 0.2 mm apart, and with dipole point separation of 0.17 mm (Freeman 1975). The solid curve is the predicted profile of potential at the surface for an activity density with standard deviation of 1.0 mm. The dots show the effect on rms amplitude of adding 'noise' (random numbers independently for each channel). Each simulated burst is a set of 64 time series 100 msec in length with noise superimposed on a sinusoidal curve with signal to noise ratio of 4: 1. The smooth curve represents noise-free amplitude and is scaled to the maximal simulated array amplitude. (From Freeman and Baird 1986.)

The procedure used is to first apply the FFT to the AEEG in order to find the peak temporal frequency of a burst; second to determine the real and imaginary parts of the complex gain at that frequency for each of the 64 traces; third to apply the spatial FFT separately to the real and imaginary parts; fourth to combine the results in the spatial frequency domain; and fifth to normalize to unit power at 0 c/mm and take the natural logarithm. Results such as that in Fig. 22 show that the spatial spectrum is essentially flat for the rabbit above 0.4-0.8 c/mm and for the cat above 1.0 c/mm (Freeman 1978, 1980). Were we digitizing a time series with its highest frequency at 1.0 c/sec, we would require a digitizing interval <0.5 sec corresponding to one-half cycle duration to avoid aliasing; correspondingly we require an interelectrode distance no greater than 0.5 mm. For an 8 X 8 array this provides for a 3.5 mm X 3.5 mm 'window' that neatly fits surgically onto the lateral aspect of the bulb. Contrary to intuition, spatial resolution depends not on the size of the individual electrode face but on the interelectrode distance. In order to minimize electrode noise the surface area of each electrode should be as large as possible. We have found that a stainless steel wire diameter of 0.25 mm is satisfactory (Eastman 1975). Further description is devoted to analysis of EEGs from these now standard arrays.

Fig. 22. Upper frame: the 1-D FFT is applied to the complex gain at the peak frequency of 20 bursts, and the means and standard errors (S.E.s) are plotted for normalized natural logarithm of power as a function of spatial frequency. For this subject the spatial frequency is asymptotic to a noise level above a spatial frequency of 0.4 c/mm. Lower frame: the 1-D FFT is applied to the complex gain of the simulated bursts in the same manner as for the EEG. The width of the focus is fixed at 1.0 mm; the 4 curves (2 of which are reproduced in the upper frames) result from shifting the origin of the granule cell Gaussian density function (Fig. 21) to points at the indicated distances from the circumference of the array along a radius of the array. A similar family is obtained by changing the width of the focus. The linear or circular array is useful for detecting high spatial frequencies (if any) with a limited member of electrodes, but it cannot be used to localize EEG foci or to measure their widths. (From Freeman and Baird 1986.)

The two-dimensional Fast Fourier Transform (2-D FFT) is applied to each 8 X 8 matrix after it is embedded in a 32 X 32 matrix of zeroes, multiplied by a Hamming window, and translated to bring the low spatial frequencies to the center of the output plane (Gonzalez and Wintz 1977). No significant departures from radial symmetry have been found in the two -dimensional spectra, so the output is radially integrated, divided by the power at 0 c/mm, and converted to the natural logarithm to give log normalized power as a function of spatial frequency in c/mm. If the phase pattern of the component is to be included as well, then we subject each of the 64 traces of EEG data to the 1-D FFT, calculate the real and imaginary parts of the complex amplitude at the desired temporal peak frequency, apply the spatial transform to the matrices of both parts, and sum the spatial spectral parts prior to derivation of the power spectrum, normalization, and logarithmic conversion.

Fig. 23 shows the characteristic shape of the spatial spectrum of the bulbar EEG in the form of the ensemble average and standard error over sets of 20 bursts from 6 rabbits. The energy decreases monotonically from near 0 c/mm with increasing spatial frequency to a plateau above 1 c/mm. The inflection occurs between 0.5 and 0.8 c/mm. A small peak commonly appears near I c/mm; this is an artifact due to the regular spacing of electrodes near 0.5 mm in the array. There is no other evidence for spatial periodicity. The rate of fall-off over the low spatial frequency range is determined by the size of an active EEG focus: a broad focus gives rapid fall-off, while a narrow or multiply peaked focus gives less rapid fall-off. This segment of the spatial spectrum invariably shows an upward convexity. This configuration proves that the energy of activity is indeed concentrated in the focus, as suggested by inspection of the EEG and of contour plots (Figs. 15 and 17). That is, the burst focus does not result from local synchronization of on-going activity with some spatially uniform temporal amplitude but with previously dispersed frequency and phase. Were that the case, this segment would conform to a straight line (Katznelson 1982). The burst and interburst segments do not manifest respectively 'synchronization' and 'desynchronization' in classical EEG parlance. This conclusion is confirmed by concomitant unit and EEG recording in the bulb. The mean rates of firing of single and multiple units wax and wane with respiration in concert with the amplitude of EEG burst activity; and spatially the local mean firing rate is proportional to the local EEG amplitude. The burst is a manifestation of augmented metabolic energy release, not of augmented coherence at a fixed energy level in local regions of activity comprising foci.

Fig. 23. The spatial spectrum of the EEG from an 8 X 8 array of electrodes on the bulb of a rabbit is shown as the normalized natural logarithm of power versus spatial frequency in c/mm. The mean and standard error are pooled for 20 bursts from each of 6 rabbits. The open dots show the spectrum of a uniform rectangular pattern subjected to the Hamming procedure; this defines the upper limit on spatial analysis that is imposed by the size of the array window. The 4 curves show the spectra predicted for pools of granule cells with activity density functions in the shape of a bivariate Gaussian distribution having standard deviation W in mm. The lower limit for this model is shown by the curve for width W = 0. (From Freeman and Baird 1986.)

In order to understand and explain the form of these spectra we construct spectra from a volume conductor model of the bulbar mechanism (Freeman 1975, pp. 212-228, 1978), with output that simulates the convex low-frequency segment of the EEG spatial spectrum. Given a pattern of potential field on the surface of a volume conductor, there is an infinity of possible source-sink distributions compatible with it (Nuñez 1981; Katznelson 1982; Ch. 13 and 14). If recordings are made at selected points throughout the volume, and if close studies are made of the cytoarchitecture and physiology of the neurons in the central region of the potential field, the possibilities become severely restricted. The approach used is to postulate a set of sources and sinks that is consistent with depth recordings and the properties of a putative neural generator, to calculate the surface potential at locations corresponding to the geometry of the recording array, and to adjust the parameters of the sources and sinks until the calculated potential optimally conforms to the observed potential. We use the electrostatic equivalent as our approach for the bulb (Freeman 1975, pp. 172-185), in which each local source or sink is represented by an equivalent point charge. This differs from the electric field approach as used by Rall and Shepherd (1968) for the bulb. The results are the same, but the scalar computation of the potential from Coulomb's law is much simpler than constructing a current network.

We postulate that the granule cells sustain the EMFs of the currents that underlie the dominant EEG oscillations, which are added with currents contributed by other cells to the raw EEG. The granule cell sources and sinks constitute a distributed current dipole in a plane parallel to the bulbar surface at an average depth of approximately 0.58 mm with point separation of ±0.17 mm (Fig. 21, lower frame). We represent the regional distribution of the dipole moment density under the surface with a bivariate Gaussian density function of standard deviation W. The potential function at the surface is bell-shaped; the width of the bell varies in near-proportion to the value of W. The rate of fall-off in spectral energy with spatial frequency varies inversely with width, and therefore with W. When W = 0, the width of the bell is minimal (for a dipole depth of 0.58 mm the half-amplitude radius at the surface is about 0.5 mm); we call this limiting case the 'point spread function' of the granule cell surface potential. Its spectrum defines one extreme for the rate of fall-off that can be expected in observing granule cell activity (Fig. 23, W = 0). The other extreme is defined by the Hamming window applied to a square with the designated size of the electrode array (the dotted curve in Fig. 23). Observed spectra lie between these limits with predicted radii between W = 1.0 and W = 1.5 mm. We confirm these values by measuring the widths of the corresponding EEG foci in the spatial coordinates.

The energy manifested in the plateau at higher spatial frequencies cannot be attributed to granule cells and should be removed irrespective of its various sources. We do this by using a low-pass exponential filter (Gonzalez and Wintz 1977) as shown in Fig. 24. The matrix of 64 amplitudes of a component from a burst is transformed by the 2-D FFT without Hamming; the filter is applied to the two-dimensional image in the. spatial frequency domain to attenuate the energy at high frequencies; the inverse 2-D FFT carries the filtered image back to the spatial domain. We estimate that this procedure removes 15-25% of the energy in bursts, and up to 40% in interburst segments, when the cut-off frequency of the filter is 0. 5 c/mm.

Fig. 24. The performance characteristics are shown for the exponential spatial low-pass frequency filter applied to spatial spectra as normalized logarithm of power derived by the 2-D FFT without Hamming. A(f) = exp [ln(2^-0.5) (f/f_o)ⁿ], where f is frequency in c/mm. The cut-off at f_o is the 3 dB fall-off point. From the results exemplified in Fig. 23 the cut-off has usually been set at f₀ = 0.5 c/mm and n = 4 (see also Fig. 29). For high pass filters n = -2. (From Freeman and Baird 1986.)

The final step is to apply spatial deconvolution in order to correct the image for the distortion caused by volume conduction (Freeman 1980). We approximate the bulbar generator by an array of 64 elements, one underlying each electrode. We represent the activity of each generator by the variable u(t) and the overlying potential by v(t). Then each potential is the sum of contributions from all 64 generators.

where the a_ij are coupling coefficients determined by the point spread function and the distance between each pair of electrodes. In matrix notation at each time t, or for each component over t,

Therefore

We multiply the 64 X 1 column vector V of amplitude values by the inverse of the 64 X 64 matrix A to recover an estimate of U, the granule cell activity density function. The procedure is repeated for each desired component in each burst of interest.

The procedure acts like a high pass filter; it accentuates contributions from current sources superficial to the granule cell generator (Martinez and Freeman 1984) that have narrower point spread functions, and also the electrode noise for which this function has zero width. Thus, low pass filtering must precede deconvolution. The activity from more distant current sources is dis-accentuated, and activity at the monopolar electrode that appears spatially uniform in amplitude is relegated to an edge effect.

The operation of deconvolution is in principle related to the process of focusing a blurred image taken at improper focal depth. For this reason we call the inverted matrix, which is a type of Laplacian operator, a 'software lens.' The optimal focal depth is predicted from the anatomical depth of the generator for an experimental evaluation (see Section 18.12, Fig. 29).

This set of techniques should find use in scalp recording of EEGs (Ch. 3 and 13). First, curvilinear or rectangular arrays of 64 electrodes at regular spacings (e.g., 1.0 cm, 7 cm X 7 cm) and the 2-D FFT at selected times (or of selected components measured by curve-fitting) will give the spatial spectrum that serves to specify the spacing and window size that is appropriate for different experiments. Second, spatial filtering may serve to minimize EMG potentials. For this purpose it will be desirable to determine their point spread functions, propagational patterns, and spatial spectra, so that appropriate spatial filters can be designed to separate them from the EEG. A low pass filter would suffice only to the extent that the EMGs could be shown to have narrow point spread functions owing to the locations of their generators just under the scalp. Spatial averages should be taken along directions parallel to the muscle fibers. Third, the contributions to 'monopolar' recordings from sources affecting the reference but not the array electrodes (whether from ECG, EEG, EMG, or eye movements), may be attenuated for the inner 6 X 6 or 4 X 4 subset of the 8 X 8 array of electrodes by use of spatial deconvolution of some alternative Laplacian operator, depending on the spatial properties of the generator of the desired components of the EEG.

18.12. Optimization of measurement; the behavioral touchstone

The procedures described in the last two sections can give access to the detailed spatial patterns of the EEG of defined neural populations in local regions of brain tissue. However, this is new and largely unexplored terrain. We have no clear bases on which to form expectations of what the patterns ought to look like, no equations to solve to provide archetypal shapes, and no analytic spatial basis functions to fit to the data (other than the journeyman bivariate gaussian). We can print out large numbers of these patterns as contour or density plots (e.g., Figs. 17 and 20) and tape them to the walls of our laboratories for inspection, but they quickly exhaust even the remarkable powers of the human visual system for detecting regularities of form. Our task now is to discuss procedures for condensing these vast quantities of data, for seeking regularities in them, and for assessing the significance of differences between subgroups of samples of EEG spatial patterns. The current availability of large arrays of electrodes and of pre-amplifiers, and the recent introduction of optical techniques for observing color changes in voltage-dependent dyes applied to neural tissues indicate that this task will be approached by increasing numbers of investigators in the near future. At present the only body of local EEG data (as distinct from scalp EEG data) that is sufficiently large to support analysis is that from the olfactory system of the rabbit. Therefore, this section consists largely of an example of spatial analysis for the rabbit olfactory bulb. Here we follow, and thereby illustrate, the outline laid down in the first several sections of this chapter.

We choose our basic model or hypothesis (Section 18.1) on the premise that the olfactory bulb is involved in the neural mechanism for olfactory discrimination. The receptors form a two-dimensional array in the nose, and their axons project to the bulb with a degree of topographic order in the primary olfactory nerve. Each receptor responds selectively to some odors and not to others. We assume that each stimulation by an odor during an inhalation causes a unique spatial pattern of receptor activity for that odor and inhalation, resulting in another unique pattern of excitation onto the bulbar surface. If an animal is able behaviorally to discriminate the odor, we infer that a spatiotemporal pattern of activity representing that odor must exist in the bulb and be transmitted to the prepyriform cortex as the basis for discrimination, whether or not the discriminatory process is brought to completion in the bulb. Further, there must be a different neural activity pattern for each discriminable odor, and some attribute of each pattern must be invariant over every occurrence in which correct discrimination of its odor takes place. Our hypothesis is that these activity patterns are manifested in some aspects of the bulbar EEG. We predict that when the EEG bursts have been properly measured, the numbers will be found to contain odor-specific information.

This hypothesis is also applicable to studies of primary sensory and motor neocortex. It is obviously more problematic in other areas, suggesting that the search for basic EEG mechanisms might best be focused at present onto the main sites of entry and exit of the forebrain (Ch. 17 and 19).

The key requirement for testing this hypothesis is a database (Section 18.2) that is recorded from the brain during performance of odor discrimination. We know that each animal generates its own generic spatial pattern of bulbar EEG activity, and that this pattern changes to a new form whenever the animal is trained to respond to a test odor (Freeman and Schneider 1982; Viana Di Prisco and Freeman 1985; Gray et al. 1986). We infer that odor-specific information may be manifested in the EEG only after the completion of training.

We condition thirsty rabbits to respond with licking (with or without sniffing) to one odor CS+ by presenting it paired with delivery of water, and to respond with sniffing and not licking to a second odor CS- by presenting it on randomly interspersed trials without water. Four sets each of 90 bursts accumulate over 3 sessions: odor bursts CS+ and CS-, and control bursts C+ and C- taken just prior to stimulation on each trial. The aim of analysis is the development of techniques that optimally separate bursts CS+ and CS- into tight clusters distinct from each other and from clusters of bursts C+ and C-, but that do not segregate the control bursts C+ from C-. We restrict consideration to those trials in which 'correct' responses occur, and to the 4 of the 5 rabbits that learn the task.

Because we lack analytic basis functions (Section 18.3), we find the rms amplitude for the EEG trace on each channel, and treat each burst as a 1 X 64 vector. In order to obtain an overview on large sets of bursts (e.g., 120) we use the non-linear mapping algorithm of Sammon (1969), which maps each 64-dimensional burst pattern into a point in the two dimensions that optimally represent the degrees of similarity among a set of bursts by the Euclidean distances between their points. Examples are shown in Fig. 25. The upper frame shows the points for 60 control bursts from two rabbits. The subject differences are reflected by the two clusters, and the background variation is reflected by the scatter in each cluster. The lower frame shows the scatter of 120 control bursts from a rabbit with a bimodal spatial rms EEG pattern. That is, there are two maxima of EEG burst amplitude P₁ and P₂ separated within the array window by about 2 trim, and the dominance of one over the other varies randomly on sequential bursts. These are labeled by open and solid dots; the bursts with highest and lowest mean rms amplitude are labeled respectively by (+) and (-).

Fig. 25. Upper frame: the results are shown of applying the non-linear mapping algorithm of Sammon (1969) to 60 control bursts from 2 rabbits, showing that the differences between subjects are easily detected by this procedure, and that the cluster for each subject is unimodal. Each burst is first expressed as a 1 X 64 vector of channel rms amplitudes after bandpass filtering (10 and 150 Hz) and then as a point in 2 selected dimensions for display of the degree of similarity to other bursts by a Euclidean distance between points. Lower frame: a set of 120 control bursts from another subject likewise reveals a unimodal distribution; see text for the meaning of the plotting symbols. (From Freeman and Baird 1986.)

Each of these features makes some contribution to the discrimination among bursts. The outstanding result is that the distributions of control bursts (and of odor bursts as well) are unimodal with a density that tends toward the gaussian. Then each cluster can be described by a centroid and a pooled standard error (S.E.) in the two dimensions (Ch. 16 and 17). We apply cluster analysis to sets of bursts including both the controls (C+ and C-) and the test odor bursts (CS+ and CS-), finding (Fig. 26) that the odor bursts differ in rms spatial pattern from the control bursts but not from each other. The standard errors (S.E.s) of the odor burst clusters exceed the distances between their two centroids. Moreover, the S.E.s of odor clusters exceed the S.E.s of the control clusters (Viana DI Prisco and Freeman 1985).

Fig. 26. The results are shown of applying the non-linear mapping algorithm to a set of 120 bursts from the olfactory bulb of a trained rabbit, Each burst is expressed as a set of 64 rms amplitudes after bandpass filtering the time series (10-150 Hz). The means and standard errors (S.E.s) are shown for the two sets of odor bursts (CS+ and CS-) and for the preceding control bursts (C+ and C-) on trials in which correct behavioral responses occurred. Each cluster is approximately bivariate gaussian. The t test using the pooled variance for each pair of clusters indicates that the centroids CS+ and CS differ significantly (P < 0.01) from the centroids C+ and C- but not from each other. Therefore the channel rms amplitude of bursts does not serve as an adequate descriptor for purposes of detecting, odor-specific information (cf., Fig. 30). (From Viana DI Prisco and Freeman 1985.)

The critical problem before us is to find methods for data analysis that separate the two centroids CS+ and CS- while not separating those for C+ and C- (Ch. 17). The most direct method is to calculate the centroid for each subset of bursts in 64-space (the means of the 64 channels over each control and odor subset after normalization by channel over the entire set), and to express each burst by its Euclidean distance from each of two centroids: C+ and C- for control bursts, and CS+ and CS- for odor bursts. The degree of clustering is expressed by the segregation of the points representing bursts as being closer either to its own centroid ('correct') or to an apposing centroid ('incorrect'), that is, to either side of the diagonal in Fig. 27. The effectiveness of a given analysis is expressed by the % correct classification of the odor bursts of both types less the % correct classification of the control bursts of both control groups. This percentage difference is our 'behavioral touchstone' for methodological optimization and evaluation (Section 18.6).

Fig. 27. The centroids of 64-space are found for the dominant components of 4 groups of bursts: C+, C-, CS+ and CS-. The Euclidean distances for each control burst are found from C- (ordinate) and C+ (abscissa, upper frame), and the Euclidean distances for each test burst are found from CS+ and CS- (lower frame). The effectiveness of burst separation by this algorithm is expressed by the percentage of bursts that lie closer to their own (i.e., 'correct') centroid than to the other. The behavioral touchstone is given by the difference between % 'correct' for odor bursts (86%) minus the % 'correct' for the control bursts (55%), in this example 31%. The bursts having peak frequencies f₁ < 55 Hz are omitted. The smaller spread of points for control burst reflects the characteristic tendency of bulbar EEG amplitude to decrease when an odor is presented to an animal. (From Freeman and Viana Di Prisco 1986.)

As predicted by the result 18 in Fig. 26 the description of bursts by their rms values does not segregate the two types of odor bursts significantly; the variance within the two clusters is too great. We use curve-fitting (Section 18.4) by the methods in Section 18.10, and spatial filtering and deconvolution (Section 18.11). Analysis of the amplitudes of the dominant components gives some reduction in the variance, but it is still insufficient to separate the odor bursts. The crucial step emerges from temporal spectral analysis (Section18.5). We find that the frequencies of the dominant components are broadly distributed across the range from 10 to 125 Hz with concentration in the range from 60 to 70 Hz, but there are substantial numbers of bursts with dominant frequencies between 20 and 60 Hz (Fig. 28, upper frame). The frequency modulation parameter FM in Eq. 17 of Section 18.10 on the average is negative (slowing of frequency during a burst) for frequencies greater than 60 Hz as predicted by the describing functions in Fig. 10, but on the average it is positive (increasing within-burst frequency) for frequencies lower than 60 Hz (Fig. 28, lower frame). Moreover, the variance of FM across bursts increases exponentially with decreasing mean burst frequency. These low frequency bursts are not accounted for or explained by the KII model as presently constituted (Section 18.9 and Appendix C). We postulate that these lower frequency bursts constitute a distinct class, that their irregularity in the temporal domain may be accompanied by large pattern variability in the spatial domain, and that they should be removed from the cluster analysis.

Fig. 28. Upper frame: the distribution is shown by the bars for the frequencies of the dominant components from 4 rabbits on trials in which a correct behavioral response occurred. The mean and S.D. are shown for each subgroup of the fraction of total burst energy incorporated into the dominant component. By this measure the coherence of bursts is greatest in the frequency of range of 50-80 Hz. The bursts with frequencies f₁ < 55 Hz during odor presentation outnumber those during control periods by 4 to 1. Lower frame: the frequency parameter FM in Eq. 1 of Section 18.10 is expressed as the % modulation of the mean frequency. A negative value means that frequency is slowing down within the burst, which is predicted during burst onset by Fig. 10; a positive value means that frequency is augmenting as during burst termination. A small part (<5%) of the variation in FM expressed here as the standard deviation about the mean (dots) may be attributed to sampling of bursts during editing. The large variance in FM at lower mean frequencies is attributed to chaos (Section 18.13). (From Freeman and Viana Di Prisco 1986.)

Fig. 29 (upper left) shows the effect on the behavioral touchstone derived from the analysis of dominant components of deleting bursts having frequencies lower than the value indicated on the abscissa. The optimal separation of bursts occurs with the cut-off at 55 Hz. An example of the Euclidean distances from one subject is shown in Fig. 27. The results of non-linear mapping applied to the reduced set of dominant components is shown in Fig. 30. The reduction of S.E.s for clusters of odor bursts suffices to show clear separation of their centroids. Comparison with Fig. 26 shows that the detection of odor-specific information in sets of higher frequency bursts is obscured by the large spatial variance contributed by the lower frequency bursts, when the latter are included in the analysis.

Fig. 29. Four examples are shown of use of the 'behavioral touchstone' for methodological optimization. Upper left: the stippled bars show the % difference (Fig. 27) as bursts with peak frequencies less than the indicated frequency are deleted from the cluster analysis. The peak occurs because the separation of odor bursts is high with the cut-off at 55 Hz, but equally important, the separation of control bursts is at a minimum. This is 'correct' for the control bursts, because the control period is the same preceding presentation of both CS+ and CS- odors (see also the legend for Fig. 30). The proportion remaining after deletion is shown by the open bars. This result supports the criterion of peak frequency f₁ < 55 Hz to identify bursts as chaotic. Also deleted are bursts with -50% < FM < +50% and bursts for which the fraction of burst energy in the dominant component is less than 0. 15, comprising together about 10% of bursts. (From Freeman and Grajski 1986.) Lower left: the % difference is used to evaluate the effects of n-th order exponential spatial filters on the EEG as in Fig. 24 (Gonzalez and Wintz 1977). The upper curve shows that maximal separation of odor bursts occurs with a 4th order low-pass filter cut-off frequency (down 3 dB at 0.5 c/mm), as predicted by spatial spectral analysis of the bulbar EEG. The lower curve shows that minimal separation occurs with a 2nd order high-pass filter set at 0. 17 c/mm, as predicted by the average of measurements of the widths of active EEG foci, which is close to the passband predicted by the width of the electrode array. The order parameter (N) is optimized in the same way. These results prove that the odor-specific information exists in the spatial passband that is characteristic of the bulbar granule cell EEG generator (Freeman 1975, Ch. 4, 1978, 1980). (From Freeman and Baird 1986.) Upper right: the % difference is used to determine optimal depth of focus for the 'software lens' (Section 18.11) to correct for the effects of volume conduction on the bulbar EEG. The theory implies that optimal focal depth is at the mean dipole depth of the granule cell dipole generator (Fig. 21, lower frame). Penetrating microelectrodes show that the site of reversal of ERP and EEG potential for active foci restricted to the near-planar lateral wall of the bulb lies on the average 0.1 mm. superficial to the mitral cell layer (Freeman 1975, pp. 219-228). The mean depth of the layer on postmortem histological examination in these rabbits is 0.71 ± 0.18 mm; the optimal depth of focus is 0.49 or 0.54 mm for the rabbits giving this data set, which is about 0.1 mm less deep than predicted but within the large standard error. (From Freeman and Baird 1986.) Lower right: the % difference is used to evaluate the optimal number of factors to be used in the Euclidean distance algorithm after factor analysis of the data. Two examples are shown from the 4 rabbits, one with maximal separation using 4 factors, the other with 9 factors. When all factors with non-zero eigenvalues are used, the subgroups all approach 100% correct classification, and the % difference measure goes to zero. (From Freeman and Grajski 1986.)

This point is confirmed by inspection of the correlation matrix of the dominant component amplitudes for each subject. The mean pair-wise correlation coefficient (by z-transform) for control bursts with frequencies f₁ > 55 Hz exceeds 0.9 within subjects; that for odor bursts exceeds 0.8; but that for bursts with frequencies f₁ < 55 Hz lies between 0.1 and 0.2. This is remarkable, in that spatial ensemble averages over bursts with f₁ > 55 Hz and over bursts with f₁, < 55 Hz both conform to the same spatial pattern characteristic for each subject. In other words, the spatial patterns of the one group are predictable from odor condition, whereas the spatial patterns of the other group are not, although both groups must lie within certain bounds, in order that their ensemble averages conform to the same norm.

This finding is important for the classification of bursts by the criteria discussed in Section 18.9; we take this up in Section 18.13 and Appendix D. Here we conclude that, as in Section 18.4 we distinguish between 'correlated' and 'uncorrelated' variance in order to remove the latter from sets of ERPs, we must likewise detect and remove the massive variance contributed by an as yet unexplained type of EEG event, before proceeding to cluster analysis and factor analysis of well-behaved EEG events. This point is obvious in respect to known sources of artifact such as the EOG in scalp-recorded ERPs. Here we deal with more subtle aberrancies, the possible nature of which we discuss in Section 18.14.

This completes one pass through the measurement of the spatiotemporal properties of the bulbar EEG. The next step is to repeat while optimizing each step by reference to the touchstone. Two examples are shown in Fig. 29 (lower left) for optimization of the depth of focus of the 'software lens' (Section 18.11), and for optimization of spatial filters (lower right) such as that illustrated in Fig. 24. These optimized data provide a basis for evaluating a confidence interval for values of the touchstone. This. is done by 'jack-knifing' (Ch. 17): the procedure is repeated N times for each subject with N bursts, while deleting a different burst on each computation. The values for each subject are normally distributed. The square root of the mean within-subject variance is 4.09, so that the 99% confidence interval for the touchstone is ± 10%. By this criterion the touchstone is above zero for each of the subjects. However, the touchstone is not a conclusive statistical test; it is a rapid and inexpensive assay of the procedures used to optimally restructure data in order to test hypotheses. The first priority in correlation of the EEG with behavior is not to explain the brain; rather it is to validate the measurements before introducing the data into statistical and explanatory models.

18.13. Spatial basis functions from factor analysis

These results of cluster analysis suggest that, during presentation of an odor as a CS to a trained animal, the spatial pattern of the dominant oscillatory component of the EEG tends to converge to a stereotypic form, such that the underlying neural activity mediates the response to the presence of that odor (Freeman and Skarda 1985). Our behavioral touchstone indicates that convergence either is abrupt but often delayed, or is gradual. Over the 3 successive bursts recorded about 0.5 sec apart on each trial after odor arrival the touchstone value rises from 8.9% on the first sample (65.2% odor -56.3% control correct classification) to 16.1% (74.1% - 58.0%) on the second sample and to 26.6% (82.4% - 55.8%) on the third sample. The different spatial patterns are established over most or all of the array window. This is shown by the fact that when the touchstone is repeatedly recalculated after deleting groups of 8 channels, beginning with those 8 having the lowest amplitude or vice versa, the goodness of separation deteriorates progressively. This conclusion is supported also by detection of widespread phase gradients such as that exemplified in Fig. 17. Our requirement stated in Section 18.12 for an invariant pattern over repeated presentations of the same odor, but a different pattern for different odors, may be met by the spatial pattern of the amplitude of the oscillatory burst over the entire array, and, by inference from intrabulbar recording (Bressler 1984) over the entire bulb. Accordingly the proper basis function for spatial analysis must be a pattern that extends over the entire array, just as the proper basis function for temporal analysis extends over the entire duration of an ERP (Section 18.4).

The spatial basis functions are empirical; they may be derived from the data by selective averaging, or principal components analysis (PCA - Ch. 5) over a set of dominant components from each subject and stage of training. Several considerations enter into the selection of a procedure, mainly involving the required conditions of linearity, normality and equality of variances. First, the physiological model to be tested is that the spatial EEG pattern tends to converge to one form during the control state and to different forms during each odor state. However, at best it appears that less than 80% of recorded bursts approach convergence in the odor states, since the degree of convergence captured in the time window of sampling may vary, and errors of measurement and classification obscure underlying regularities. Apart from removal of bursts labeled as 'incoherent,' there is no a priori test as yet to determine which recorded bursts optimally represent convergence. Therefore, selective averaging is not feasible.

Second, a number of findings all indicate that the spatial patterns are not restricted to high-amplitude foci or to the area of the array window; probably they cover the entire bulb. These findings include: the relatively equal contributions of all channels to burst separation by the touchstone, the common instantaneous frequency, and the common phase gradient approximating a planar or conic surface over the array. There is evidence from anatomical analysis that the EEG foci occur over bulbar sites that have higher densities of mitral cells and glomeruli, perhaps corresponding to 'olfactory foveas' (Meisami and Emamian 1985). This means that the high EEG amplitudes in the focus characteristic of most rabbits have no significance with respect to odor discrimination apart from the degree of acuity; they are accidents from the anatomy. Then the procedure we adopt of normalization by channel across sets of bursts, which we justify on the basis of optimal burst separation according to the touchstone, can also be justified in physiological terms, that each local part of the bulb contributes to discrimination irrespective of the mean local EEG amplitude, the latter being due to anatomical' factors unrelated to specific odors (Freeman and Van Dijk 1987).

Third, in all odor subsets the variations in burst amplitude from the means greatly exceed the differences between subset means. The amplitudes of control bursts generally exceed those of odor bursts (Fig. 27). Although this aspect of bursts does not enter into CS+ from CS- odor discrimination, the groups should be normalized to equal variance prior to PCA. The magnitudes of the variances are approximately equal across sets of odor bursts after removal of the 'incoherent' bursts. Finally, the distributions of channel amplitudes typically are skewed owing to the relatively small numbers of very high values contributed by channels at or near the peaks of foci, but the distributions are normal after the normalization by channel across entire sets. For these several reasons we use PCA of the correlation matrix rather than of the variance-covariance matrix (Johnson and Wichern 1982).

The statistical model is

where u_1,i,j is the amplitude of the normalized and filtered dominant component (Eq. 17 in Section 18. 10 and Eq. 18 in Section 18. 11) on the i-th channel (i = 1, 64) of the j-th burst (j = 1, J bursts entered into the analysis); X_i,k is the i-th loading of the k-th factor derived from the eigenvector matrix for k_ee = 1, K factors; a_j,k is the k-th factor score for the j-th bursts; and e_i,j is an error term assuming independent 'noise' over channels and bursts. In matrix notation

where U, is dimensioned (64 X J), A is (J x K), X is (64 X K), and by assumption the residual E is (64 x J). The latter assumption must be tested by examining the residuals.

PCA consists simply in rotation and translation of the data set to new coordinate axes; these are defined by the largest few components of the variance that are built into the data set by odor testing. The main difficulty is that these axes are very sensitive to distortion by outliers resulting from mistakes in editing and preliminary classification, by errors of measurement, or by use of overly small samples. Different subsets of data might give different factor loadings X as estimates of what we think of as stereotypic underlying neural activity patterns; what we seek is factorial invariance. The approach we use is to structure the analysis so as to optimize the number of factors K. Fig. 29 (upper right) shows use of the behavioral touchstone for this purpose on the factor scores instead of the channel amplitudes. Two examples are shown of the effect on the % difference of increasing K from 2 to 12. Two rabbits have maxima at K = 4 and the other two at -K = 8 and K = 9. There is no significant change in the % difference with the use of factor scores at these maxima from the use of channel amplitudes. The touchstone serves also to evaluate use of Varimax rotation. Factors are selected according to the Kaiser-Gutman rule (those with eigenvalues these incorporate on the average 92% of the variance in 4-10 factors. After rotation the touchstone is again used to test the factor scores. The maximum % difference is at K = 5 for two rabbits, at K = 4 for another, and at K = 3 for the last, with mean % difference is equal to 24.8% (78.4% - 53.6%). Oblique rotation deteriorates separation by the touchstone, suggesting that the factors are orthogonal. Expanding the data set to include bursts from all trials instead of only those trials on which a correct response occurred results in K = 4 for all 4 rabbits, but the touchstone is reduced from 24.8% to 17.9% (70.5% - 52.6%). In any case the number of variables used in the separation of bursts is less than the number of bursts in every subgroup to be separated (Freeman and Viana Di Prisco 1986), so that indeterminacy is avoided in equations used to test statistical models.

An important check on the procedure is to reconstruct the patterns in accordance with Eq. 21, and by correlation and comparison of sums of squares, to verify that the reconstituted patterns conform to the data patterns less the residuals. The factor scores must be computed by the regression method (Johnson and Wichern 1982); use of factor scores by the ordinary (unweighted) least squares method is not satisfactory. The sums of squares of the differences between the data (normalized channel amplitudes) and the reconstructed channel values from the sum of the product factor scores and factor loadings in Eq. 21 are divided by the sums of squares of the data. In the example under discussion the average value for the correctly classified bursts for the 4 subjects responding behaviorally was 0.039, and for the incorrectly classified bursts 0.136. For the subject without -behavioral evidence for discrimination the ratio was 13.77, showing that the results of factor analysis were meaningless.

The next step consists in testing for factorial invariance. We note in Section 18.6 that success in fitting a sum of basis functions to an AERP does not suffice to define an appropriate set; the same basis functions must hold over sets of ERPs from differing behavioral states, and the coefficients must be well behaved. Analogously, the factor patterns (spatial basis functions) must be shown to hold over independent subsets each with its own covariance matrix. The data set is divided into a 'learning' subset used to define the factors and a 'validation' subset used to test their efficacy in classification. The multiple. factor patterns are compared by correlation (Table 1). Optimally we must show that one set of spatial basis functions derived from a 'learning' subset from an animal within a behavioral state, such as responding to a certain set of odors, holds for all subsets over the duration of that state and does not hold for subsets taken before or after that state.

PCA, like spectral analysis, merely puts the data into a new and reduced coordinate space; numerous procedural safeguards are needed. These statistical procedures, together with various forms of discriminant analysis of factor scores derived by any of half a dozen methods, go beyond the scope of this chapter and are discussed in Ch. 5, 16 and 17. Factor components, like spectral components, are not ends in themselves; they are the raw materials of model-building. This aspect is taken up in the next section.

TABLE I: RESULTS OF USING CORRELATION OVER FACTOR LOADINGS DERIVED FROM ANALYSIS OF SUBSETS OF THE DATA TO TEST FOR FACTORIAL INVARIANCE

In A the data are divided into 2 subsets by partitioning the entire data set into 2 parts, one of odd numbered bursts and one of even-numbered bursts. In B the factor analysis is done on the data from each of the consecutive 3 behavioral sessions yielding the bursts. The asterisks denote the correlation coefficients of the first 4 factors reversed in sign when necessary, and (re-ranked in the few instances where the order by rank of the eigenvalues differed over sessions) averaged (by Fisher's z-transform) over subjects and also (in B) inter-session comparisons. These high coefficients reflect the stability of the factor patterns over sessions. As might be expected, each set of factor loadings from a 'learning' subset serves approximately as well to classify bursts from a 'validation' subset. (From Freeman 1985c and Grajski 1986.)

18.14. Limit cycle and chaotic attractors manifested in the EEG

These optimized temporal and spatial measurements give us the essential findings that we need in order to choose criteria for characterizing the non-linear modes of the EEG. The background activity, best seen in resting but awake, unmotivated animals without bursts, consists of long stretches of irregular fluctuations that are unpredictable over segments longer than a few hundred milliseconds and unreproducible from one segment to the next. So also are their spatial and spectral patterns. Ensemble averages of successive EEG time segments tend to approach a straight line because of EEG phase fluctuations, which is the basis for extracting ERPs that are phase locked to an event. Yet the activity is clearly bounded in its distributions of amplitude and frequency, because ensemble averages of spatial and spectral patterns converge to norms characteristic for each subject. Most importantly, no matter how erratic the wave-form of any time segment may appear, the same time series occurs concomitantly over the entire array, although with unpredictable local differences in EEG amplitude. This means that either each sequence emerges from within the bulb by cooperative interaction, or that the sequence is imposed from outside. The olfactory receptors have no known mechanism by which to generate such patterns, nor is there any known source of centrifugal input with the required high-pass transmission line characteristic. In any case when centrifugal input is surgically blocked, the EEG continues with little change in its time-space characteristics, mainly in the loss of flexibility in relation to behavior (Freeman 1975, pp. 407-414).

The bulbar EEG appears to be random (or chaotic) by several criteria: (a) the EEG amplitude histogram from segments one or more minutes in length closely approximates the gaussian density function, though with truncated tails; (b) the form of EEG fluctuations resembles 'shot noise' generated by passing white noise through a bandpass filter; and (c) the envelope of the EEG conforms to the Rayleigh distribution (Freeman 1975, p. 148). However, the only known source of apparently random action within the bulb is in pulse trains of individual neurons; interval histograms and autocorrelation functions show the form of a Poisson process with a dead time (Freeman 1975, pp. 150-154). That is, following each pulse there is a refractory period with zero pulse probability, but thereafter the time of occurrence of the next pulse cannot be predicted from preceding pulses. The cross-correlations of pulse trains from pairs of neurons show no evidence for covariation apart from their participation in cooperative mass action. No mechanism has been quantitatively proven whereby statistical regularities might cause an equivalent to brownian motion to appear locally in the field potentials of clusters of granule cells. Even if there were such a mechanism, it could not account for the widespread spatial coherence of the temporally unpredictable basal activity. We infer that the basal EEG activity may manifest the existence for the bulbar mechanism of a chaotic attractor, for which the output is a pseudo-random process, and that the non-zero equilibrium basin of the KII set is inaccessible to the in vivo KII mass. Our criteria for this conjecture are the boundedness in temporal, spatial, and spectral distributions, the unpredictability of temporal wave-form and spatial amplitude pattern within those bounds, and most importantly the widespread spatial coherence manifested in the common instantaneous frequency over the surface. Appendix D gives an example of the quantitative treatment of this distinction between noise and chaos.

The EEG burst with inhalation clearly manifests augmented metabolic energy release with the emergence of order from chaos. It discloses a bifurcation to a different basin of a different attractor, the result of the input and the attendant non-linear parametric changes in the bulb (Freeman 1979b). The EEG tends toward sinusoidal oscillation, the spectrum becomes more sharply peaked, and at high temporal frequencies the spatial pattern of amplitude tends to converge to a reproducible pattern. By these criteria we propose that this burst activity represents convergence toward a limit cycle (Appendix D), and that the convergence aborts after the termination of the input. We speculate that a limit cycle attractor is formed during learning to respond to each distinguishable odor given as a CS (including the background odor complexes). These criteria appear to hold only when the frequency of burst oscillation is above a certain level or range at or near 55 Hz (Fig. 30). Clearly for bursts with lower frequencies there is spatial incoherence and temporal incoherence as well, manifested in the high values for the FM parameter in the curves fitted to their dominant components and in greater spread of burst energy across components (Fig. 28). We believe that these events manifest the existence in the bulb of one or more chaotic attractors other than that manifested in the interburst periods and at rest (Freeman 1986, 1987).

Fig. 30. The map shown is derived by use of the non-linear mapping algorithm of Sammon (1969) applied to bursts having frequencies f, greater than 55 Hz (compare with Fig. 26). The means and standard errors are shown for 4 sets of bursts with the indicated numbers of bursts in each set: 2 sets of odor bursts (CS+ and CS-) and 2 sets of bursts taken in the control periods immediately preceding odor presentation (C+ and C-). Because the odors are presented on randomly interspersed trials in the same session, the control bursts are considered to come from the same state and should not differ significantly. (From Freeman and Viana Di Prisco 1986.)

Now we consider how chaotic activity can arise from within a KII mass such as the bulb. Inspection of each frame in Fig. 10 will show a diagonal line to the left of the vertical axis starting near w = 400 radians/sec. This is the approximate location of the boundary of the stable domain below and an unstable domain above. The instability arises in the KII linearized equations, because, with certain combinations of the feedback gain coefficients that place a complex pole above the diagonal line, another pole moves on the real axis to the right of the origin. In the amplitude range of the normal EEG we infer that this pole is due to mutual excitation in the component KI, set. The non-linear model is stable, because with increasing output amplitude (and therefore feedback gain), the pole moves to the left toward and beyond the origin. This analysis indicates that the source of bulbar chaotic activity lies in the tendency for the KI_m mass (the mitral cells) to generate runaway excitatory activity that is curbed (with an intrinsic delay) by the recurrent activity of the KI_G mass (the granule cells). This curbing effect occurs only after, and to the extent that, excitation has- already grown. On each recurrent pass through the negative feedback loop, the gain and therefore the frequency is different. That is, in 2 degrees of freedom afforded by feedback interaction in the KII mass, activity is divergent in respect to mutual excitation, but convergent in respect to negative feedback. Then an external 'random generator' need not to be postulated to account for the bulbar background activity; it manifests a chaotic attractor of the KII_MG mass. Its behavioral significance is discussed in Section 18.14.

We conclude that a stable equilibrium point is not accessible in the state space of the bulbar mechanism, other than that which becomes apparent under deep anesthesia. That is, the background EEG does not reveal an equilibrium state under external random perturbation. Most important at this stage is that the concept of chaos gives us a basis for partitioning the variance among sets of bursts into 'correlated' (pattern differences between odors) and 'uncorrelated' (chaotic bursts) fractions. In physiological terms it is reasonable for us to postulate that over repeated inhalations the bulbar mechanism may often fail to converge to a coherent spatiotemporal activity pattern. Thus, it makes sense for us to look for criteria by which will allow us to label such failures as members of an aberrant class or set of classes. This step will be crucial in all future analyses of EEG spatiotemporal patterns.

18.15. Conclusion. The KII set and perception

The key to understanding is measurement, and the key to measurement is understanding. We raise ourselves from level to level of explanation only by building on the foundation of prior constructions. We began this discussion with waves, and will end it here with attractors, but we have explained only a small sample of the EEG. We have before us a vast area to explore. It is a terrain of specialized and adaptive generators, some of which contribute directly to the EEG by accidents of the geometries of their electric fields, others clandestinely by their interactions with those we can observe. They are all interconnected within and across different levels of a hierarchy. Some groups of neurons exchange inputs and outputs, and some serve by their outputs to vary the parameters of others so as to induce bifurcations. Overall, the nervous system displays deterministic evolution that is guided by its own internal development structure and is modulated by the environment. We now have the conceptual tools that in principle are adequate for understanding the basics of this complex deterministic machine, but to use these tools effectively we must first employ the classical tools of linear systems and multivariate statistical analysis. These standard tools of physics and engineering provide our bedrock of comprehension and orientation. To neglect them is to trust in luck and dead reckoning, and to risk wandering aimlessly in the endless deserts of the EEG while hoping always for a waterhole.

This essay is concerned mainly with procedures used at the interface between raw data and hypotheses, and therefore with measurement procedures and the validation of those procedures by use of behavioral touchstones. The examples have been selected and simplified in order to illustrate the main outline of the procedures: pre-structuring of the experimental conditions; organization of the data base, selection of basis functions, measurement by curve-fitting, synthesis of a flow diagram; formulation of the dynamic equations, and interpretation of the coefficients. 'Emphasis has been given to the use of electric stimulation in lumped, linear domains, and to ordinary linear differential equations, because this is an optimal place of embarkation for further studies. Sensory stimulation is feasible but requires an understanding of the transformations occurring between the stimulus and the activity arriving at the locus of interest. Analysis of temporal dispersion, propagation delay, diffusion, spatial divergence, and distributed feed back require partial and integrodifferential equations that fall within the same procedures but are unsuitable for an introduction (Freeman 1975, Ch. 2 and 4). Various types of non-linearity require their own special techniques for analysis, but many of them can be linearized.

Modeling the EEG involves greater complexity. The 'inputs' may appear as random variables, and changes in wave shapes may derive from state changes described as bifurcations in large-scale dynamic systems. Yet the principles are the same. Any analytic basis function for which the derivative is calculable can be used for measurement. If it is inferred from the' data, it can lead to a generic differential equation; if it results from a model, it may lead to new precision of measurement and insight of interpretation. In either case it can lead to better understanding of brain function and behavior.

We have not yet completed the systematic analysis of chaos and of the bifurcation between chaos and limit cycle in the KIII set. When we first encountered chaos while developing the KII model, we were dismayed because unprepared to cope with it, and we firmly shut the door on it by imposing the necessary constraints on the design. Now mathematicians have provided biologists and electroencephalographers with the conceptual tools needed to cope with chaos (Helleman 1980; Abraham and Shaw 1983; Garfinkel 1983), and we should exploit them. In retrospect it is all too obvious that the EEG manifests the activities of chaotic and limit cycle generators. We should attempt to interest mathematicians and physicists in our work. They will need our help in the design of realistic and testable models; we will need their help to set up the state spaces of the parameters in Lagrangian coordinates, to establish the requisite hyperbolic metrics, to measure the fractal dimension of chaos (Appendix D) and to define the. trajectories on which bifurcations occur. The KII set provides an excellent starting point for understanding the cellular biophysics of the EEG. It is, after all, an ultimately simple, deterministic, generalized descriptor of the dynamics of a mass of non-linear excitatory and inhibitory neurons. We predict that it will become a prototypic and essential module for all parts of the cerebral cortex (Rössler 1983).

The KII set has equivalent value as a starting point for psychological analysis of the brain states and operations that underlie goal-directed behavior. These include expectancy, attention, and the capabilities of animals to establish and use 'search images' (Freeman 1983; Freeman and Skarda 1985) that constitute the selective neural control of input based on past experience. Over the past 3 decades consensus has grown that the element of cognitive function is not the single neuron but the nerve cell assembly (Hebb 1949). Von Neumann adopted this view in his last work (1958). Turing (1952), having dropped his interest in logical machines, similarly turned his remarkable intelligence to the study of the morphogenetic properties of large systems of interactive particles, thereby establishing a new field of mathematical analysis of living systems. Numerous works on the theory of self-organizing ensembles have since appeared in the past 5 years.

If the nerve cell assembly is the basis for cognitive function, then it is essential that we have first the tools to make direct observations on the spatial activity patterns of assemblies with arrays of electrodes for simultaneous recording from waking animals engaged in cognitive behavior, and second that we have models that explain how the observed patterns are generated and how they help to determine behavior.

Our present interpretation is that the KII nerve cell assembly in the bulb can sustain an indefinite number of stable attractors. We suggest that in the background state it has a chaotic attractor, which holds it in a state of sensitivity and unbiased readiness to respond quickly to any odor that might arrive. With inhalation the rise in level of receptor input increases the non-linear feedback gain within the assembly, so that it is driven away from this basal attractor. The odor input also determines the initial conditions for a state change represented by a bifurcation to a different attractor. For each odor that the rabbit has learned to identify there may be a limit cycle attractor based on synaptic changes during the learning process. The formation of an attractor, we propose, is revealed by changes in EEG spatial pattern with olfactory conditioning (Freeman and Schneider 1982; Viana Di Prisco and Freeman 1985; Gray et al. 1986). The domain of input during the learning period establishes the basin for that attractor. If on subsequent inhalations the inspired air contains that odor, the assembly is placed in the basin for that attractor, and the stereotypic spatial pattern that emerges may suffice to encode that odor, irrespective of which among the set of equivalent receptors received the odor on a given trial. This operation constitutes figure completion and generalization over equivalent stimuli. In this view sensory selection consists of a learned limit cycle attractor and its basin for receptor input (Freeman and Skarda 1985).

In the control state when the animal is breathing air without conditioned stimulus odors, a reproducible spatial pattern results that may suffice to encode the background odor complex. This stereotypic control burst may represent confirmation by the olfactory apparatus of the persistence of an unchanging background. When novel or unexpected odors are given, there is evidence that the assembly is driven from low-level chaos to a high-level attractor that lacks a stable spatial form. Such 'chaotic bursts' may signal that the bulbar assembly has failed to converge to a learned attractor of a known odor, but that an unidentified odor may be present that will require orienting or re-testing. In other words, the bulb may have the capability for reporting: 'test is inconclusive, re-testing is imperative.' According to this interpretation inconclusive outcomes also occur sporadically during control periods, and more frequently during testing with conditioned odors.

Success in the description of these self-organizing properties of nerve cell assemblies in hierarchies will have major impact on several fields of study. Most obviously, the neurological basis for the study of perceptual and movement disorders will be improved. Understanding of mechanisms controlling fetal development of the brain will be enhanced; it is obvious from fetal movements in utero that the brain is highly active during its formation, and that the activity patterns arising from within must contribute to formation of patterns of organization of connections. The structure of the adult brain may in large measure reflect the stored patterns of past activity, just as we believe that the patterns of neural activity during learning leave their traces in synaptic changes that guide future activity. Studies by psychologists of mental images, attention, and voluntary motor activity will be greatly assisted by detailed understanding of the physiological bases for these properties underlying human and animal behavior (Freeman 1983). And just as the design and early development of digital computers and automata took place in the context of theories of the nervous system in the 1940s, the next major developments in machine intelligence may occur in the context of the mathematics of self-organizing nerve cell assemblies.

Appendix A. Linear differential equations for ERPs and PSTHs.

The procedure for constructing equations from data lies at the cutting edge of analysis. We illustrate it first for the passive membrane in order to draw the principles from the most widely known instance.

Suppose that a pair of electrodes is placed on each side of the membrane, one for passing a transmembrane current impulse (the input, I(t) at t = 0), and the other pair nearby on each side of the membrane, for measuring the transmembrane potential difference (the output, v(t) or simply v). Within the linear range v jumps to a new value v₀ with the pulse and then decays monotonically back to the baseline that existed before the pulse, say v = 0. The response wave-form seems to conform to an exponential curve:

This equation is the solution to a first order linear differential equation, using _t HETA to denote dv/dt:

where k_ee and a are constants to be determined. From knowledge of physical systems in general, and of the chemical and anatomical structure of membrane in particular, a model is devised in which the membrane is represented by a resistance, R, in parallel with a capacitance, C. A time constant is defined, a = 1/RC, and v₀ = kI = I/C. Then the theoretical curve for a pulse input,

is fitted to the experimental data by linear regression in semi-logarithmic coordinates in order to measure R and C. Finally, we interpret R and C to mean the existence in the membrane, respectively, of pores through which ions may pass and of stored charge on or in the lipid layer acting as a dielectric.

In review, the 3 steps are: (1) making an observation of an event (the wave shape); (2) constructing, evaluating, and solving a differential equation; and (3) interpreting the parameters of the equation in terms of subsystems. The 3 steps are used in each of the following 3 illustrations.

The effect of axonal, synaptic and cable delays is illustrated for a KO_i set. Suppose that a stimulating electrode is placed on the lateral olfactory tract (LOT) to stimulate antidromically the axons of the mitral cells so as directly to excite the bulbar Kl_i mass comprised by granule cells. In the appropriate experimental conditions (see Ch. 4 in Freeman 1975) the response of the k_ij set is recorded from a monopolar electrode on the bulbar surface with respect to a distant reference electrode. In waking or lightly anesthetized animals the impulse input initiates prolonged interactions, but under very deep anesthesia the interactions are suppressed, and the KI_i mass is reduced to the KO_i state.

With the impulse input, I(t), a rapid increase is observed in surface-negative potential followed by a less rapid decay (Fig. 31, upper curve). For simplicity, the factors are omitted for afferent temporal dispersion in the input tract and the late positivity or dendritic after-potential. By appropriate stimulation with pairs of stimuli to test for proportionality and additivity, a linear range for observation is established. The response can be fitted with an equation in the form

Fig. 31. The data show predicted (curves) and observed (symbols) bulbar responses (averaged evoked potentials or ERPs) under deep pentobarbital anesthesia (the open-loop state) on single-shock LOT (a) or single-shock PON (b) electrical stimulation. (From Freeman 1975, pp. 309-314.)

The first term represents a down exponential curve, v = v₀ at t = 0, and conforms to the decay of the response. The decay rate is a. The second term represents an upward exponential curve, v = -v₀ at t = 0, and the coefficient b determines the rate of rise. When the second curve is subtracted from the first, the amplitude at t = 0 is zero, and the peak latency is seen to depend on both a and b. This is the characteristic form, of an EPSP as well as the simplest form of cortical evoked potentials. The curve must be fitted to the data using non-linear regression, because the form corresponds to the sum of two exponentials.

Eq. 4 is the solution to a linear second-order differential equation,

where k_ee is a gain coefficient. A KO model is devised to represent the KO_i mass as a large number of neurons in parallel, each of which performs 3 operations on its input: input connection strength represented by k_j; combined axonal, synaptic and cable delay represented by b_j; and passive membrane decay represented by a. The KO set represents an average neuron in which the passive membrane model is one of its parts: a = 220/sec. - 1/4.5 msec (the passive membrane decay), and b = 720/sec = 1/1.4 msec (lumped axonal, synaptic and dendritic cable delay) represent the average rate constants. When the stimulus intensity is changed over a wide range, the response amplitude changes proportionately, but the rate coefficients do not change. This is evidence for function in a linear dynamic range. When extrapolated to the waking state it is a basis for fixing the rate constants in higher order equations and allowing only the gain coefficients to change. The interpretation is that the time constants of neurons are determined by the biophysical properties of cell membranes and by dendritic architecture, that these do not change rapidly with behavior, and that functional changes in relation to behavior are based on changes in strengths of interaction within neural masses represented by the k's.

The procedure is repeated under deep anesthesia (what is called the 'open-loop state,' because interactions are suppressed), but this time with stimulation of the KO_i set orthodromically by way of the primary olfactory nerve (PON) and the mitral cells (the bulbar KI_e mass). The output is still that of the granule cells, but the rise time of the response to the impulse is slower (Fig. 31, lower curve). This reflects the delay contributed by the synaptic, cable, and passive membrane properties of the KO_e mass in series with the KO_I mass. We construct a system of 2 second order linear differential equations:

Eq. A6a says that the input k_PI(t) to the KO, set is an impulse, and the output of the KO_e set V_e has the form specified by Eq. A4. Eq. A6b says that the input k_ev_e to the KO_i set is the output of the KO _i set v, times its gain coefficient k_e, and the output of the KO_i set v_i can be predicted by solving Eqs. A6a and A6b for v _i.

The solution gives a curve, which is fitted by non-linear regression to the observed wave-form to evaluate a and b. Experimentally a, = a_i and b, = b_i, so that the equation for the fitted curve (Fig. 31, lower curve) is:

Further details are given elsewhere (Freeman 1975).

When the values are determined for the open loop rate constants a and b, the form of the response of the KII set can be predicted for the waking 'closed-loop state' when interactions are present. For simplicity, the model is restricted to the 'reduced' KII set, in which the interactions within the Kl_e and KI, sets cancel each other. The proof that this occurs and the relevant experimental conditions are described elsewhere (Freeman 1975). The outputs are taken from the KI_i set to simulate the bulbar ERP, and from the Kl_e set to simulate a bulbar poststimulus time histogram (PSTH). For PON stimulation, the differential equations are

Eq. A8a says that the input to the KI, set is the sum of an impulse k_pi(t), minus the output of the Kl_i set v_i times its gain coefficient k_ee,. The KI_i input is inhibitory, so the sign of k_i is negative. Eq. AM says that the input to the KI_i set is given by the output of the Kl_e Set times its gain coefficient k_e. The sign is positive because the Kl_e output is excitatory. For LOT stimulation the input term k₁I(t) is added to Eq. A8b and the input term is deleted in Eq. A8a.

This pair of equations describes a negative feedback loop (Fig. 9). The solutions consist of damped (exponentially decaying) sine and cosine waves with frequency w in radians per second, phase f in radians, and decay rate cc in seconds:

Eq. A9a predicts the closed-loop output of the Kl_e set (the PSTH), and Eq. A9b, that of the KI_i set (the ERP). These predict that the activities of the KI_e and k_ij masses should oscillate at the same frequencies and decay at the same rates but with 90° of phase lead for the output of the KI_e set over that of the k_ij set. This should hold both for PON and for LOT input. The values for w and a should depend on the site and intensity of input (Fig. 32).

Fig. 32. The data show predicted (curves) and observed (symbols) bulbar responses (ERPs) from granule cells (KI_G mass) to single-shock electrical stimulation of LOT (left) or PON (right) under light anesthesia. A: poststimulus time histograms (PSTHs) from mitral cells (KI_m mass) derived concomitantly are fitted with curves simulated with the K11_MG set. (From Freeman 1975, pp. 314-321.)

The basis for the oscillation and for the KI_e phase lead is intuitively explained as follows. Excitation of the KI_e mass excites the KI_i mass, which then inhibits the KI_e mass. Inhibition of the KI_e mass disexcites (inhibits) the k_i mass, which disinhibits (excites) the KI, mass. The cycle repeats with a time period that depends on the rate constants, a and b, and on the product of k_ek_i, which is called the negative feedback gain, K_n. If we know a, b, w₁, and a₁, then we can evaluate K_n, which is the numerical estimate or measurement of the strength of interaction between the KI_e and KI_i masses (see Fig. 10).

For contrast, the form of the response is predicted for the KIP set (Fig. 9) upon impulse stimulation of the PON. The neurons are mutually excitatory (Freeman 1975; Martinez and Freeman 1984), so this is a case of positive feedback. The differential equations are

Eq. A10a says that the input k_PI(t) to the KI_P set is received by a subset with active state p_r that is not transmitting at the time of I(t) at t = 0, and that this subset also receives from the transmitting subset with active state p_t times k_t . Its output is p_r times its gain k_r, Eq. A10b says that the input to the transmitting subset p_t is the output p_r of the receiving subset times its gain coefficient k_r

The solutions for the outputs of the two subsets are

where p_rj, p_tj, and a_j(j = 1, ..., 4) are real or complex. The 4 rate coefficients result from the combination of 2 second order differential equations, and they are identical in Eqs. A11a and A11b. The total output of the KI_P set is p = p_t + p_r.

The predicted output (Fig. 33) consists of a sharp increase in pulse density followed by a slow decay to the prestimulus baseline without overshoot. The rate of increase depends on whether or not the subset under observation receives the PON input volley. The decay rate to the baseline, al, depends on a, b, and the product of k_rk_t , which is the positive excitatory feedback gain, K_P. From the values of a_j, p_rj, and p_tj, j= 1, 4, the strength of interaction, K_P.) can be evaluated in the KI_P mass (Appendix B).

Fig. 33. The plotting symbols show the digital values of the PSTHs taken from a periglomerular cell in the olfactory bulb at 6 successively increased stimulus intensities. With increased input strength and output magnitude the response decays more rapidly to the baseline without overshoot. The dashed curve is the best fit of an empirical set of basis functions by the least squares criterion. The two solid curves are the predicted outputs of the forward and feedback limbs of the KI_P set as shown in schematic form in Fig. 10. It is apparent that at low stimulus intensity this unit is not directly excited by the primary olfactory nerve (PON) input but secondarily by its neighbors. At higher stimulus intensity it is directly excited. For analysis see Fig. 34. (From Freeman 1975, pp. 285-291.)

Appendix B. The root locus method applied to ERPs and PSTHs

The results of measurements of ERPs with linear basis functions consist of the numbers of amplitudes, frequencies, phases and decay rates of the sum of terms in an equation containing sines, cosines, exponentials, and so forth. We need displays of the sets of coefficients from the ERPs comprising each AERP. The most useful display is called the s-plane; it derives from the fact that the equation of the curve fitted to the data is the solution of an ordinary linear differential equation with an impulse input. This differential equation can be written several ways: as an n-th order set of derivatives with respect to time; as a set of n first order time derivatives; or in operational form as an algebraic equation by use of the Laplace transform, in which the differential operator d/dt is replaced by the Laplacian operator s. In fact the simplest and easiest way to construct the differential equation is to look up in a table the Laplace transform of each of the sum of terms in the equation used for curve-fitting and assemble them in the form of a factored polynomial.

The general form, of the differential equation in s is the ratio of two polynomials. Where v(t) represents an EEG, ERP or PSTH amplitude as a function of time, and V(s) represents the same as a function of s, for input I(s).

For ERPs and poststimulus time histograms (PSTHs) I(s) = I times the stimulus intensity. The coefficient K is a real number representing amplification or gain. The coefficients p_i represent the characteristic frequencies and decay rates of the neural system that are called the poles. The coefficients z_j represent the frequencies and decay rates at which activity cannot occur, as already noted in spectral analysis of ERPs (Section 18.5). The bases for these manifestations of 'forbidden zones' are called the zeroes of the system. The poles and zeroes together comprise the roots of the system. The zeroes are not explicit in rate coefficients from ERPs but are implicit in the relative amplitudes and phases of the several terms in the equation of the fitted curve. Identification and evaluation of the zeroes may become a major part of the task of synthesizing a model.

Zeroes appear in the equations for neural systems mainly for 3 reasons. First, the poles of the transfer function of the feedback limb in a loop always re-appear as zeroes in the closed loop transfer function (Fig. 34). Second, rational approximations for distributed delay and dead time always introduce zeroes into the right half of the s-plane. Third, summation at convergence points of two or more pathways introduces zeroes. For example, the damped cosine wave components of the currents from evoked potentials of two neural populations sum in the volume conductor; they may beat together and partially cancel each other. The sum of their transfer functions, each with a complex pole pair, contains 4 complex poles and 2 complex zeroes (Freeman 1975, Ch. 2 and 5).

Fig. 34. Two plots of the s-plane are shown in order to demonstrate the procedure for estimating closed loop gain from measurements of the PSTHs from a Kl_e set stimulated electrically at 6 different intensities. The upper frame shows the poles derived by fitting sums of two exponentials and a damped cosine to the PSTHs; the zeroes cannot usually be evaluated explicitly. The dark curves show the estimated locations of the closed loop root loci for the KI set. The lower frame shows the full root locus diagram after the locations have been estimated for the open loop poles, which become the closed loop zeroes with non-zero loop gain. Selected gain contours are shown, including those for the minimal and maximal gains determined from the PSTHs and for unity gain. The loop gain increases with decreasing response amplitude and by extrapolation it reaches unity at zero amplitude. This is evidence that the KI_e, set is stabilized by its intrinsic non-linear gain, and that an inhibitory interneuron population is not required, although its presence is not thereby ruled out. (From Freeman 1975, pp. 291-305.)

Each of the roots is either a real number in the case of exponential decays, or a complex number in the case of damped sine or cosine waves. In the latter case the frequency w and decay rate a are transformed by Euler's theorem to a pair of complex conjugate poles or zeroes, for example, p_i = a + jw and p_i+1 = a - jw, where w = 2πf, and a = 1/t, a time constant. The s-plane for display has the coordinates a and jw (Figs. 8, 10 and 11). The real numbers for s = -z_i and s = -p_i lie on the a-axis, negative to the left. The complex numbers lie in the plane at s = -a ± jw. Owing to mirror image symmetry only the upper half of the s-plane need be shown. Then the results of measurement for each ERP take the form of a constellation of roots (open loop poles marked by X, closed loop poles by ∆, open loop zeroes by 0, closed loop zeroes by ) in the s-plane.

We will now illustrate the method for evaluating loop gain K_P. The set of frequencies and decay rates of the dominant component are shown as triangles; the s-plane in Fig. 34. The upper frame shows the locations of the poles derived by fitting a curve to each of the PSTHs shown in Fig. 33. The zeroes are not yet defined or evaluated. The curves represent the root loci that are to be defined. A differential equation is formulated with the two identical KO sets arranged in mutual excitation. This is solved for different values of the open loop rate constants a and b and for a set of values of feedback gain K_P. The rate constants are determined by trial and error in order to place the root locus through the complex poles. The gain is evaluated by approximation to force the contour for a gain value through the real axis at the location of a real-valued pole nearest the origin. The same gain contour should then pass through the complex-valued poles from the same PSTH. The gain contours for unity gain and for the minimal and maximal gains of the PSTHs are shown as the closed curves in Fig. 34, lower frame. The open loop poles become the closed loop zeroes of the non-zero feedback system. These are shown as small squares on the negative real axis. Further details with more examples are to be found elsewhere (Freeman 1975, Ch. 2, 5 and 6).

Appendix C. Non-linear integrodifferential equations for EEGs

The circuit diagram for the KII subset is shown in Fig. 9. As described in Section 18.8 each KI subset is replaced by 2 KO subsets. Each component KO subset has 2 state variables: u(t) in the dendritic current or wave mode and p(t) in the axonal impulse density or pulse mode. The conversion from pulse to wave at synapses is linear and is represented by a synaptic coefficient. The conversion from wave to pulse is a non-linear operation that is represented by the sigmoid curve in Fig. 35 (left frame).

The equations are expressed in normalized pulse density Q and in normalized wave amplitude v with respect to a basal level of wave bias u_o representing a scalar quantity that varies with the level of arousal or motivation.

By appropriate choice of units of measurement, p_o = u_o and

The normalized equations for Fig. 35 are

Fig. 35. The left frame shows the non-linear relation between dendritic wave amplitude and impulse density at trigger zones of neurons in a KO set from Eq. 14. The 3 triangles show the rest state at v = u_o and p = p_o, the mean background pulse density. The lowest curve corresponds to the state under light anesthesia, the middle to a waking but unmotivated state, and the higher curve to a mild degree of motivation, such as that induced by a few hours of food deprivation. The right frame shows the corresponding gain curves from Eq. 16. Maximal gains may exceed values of 20; the sensitivity of the KI_e set may vary over a range of 1-40,000. (From Freeman 1979a.)

The lower asymptote is determined by the distribution of thresholds in the KO subset. The exponential increase is determined by the properties of the axonal trigger zones just below threshold in the Hodgkin-Huxley system. The upper asymptote is determined by the refractory periods and long-lasting recovery processes for all neurons in the subset (Freeman 1979a). The maximum Q_m is also a function of level of background bias u_o, so that the form of the curve depends on the value of a single parameter (Fig. 35). It is evaluated by fitting the curve to experimental data (Freeman 1979a). From the fact that pulse densities cannot be negative,

The non-linear gain (Fig. 35, right frame) is given by

The wave amplitude v_max at which the maximal gain occurs is displaced to the excitatory side.

It is this property that underlies the parametric change in KII feedback gains with changes in level of input, because excitatory activity in the KI_e set is self-regenerative.

The linear operations of axonal, synaptic and dendritic cable delays and the passive membrane decay are approximated by a linear second order ordinary differential equation.

The 4 elements in each subset are linked with other elements according to Fig. 9 in the following equations.

where N is determined by measurements of spatial divergence in the nerve cell assembly (Freeman 1975, Ch. 4), and by the boundary conditions, either circular for the 1-D case or toroidal for the 2-D case (Freeman 1979c). The strengths of action k_ij are defined and described in Section 18.7. The superscripted excitatory gains k_ee^jk are modifiable in order to represent synaptic changes during associative learning. The others are fixed. The input from receptors is designated I_j for the j-th element. The superscripted parameters z_e can be modifiable in order to represent the synaptic operation of habituation.

Examples of the solutions to these equations for input that simulates the time course of receptor input density to the bulb during inhalation are shown in Fig. 14. Simulation of chaotic activity requires a structural change in the model, that can be described generically as an increase in the positive and negative feedback within the KIII set (Section 18.1). This aspect of the model has not yet been adequately developed (Freeman 1986, 1987).

Appendix D. Estimation of dimensions of chaos in EEGs (with K.A. Grajski)

Techniques for demonstrating and characterizing chaos in experimental preparations arise from the study of theoretical dynamical systems, that is, models with solutions that describe the time evolution of the states of the systems based upon a knowledge of initial conditions and state transition laws. Associated with each model is a geometrical 'phase space' representing the entirety of states which the modeled system may occupy. (The dimension of phase space is equal to the number of measures required to specify an arbitrary initial state for the modeled system.) The time evolutions of the states for this system are represented by trajectories from an initial point in phase space. Deterministic chaos or turbulent behavior of trajectories has been described for many dynamical systems and shown to differ fundamentally from random behavior (Gurel and Rössler 1979; Helleman 1980). The techniques developed to characterize chaotic behavior and differentiate it from random activity in dynamical systems can, in principle, be applied to experimental data (Rössler 1983).

The behavior of trajectories is described both temporally and geometrically. Their relationship is analogous to that between temporal and spectral time series analysis; they are fundamentally related through the Fourier transform, and each captures a unique aspect of a time series. With each type of asymptotic temporal behavior there is associated a corresponding subspace of phase termed the 'basin of attraction,' that has the property that once a trajectory falls within it, it will, barring changes in the parameters of the system, remain there. Asymptotically, the temporal behavior of the trajectory is that of the associated attractor.

The basis for a temporal description of trajectories is an analysis of the stability characteristics of local regions of phase space. A trajectory may diverge from, remain fixed on, or contract along each of the N dimensions corresponding to any subregion of phase space. The spectrum of N Lyapunov characteristic exponents quantifies these tendencies. The relationships between these exponents and various measures of the dimension of attractors has received much study (Farmer 1982). However, the nature of these relationships is in question; the exponents are, difficult to compute and do not yield as much information as do geometrical approaches (Grassberger and Procaccia 1983).

Geometrical (phase space) descriptions of random perturbations from a stable state differ fundamentally from those for chaotic activity. For a stable state the phase space volume contracts within the basin towards the attractor. That is, local neighborhoods corresponding to 'close' initial conditions give the same asymptotic behavior and are bounded within the structure of the attractor. The simplest stable state is associated with a point attractor of dimension 0. A limit cycle attractor is dimension 1. Chaotic -activity, on the other hand, may stretch and fold as well as contract phase space volume. The stretching corresponds to divergence along particular dimensions, i.e., instabilities. The folding is owing to the fact that the activity is bounded within the attractor by the presence of convergence along other dimensions, i.e., stability. While chaotic activity may appear random, it is strictly determined by the interplay of these convergence and divergence properties. The structure of chaotic attractors is described with a variety of measures of 'dimension' (Farmer 1982). The most commonly applied is the Hausdorff-Besikovich or fractal dimension (Mandelbrot 1982). Experimentally, a fractal dimension lower than the dimension of phase space suggests that the dynamics of the system behavior can be modeled within the subspace belonging to a chaotic or 'strange' attractor (Farmer 1982).

Before we apply any 'dimension seeking' algorithm, we must specify the relationship between experimental data and the phase space of a corresponding dynamical system. In this appendix we use field potentials to construct phase spaces by fitting to them the solutions of differential equations. Similarly, with the array recording technique each channel is interpreted as sampling the activity along one dimension of phase space. The corresponding dynamical system is the distributed KII set (Sections 18.1 and 18.7). At each time the spatial sample is a measurement of a point in a 64-dimensional space. The search for chaos in such data is a search for a spatial chaotic attractor. An alternative approach reconstructs phase space from single time series. (Packard et al. 1980; Grassberger and Procaccia 1983). The main requirement is several thousand samples of activity associated with the attractor. However, this approach is less appropriate for the olfactory EEG, because its activity appears to shift repeatedly between limit cycle and chaotic attractors. Since an unbroken stretch of EEG reflects the contributions of several attractors, the resulting measures of dimension may be misleading. The advantage of a multi-channel approach is that phase space can be sampled spatially during presumed activity owing to a particular attractor. In either case the time series is segmented with separate evaluation of burst and interburst states. A limitation shared by both approaches is the 'curse of dimensionality.' As the presumed dimension of phase space increases, the number of samples required to attain a robust result increases exponentially.

One method for assessing the fractal dimension of an attractor is the algorithm due to Grassberger and Procaccia (1983). This algorithm exploits the geometrical difference between random perturbations and chaotic activity to calculate an estimate of dimension. Due to the stretching of phase space in an attractor, points become temporally uncorrelated. Owing to the folding, however, points are spatially correlated. The spatial correlation is measured by:

where Y is the number of pairs (i,j) such that |_i - _j| < r; that is,

Using each available point X_i, once as a reference, the number of points which fall within a distance r is counted. Selecting an appropriate distance range is a 2-step process. An initial 'broad band' analysis with large distance increments is used to assess the average distance between points. With these averages as the middle values, a 'fine tuned' analysis (100-1000 distance increments) is applied. The result is averaged over all points; the procedure is repeated incrementally over a range of distances.

Grassberger and Procaccia show that for small r,

The exponent v is shown to be a lower limit for the fractal dimension and is termed the correlation exponent (CE). Random noise is shown to give CE = M in M dimensional space (M = 1, ..., 64). If the value of the exponent calculated is less than M and is fractal, then the data stems from deterministic chaos. With increasing distance r the quantity C(r) increases from zero to some maximal value <M. The slope of this curve is approximated in log • log coordinates by linear regression on the points in the range of r wherein the slope is maximal and near linear.

As an example 360 bursts (38 x 64) recorded from a rabbit over 3 sessions are divided into 3 groups by the criteria described in Section 18.12 (Fig. 29, upper left): coherent control bursts, coherent odor bursts (both presumed limit cycle activity), and incoherent bursts (presumed to be chaotic activity). These criteria are the frequency of oscillation, and spatial and temporal coherence measures discussed earlier. Assessment of the dimension of presumed limit cycle activity serves as one control for the application of the Grassberger-Procaccia algorithm. We predict that limit cycle activity will have an integer correlation exponent (Farmer 1982; Grassberger and Procaccia 1983). However, the best experimental data are subject to random disturbance which may affect measurement with the algorithm (Ben-Mizrachl et al. 1984). The degree of disparity of the CE from an integer value may suggest a 'confidence interval' for subsequent measures. Another control is to generate random vectors in 64-space where the total number of points is similar to that for the 3 groups. The degree to which the correlation exponent differs from 64 indicates the strength of the 'curse of dimensionality.' If there is a breakdown in equality at some value of dimension lower than 64, then that value defines the upper limit of resolution for that sample size.

Preliminary results with N = 6150 samples from coherent control bursts under air stimulation yields CE = 4.096; 2770 samples from coherent bursts under odor stimulation yields CE = 4.881; 4750 samples from incoherent bursts gives CE = 4.578. The presumed limit cycle activity has values 'close' to the integers 4.0 and 5.0, suggesting a 'confidence interval' of 0.12. For 4000 random vectors generated with a Normal (0, 1) generator embedded in M dimensions, the Grassberger-Procaccia algorithm gives CE = M (within ±0.05) up to at least 8 dimensions. Hence, the present range of correlation exponents calculated for this EEG data appears well within the limits of confidence in the performance of the algorithm. By these criteria the CE for the incoherent bursts is fractal, and we conclude that the incoherent bursts are chaotic. The validity of these estimates of the dimensions and their physiological meaning remain to be determined.

Acknowledgements

Supported by Grant MH06686 from the National Institute of Mental Health, and by Grant NS16559 from the National Institutes of Health, United States Public Health Service.

References

Abraham, R.H. and Shaw, C.D. (1983) Dynamics, the Geometry of Behavior, Parts I and 2. Ariel Press, Santa Cruz, CA.

Babloyantz, A. and Kaczmarek, L.K. (1979) Self-organization in biological systems with multiple cellular contacts. Bull. Math. Biol., 41, 193-201.

Basar, E. (1976) Biophysical and Physiological Systems Analysis. Addison-Wesley, London.

Ben-Mizrachi, A., Procaccia, 1. and Grassberger, F. (1984) Characterization of experimental (noisy) strange attractors. Phys. Rev. A, 29, 975-977.

Bressler, S.L. (1984) Spatial organization of EEGs from olfactory bulb and cortex. Electroenceph. clin. Neurophysiol., 57, 270-276.

Childers, D.G. (1978) Modern Spectrum Analysis. IEEE Press, New York.

Eastman, C. (1975) Construction of miniature electrode arrays for recording cortical surface potentials. J. Electrophysiol. Technol., 5, 28-30.

Emery, J.D. and Freeman, W.J. (1969) Pattern analysis of cortical evoked potential parameters during attention changes. Physiol. Behav., 4, 69-77.

Farmer, J.D. (1982) Chaotic attractors in an infinite-dimensional dynamical system. Physica, 9D, 366-393

Freeman, W.J. (1962a) Comparison of thresholds for behavioral and electrical responses to cortical electrical stimulation in cats. Exp. Neurol., 6, 315-331.

Freeman, W.J. (1962b) An ergometer for measuring rate of work from cats as an index for drive. J. appl. Physiol., 14, 1071-1072.

Freeman, W.J. (1964) Correlation of goal-directed work with sensory cortical excitability. Recent Adv. Biol. Psychiat., 7, 243-250.

Freeman, W.J. (1975) Mass Action in the Nervous System. Academic Press, New York.

Freeman, W.J. (1976) Quantitative patterns of integrated neural activity. In: J.C. Fentress (Ed.), Simpler Network and Behavior. Sinauer Ass." Sunderland, MA, pp. 280-296.

Freeman, W.J. (1978) Spatial frequency analysis of an EEG event in the olfactory bulb. In: D.A. Otto (Ed.), Multidisciplinary Perspectives in Event-Related Brain Potential Research. U.S. Government Printing Office, Washington, DC, EPA-600/9-77-043, pp. 531-546.

Freeman, W.J. (1979a) Nonlinear gain mediating 'cortical stimulus-response relations. Biol. Cybernet., 33, 237-247.

Freeman, W.J. (1979b) Nonlinear dynamics of paleocortex manifested in the olfactory EEG. Biol. Cybernet., 35, 21-34.

Freeman, W.J. (1979c) EEG analysis gives model of neuronal template-matching mechanism for sensory search with olfactory bulb. Biol. Cybernet., 35, 221-234.

Freeman, W.J. (1980) Use of spatial deconvolution to compensate for distortion of EEG by volume conduction. IEEE Trans. biomed. Engng, 27, 421-429.

Freeman, W.J. (1983) The physiological basis of mental images. Biol. Psychiat., 18, 1107-1125.

Freeman, W.J. (1986) Petit mal seizure spikes in olfactory bulb and cortex caused by runaway inhibition after exhaustion of excitation. Brain Res. Rev., 10, 259-284.

Freeman, W.J. (1987) Simulation of chaotic EEG patterns with a dynamic model of the olfactory system. Biol. Cybernet., in press.

Freeman, W.J. and Baird, B. (1986) Correlation of olfactory EEG with behavior: spatial analysis. Behav. Neurosci., in press.

Freeman, W.J. and Grajski, K.A. (1986) Correlation of olfactory EEG with behavior: factor analysis. Behav. Neurosci., Submitted.

Freeman, W.J. and Schneider, W.S. (1982) Changes in spatial patterns of rabbit olfactory EEG with conditioning to odors. Psychophysiology, 19, 44-56.

Freeman, W.J. and Skarda, C.S. (1985) Spatial analysis, nonlinear dynamics, and perception: the neo-Sherringtonian view. Brain Res. Rev., 10, 147-175.

Freeman, W.J. and Van Dijk, B. (1987) Spatial patterns of visual cortical fast EEG during conditioned reflex in a rhesus working. Brain Res., Submitted.

Freeman, W.J. and Viana Di Prisco, G. (1986a) EEG spatial pattern differences with discriminated odors manifest chaotic and limit cycle attractors in olfactory bulb of rabbits. In: G. Palm and A. Aertsen (Eds.), Proc. Conf. Brain Theory, Trieste, 1984. Springer, Berlin, pp. 97-119.

Freeman, W.J. and Viana Di Prisco, G. (1986b) Correlation of olfactory EEG with behavior: time series analysis. Behav. Neurosci., in press.

Garfinkel, A. (1983) A mathematics for physiology. Amer. J. Physiol., 245, R455-R466.

Golub, G.H. and Pereyra, V. (1973) The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate. J. num. Anal., 10, 413-432.

Gonzalez-Estrada, M.T. and Freeman, W.J. (1980) Effects of carnosine on olfactory bulb EEG, evoked potentials and DC potentials. Brain Res., 202, 373-386.

Gonzales, R.C. and Wintz, P. (1977) Digital Image Processing. Addison-Wesley, Reading, MA.

Grajski, K.A., Breiman, L., Viana Di Prisco, G. and Freeman, W.J. (1987) Classification of EEG spatial patterns with tree-structured methodology. IEEE Trans. biomed. Engng, in press.

Grassberger and Procaccia, 1. (1983) Measuring the strangeness of strange attractors. Physica, 9D, 189-208.

Gray, C.M., Freeman, W.J. and Skinner, J.E. (1986) Chemical dependencies of learning in the rabbit olfactory bulb: acquisition of the transient spatial-pattern change depends on norepinephrine. Behav. Neurosci., in press.

Griffith, J.S. (1971) Mathematical Neurobiology. Academic Press, New York.

Gurel, 0. and Rössler, O.E. (Eds.) (1979) Bifurcation theory and its applications in scientific disciplines. Ann. N.Y. Acad. Sci., 316 pp.

Hassard, B.D., Kazarinoff, N.D. and Wan, Y.-H. (1981) Theory and Applications of Hopf Bifurcations. Cambridge University Press, New York.

Hebb, D.O. (1949) The Organization of Behavior. Wiley, New York.

Helleman, R.H.G. (Ed.) (1980) Nonlinear dynamics. Ann. N.Y. Acad. Sci., 357.

Horowitz, J.M. (1972) Evoked activity of single units and neural populations in the hippocampus of the cat. Electroenceph. clin. Neurophysiol., 32, 227-240.

Ingber, L. (1982) Statistical mechanics of neocortical interactions. 1. Basic formulation. Physica, 5D, 83-107.

Johnson, R.A. and Wichern-, D.W. (1982) Applied Multivariate Statistics. Prentice Hall, Englewood Cliffs, NJ.

Katznetson, R.D. (1982) Deterministic and Stochastic Field Models in the Neurophysics of the EEG, Ph.D. Thesis. Department of Electrical Engineering, University of California, San Diego, CA.

Kaufman, L. (1975) A variable projection method for solving separable nonlinear least squares problems. B.I.T., 15, 49-57.

Lin, C.C. and Segal, L.A. (1974) Mathematics Applied to Deterministic Problems in the Natural Sciences. MacMillan, New York.

Lopes da Silva, F.H., van Rotterdam, A., Barts, P., van Heusden, E. and Burr, W. (1976) Models of neuronal populations: the basic mechanisms of rhythmicity. In: M.A. Corner and D.F. Swaab (Eds.), Perspectives in Brain Research. Progr. in Brain Res., Vol. 45. Elsevier, Amsterdam, pp. 281-308.

MacGregor, R.J. and Lewis, E.R. (1977) Neural Modelling. Plenum, New York.

Mandelbrot, B. (1982) The Fractal Geometry of Nature. Freeman, San Francisco, CA.

Martinez, D.M. and Freeman, W.J. (1984) Periglomerular cell action on mitral cell in olfactory bulb shown by current source density analysis. Brain Res., 308, 223-233.

Meisami, E. and Emamian, S. (1985) Cytoarchitectural diversity in the main olfactory bulb: are there olfactory foveas? Amer. chemosens. Soc., VII, 130 (Abstr.).

Nufiez, P.L. (1981) Electric Fields of the Brain. Oxford University Press, New York.

Packard, N.A., Crutchfield, J.P., Farmer, J.D. and Shaw, R.S. (1980) Geometry from a time series. Phys. Rev. Lett., 45, 712-715.

Pickering, S.G. and Freeman, W.J. (1967) Variations of the superior colliculus evoked response in cats. Exp. Neurol., 19, 127-139.

Rall, W. and Shepherd, G.M. (1968) Theoretical reconstruction of field potentials and dendrodendritic synaptic interactions in olfactory bulb. J. Neurophysiol., 31, 884-915.

Rössler, O.E. (1983) The chaotic hierarchy. Z. Naturforsch., 38a, 788-801.

Sammon, J.W. (1969) A nonlinear mapping for data structure analysis. IEEE Trans. Comput., 18, 401-409.

Simon, W. (1972) Mathematical Techniques for Physiology and Medicine. Academic Press, New York.

Thatcher, R.W. and John, E.R. (1977) Functional Neuroscience. Vol. I. Foundations of Cognitive Processes. Erlbaum, New York.

Turing, A.M. (1952) The chemical basis of morphogenesis. Phil. Trans. roy. Soc. B, 237, 37-72.

Viana Di Prisco, G. and Freeman, W.J. (1985) Odor-related bulbar EEG spatial pattern analysis during appetitive conditioning in rabbits. Behav. Neurosci., 99, 964-978.

Von Neumann, J. (1958) The Computer and the Brain. Yale University Press, New Haven, CT.

Wilson, H.R. and Cowan, J.D. (1972) Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J., 12, 1-24.

Wilson, H.R. and Cowan, J.D. (1973) A mathematical theory of functional dynamics of cortical and thalamic nervous tissue. Kybernetik, 13, 55-80.

______________

END