- Email: [email protected]
Cortical transients preceding voluntary movement
Mathematical Biosciences 222 (2009) 86–91
Contents lists available at ScienceDirect
Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs
Cortical transients preceding voluntary movement J.W. Hartwell * JW Hartwell & Associates, 3001 Hartwell Pond Drive, Hillsborough, NC 27278, USA
a r t i c l e
i n f o
Article history: Received 21 May 2009 Received in revised form 4 September 2009 Accepted 17 September 2009 Available online 25 September 2009 Keywords: Cortical dynamics Cortical transients Voluntary movement Readiness potential Contingent negative variation
a b s t r a c t The process of initiating a voluntary muscular movement evidently involves a focusing of diffuse brain activity onto a highly specific location in the primary motor cortex. Even the very simple stereotypic movements used to study the ‘contingent negative variation’ and the ‘readiness potential’ begin with EEG indicative of widely distributed brain activity. In natural settings the involvement of diffuse cortical networks is undoubtedly even more important. Eventually, however, activity must coalesce onto specific neurons for the intended movement to ensue. Here we examine that focusing process from a mathematical point of view. Using a digital simulation, we solve the global equations for cortical dynamics and model the flow from diffuse onset to localized spike. From this perspective the interplay between global and local effects is seen as a necessary consequence of a basic cortical architecture which supports wave propagation. Watching the process evolve over time allows us to estimate some characteristic amplitudes and delays. Ó 2009 Elsevier Inc. All rights reserved.
1. Introduction The process of initiating a voluntary muscular movement evidently involves a focusing of diffuse brain activity onto a highly specific location in the primary motor cortex. Even the very simple stereotypic movements used to study the ‘contingent negative variation’ (CNV) [1,2] and the ‘readiness potential’ [3,4], begin with EEG apparently indicative of widely distributed brain dynamics. In natural settings the involvement of diffuse cortical networks is undoubtedly even more important [5]. Eventually, however, activity must coalesce onto specific neurons for the intended bodily movement to ensue. Here we examine that focusing process from a mathematical point of view. Using a digital simulation, we solve the global equations for cortical dynamics [6] and model the flow from diffuse onset to localized spike. From this perspective the interplay between global and local effects is seen as a necessary consequence of the basic cortical architecture which supports wave propagation. An understanding of brain activity preceding normal movements could have important practical implications. As prosthetic devices become more fully functional imitations of the limbs and organs they replace, the interfaces for their control take on added significance. To achieve smooth activation of these complex devices, engineers seek ‘intuitive’ control mechanisms, mimicking the natural flow insofar as they are able and interfacing to the body’s neural system as early as understanding and technology permit. ‘Mind reading’ might summarize the goal. If a user’s inten* Tel.: +1 919 732 7951. E-mail address: [email protected] 0025-5564/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2009.09.003
tions for fine motor movements could be reliably detected, then mechanical devices might become much more useful to him. And if that detection took place at the cortical level, then response times could be shortened and prosthetic devices used even by those with severe injuries to the peripheral nervous system. Advances in this technology might prove especially beneficial to socalled ‘locked-in’ patients suffering (mainly) from ALS and brainstem strokes. Already some success with brain-computer interfaces has been achieved, enabling patients to move a cursor on a computer screen [7] and to control an artificial hand [8]. Regrettably little is known about the specific cortical events that precede voluntary muscle movements. The ‘motor strip’ itself is well mapped, and it is clear that highly specific neurons in it must fire to initiate the contraction of particular skeletal muscles. The necessary cortical activity must have a distinctly localized character, focusing on certain cells and excluding others located nearby. But these same neurons are richly interconnected with others in a structure that supports complex wave activity and precludes their firing as isolated entities. Indeed, EEG patterns recorded preceding voluntary movements are characterized by very diffuse waves at the onset. The mathematical simulation reported here represents an attempt to understand and model how this diffuse beginning leads to the eventual sharply localized result. Watching the process evolve over time allows us to estimate some characteristic amplitudes and delays. Both the ‘contingent negative variation’ (CNV) [1,2] and the ‘readiness potential’ [3,4] are EEG waves marked by a negative voltage spread over a wide region at the top of the head. This voltage grows in magnitude over a period of about half a second prior to a muscular response. In the case of the CNV, larger voltages are
J.W. Hartwell / Mathematical Biosciences 222 (2009) 86–91
associated with more intense involvement in the task and with shorter reaction times [2]. Occasionally the CNV is sufficiently large as to be observed in the raw EEG, but its study usually requires time-locked averaging. The readiness potential is similarly small, but somewhat differently distributed. In particular, it is more pronounced over the hemisphere contra-lateral to the responding hand. Alan Givens writes, ‘About 200 milliseconds after the second stimulus, the CNV ends in a sharp positive-going (10 Hz) wave termed the CNV resolution’ [2]. The readiness potential appears to end less suddenly [4]. In any event, this marks the onset of muscle movement, and concomitant inputs to the adjoining somatosensory cortex contribute to the EEG from this point forward. To see how a sharply localized response can evolve in this richly interconnected system, we have examined the global equations for ‘neocortical dynamics’, used successfully by Nunez [6,10,11], Katznelson [9], Nunez and Srinivasan [11] and others to elucidate widespread cortical rhythms as observed in EEG. These equations are based on a columnar model of the cortex with inhibitory interconnections that are predominantly short-range and with excitatory connections that at their longest span a cortical hemisphere. The density of these long fibers decreases with distance, with a characteristic length of about 10 cm. Their finite conduction velocity gives rise to the propagation of waves in synaptic activity. Both Nunez [10] and Jirsa and Haken [12] have joined the global equations to the contributions of local networks in the generation of EEG. The current work stops short of this goal, considering only global relations. The equations are normally written in integral form, with the activity at a particular cortical column being determined by the subcortical and commissural inputs to that location plus an integral over connections from the surrounding cortex, taking activity at earlier times as required by the intervening distance and the conduction velocity. About 92% of the long axon nerve fibers arriving at a cortical column are from other locations in the same hemisphere [9]. Although the inhibitory connections are often associated with a substantial delay, we follow treatments by Nunez in which their effects, along with those of the short-range excitatory fibers, are considered as ‘instantaneous’ (i.e. as comparable to post-synaptic delays within a cortical column, which are also neglected). The dynamics of the system depend in part on a ratio of overall positive to negative feedback as captured in a single parameter, B, which controls the gain for inputs of relatively long wavelength. We have taken B ¼ 0:95 which preserves stability of the linearized equations. Larger values of B lead to more labile behavior and are considered more realistic by some workers, but they require constant subcortical inhibitory inputs to keep the system operating in a quasi-linear region [13]. In the work reported here, the equations have been transformed from their usual integral representation to an equivalent partial differential form better suited to digital simulation (see Appendix A). They are solved on a closed spherical surface with a circumference of 80 cm, corresponding approximately to a cortical hemisphere ‘inflated’ so as to smooth its fissures and undulations. The conduction velocity of the long fibers is taken to be 800 cm/s. The digital simulation used discrete cortical elements with centers separated by no more than 0.8 mm and time steps of 1.0 ms.
2. Initial priming Our simulation begins with a diffuse, low-level input held constant for half a second. The input has a broad maximum over the eventual focus point in one hemisphere and tapers off to zero on its opposite side. The resulting cortical activity takes on a similar pattern and reaches somewhat larger amplitude. The response rises
87
slowly with a time constant of about 100 ms and levels off with a relative gain of about 1.8. (The absolute gain is considerably greater, but we wish to demonstrate a benefit of priming with long wavelengths. As mentioned in Appendix A, we have reduced to unity a multiplicative constant that affects all wavelengths equally.) Fig. 1 shows the steady input (green) and increasing response (red), now nearing its final amplitude at the peak location. This response represents increased post-synaptic activity in pyramidal neurons across a large region. These cells are aligned perpendicular to the cortical surface, and their collective activity is reflected in a similarly widespread negative voltage on the scalp, corresponding to the initial phases of the CNV and readiness potential. 3. Trigger onset At any time after initial priming, motion can be triggered by changing the input over a disc surrounding the focus point. Here a small negative input (i.e. a reduction from prevailing ambient) is applied over a disc of radius 35 mm as shown in Fig. 2. Responses corresponding to inputs of short wavelength are oscillatory and under-damped. A small wave front begins to propagate in both directions from the discontinuity. The out-bound wave decreases in amplitude as it propagates both because of the overall preponderance of negative feedback and because the wave front spreads over an increasingly larger perimeter. The in-bound wave, however, grows in amplitude as it concentrates over a decreasing perimeter, this effect out-weighing that of the negative feedback. A disc with radius 35 mm is near the optimum for the ratio of positive to negative feedback used in the simulation ðB ¼ 0:95Þ. For larger discs the inbound wave initially loses amplitude, while smaller ones would benefit from a greater circumference of initial activity. 4. A big bounce The wave of reduced activity grows in amplitude as it propagates toward the disc center as shown in Figs. 3 and 4 (red1 traces) where it converges to a trough after 62 ms. The locally reduced activity is reflected in a momentarily positive scalp voltage. The system will rebound in a few milliseconds. All input was removed 5 ms after the trigger pulse was applied. (Green traces in Figs. 3 and 4 show zero.) The consequent second wave of opposite polarity is also converging and will reinforce the oscillatory maximum. 5. Focusing pulse As the effects of the trigger pulse and its removal converge, a small focusing pulse is applied to select the particular motor neurons necessary for the desired muscle movement. Being of very short spatial extent, the effect of this pulse is almost immediate. It combines with the waves initiated earlier to drive the activity to a large value at a highly specific location. This is shown in Fig. 5 where the overall pulse amplification is more than 13:1. The narrow width of the resulting activity places it beyond the validity of the global wave equations. A more detailed accounting would require considering interactions with the local motor network [10]. Although the trigger and focusing pulses are of low amplitude, their relative timing is crucial. For the 35 mm trigger disc used in our simulation, the focusing pulse must follow the trigger by 69 ms to achieve its greatest effect. This constitutes an appreciable fraction of the very fast reaction times often recorded in experiments in which the CNV is prominently seen. 1 For interpretation of the references to color in Fig. 4, the reader is referred to the web version of this paper.
88
J.W. Hartwell / Mathematical Biosciences 222 (2009) 86–91
Fig. 1. Simulation through 503 ms. Initial priming and trigger at focus. Green trace shows input, constant for half a second, with sudden trigger onset. Red trace shows gradually increasing cortical activity now nearing its final amplitude. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)
Fig. 2. Cortical slice at 503 ms. Green trace shows input with trigger extending over a disc of radius 35 mm. Red trace shows activity. Small waves begin to propagate away from the discontinuity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)
Fig. 3. Simulation at 532 ms. Inbound wave increases in amplitude; outbound wave decreases. All input was removed at 505 ms (green trace shows zero). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)
Fig. 4. Simulation at 562 ms. Inbound activity has converged to a deep trough. Input remains zero.
Fig. 5. Simulation at 569 ms. Extremely narrow pulse of activity follows a small focusing input. Green trace shows spatial distribution of input and red the resultant activity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)
89
J.W. Hartwell / Mathematical Biosciences 222 (2009) 86–91
Fig. 6. Recovery at focus. Response to trigger and focusing pulse decays with rapid oscillations. Initial priming disappears with a long time constant.
Fig. 7. Simulation at 1000 ms. Only a portion of the initial priming remains.
6. Recovery The duration of the focusing pulse is not critical. The spike in cortical activity is extremely narrow in time as well as in space, persisting for only a few milliseconds, less than the refraction time of the motor neurons at its focus. The dynamics during this short interval must in fact be governed more by the local motor network than by the global wave theory used here. Recovery from the trigger and focusing pulses again involves global wave propagation, however. The global recovery requires about 80 ms during which the decay shows oscillations at approximately 150 Hz (too high to be seen in normal EEG). See Fig. 6. This, together with the 69 ms lead time required by the trigger pulse, gives a duration of around 150 ms for the firing process. Because the system is operating in a quasi-linear mode, there is no corresponding time delay required between successive reenactments. A new sequence could be initiated at any point. The long wavelength initial priming disappears in about half a second after its input is removed, decaying with the same long time constant that governed its build up. The state after 1000 ms is shown in Fig. 7. 7. Summary The initial diffuse wave provides a platform of coordinated activity that allows more sharply defined inputs to be magnified in effect. The mechanism benefits from the high gain associated with static inputs of long wavelength. The transition from this global preparation to the highly discrete final result utilizes converging waves initiated by a small but more localized trigger pulse. The in-bound pattern may be visualized as similar to what would result from throwing a rock into a pond and then playing the movie of the resulting waves backwards. Finally, the precise focus is determined by local interactions, represented here by a small short pulse whose effect is greatly amplified by the waves preceding it. The complete sequence utilizes the dynamic character of the global cortical interconnections somewhat like a gymnast exploits wave motion in the elastic medium of a trampoline, using preparation and precise timing to achieve notable effect. 8. Some consequences The good news for prosthetics engineers in this simulation is that transients associated with the trigger for a voluntary movement be-
gin about 69 ms prior the firing of the specific motor neurons involved. This gives an opportunity to detect and respond. The bad news is that the waves begin small and distributed over an area that covers much of the motor strip. Most easily detected is the widespread excitatory priming activity which seems only to facilitate motor output generally. The trigger, which lends the first hint of specificity, consists in a broad disc which gives rise to a slowly collapsing annulus. The specific point of focus stands out only very late in the sequence. Sensors and circuitry to detect such a developing pattern might have to resemble the cortical wiring itself. Acknowledgments Thanks to Paul Nunez and the reviewers for detailed and helpful commentary. Appendix A Following Nunez’ ‘Solutions to the Equations for Global Neocortical Dynamics’ as presented in Appendix A to Ref. [6], we adopt similar notation and write the linearized (small signal) equations for excitatory and inhibitory activity at cortical location ðx; yÞ and at time t as
HE ðx; y; tÞ ¼ Uðx; y; tÞ þ
N X
ðq00 k200 =2pÞ
n¼1
Z Z S0
½Q E HE ðx0 ; y0 ; t r=v Þ 0
0
0
0
Q I HI ðx0 ; y0 ; t r=v Þ expfkn rg dx dy ; Z Z HI ðx; y; tÞ ¼ ðqI k2I =2pÞ ½Q E HE ðx0 ; y0 ; t r=v Þ S0
Q I HI ðx0 ; y0 ; t r=v Þ expfkI rg dx dy ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ðx x0 Þ2 þ ðy y0 Þ2 : ð1Þ Here Uðx; y; tÞ denotes excitatory inputs arriving from the opposite hemisphere or from sub-cortical structures at the same time and place. The integrals sum the contributions of cortico-cortical fibers from the entire surface of one hemisphere, weighting each region by the density of fibers from it (which decrease exponentially with distance), and using activity from the earlier time t r=v to account for the finite neural propagation velocity v. The coefficients Q E and Q I give the sensitivity for small increases or decreases in excitatory and inhibitory inputs, respectively. kI is the reciprocal
90
J.W. Hartwell / Mathematical Biosciences 222 (2009) 86–91
of the characteristic length of the inhibitory axons. The excitatory fibers are parceled into N groups, with relative densities qn and characteristic lengths k1 n . Still following Nunez, we note the marked disparity in scale between inhibitory axons and the longest excitatory fibers, the latter having a characteristic length on the order of 10 cm, while the former are much more local, being generally confined to the 1 mm size of a cortical macro-column. We consequently treat the activity as sensibly constant over their range in carrying out the inhibitory integral. Using cylindrical coordinates centered at ðx; yÞ leads to an easy calculation.
tial differential form. This allows its evolution in time to be calculated from a much more limited history than would be required by a direct attack on the integral formulation. To that end we form the surface LaPlacian, noting the necessary involvement of derivatives with respect to time when taking the partials with respect to either x or y.
h i ðo=ot þ kE v Þ2 v 2 r2 ðHE GUÞ Z Z ðv =rÞHE ðx0 ; y0 ; t r=v Þ ¼ Bk2E =2p ðo=ot þ kE v Þ S0
0
HI ðx; y; tÞ ¼ ½Q E HE ðx; y; tÞ Q I HI ðx; y; tÞðq pÞ Z Z expfkI rgr dr dh ¼ qI ½Q E HE ðx; y; tÞ Q I HI ðx; y; tÞ;
ð2Þ
S0
HI ðx; y; tÞ ¼ Q E qI HE ðx; y; tÞ=½1 þ Q I qI : This fixes the relationship between excitatory and inhibitory activity, and allows the former to be determined from a single equation:
Z Z N X HE ðx; y; tÞ ¼ Uðx; y; tÞ þ ðqn k2n =2pÞ ½Q E HE ðx0 ; y0 ; t r=v Þ S0
n¼1
0
0
Q I HI ðx0 ; y0 ; t r=v Þ expfkn rg dx dy ¼ Uðx; y; tÞ þ Q E =½1 þ Q I qI
0
HE ðx0 ; y0 ; t r=v Þ expfkn rg dx dy :
ð3Þ
HE ðx; y; tÞ ¼ Uðx; y; tÞ þ Q E =ð1 þ Q I qI Þ " X qn HE ðx; y; tÞ þ ðqE k2E =2pÞ n–E
0 0 HE ðx0 ; y0 ; t r=v Þ expfkE rg dx dy ; S0 X qn HE ðx; y; tÞ½1 þ Q I qI Q E Z Z
n–E
ð4Þ
¼ ð1 þ Q I qI ÞUðx; y; tÞ þ ðQ E qE k2E =2pÞ Z Z 0 0 HE ðx0 ; y0 ; t r=v Þ expfkE rg dx dy ; S0
HE ðx; y; tÞ ¼ GUðx; y; tÞ þ ðBk2E =2pÞ Z Z 0 0 HE ðx0 ; y0 ; t r=v Þ expfkE rg dx dy ; S0
where
," G ð1 þ Q I qI Þ
1 þ Q I qI Q E
," B Q E qE
1 þ Q I qI Q E
X
#
qn
#
qn :
n–E
G will soon be seen to be the gain for static inputs of short wavelength. B is a dimensionless control parameter identified by Nunez. It reflects the ratio of positive to negative feedback, and governs the stability of the system. Its value is thought to change depending on the brain’s state, an understanding that is reinforced by comparing steady-state solutions to EEG spectra as recorded under varying levels of Halothane [2]. Nunez and Katznelson have solved this equation in several geometric configurations. We contemplate a digital simulation, so it is advantageous to render the remaining integral equation into a par-
ð6Þ
The partial differential equation is solved by setting D ¼ HEGU and grouping like terms.
h i ð1 þ BÞo2 =ot 2 þ 2kE v o=ot þ ð1 BÞk2E v 2 v 2 r2 D ¼ BG k2E v 2 o2 ot 2 U:
ð7Þ
We solved this equation digitally, estimating the Laplacian from the nearest neighbors and calculating the new value at each point on the surface from the current distribution and the next most recent value at that point. We used a spherical surface with a circumference of 80 cm, corresponding roughly to the extent of a cortical hemisphere when its convolutions have been expanded. Some intuition for the behavior to be expected can be gained by examining a few limiting cases. Suppose the input is spatially disrffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 tributed in a sinusoidal pattern with wave number k ¼ kx þ ky , and sufficient time is allowed for the cortical activity to come fully to equilibrium. At that point all derivatives with respect to time vanish, and we have
.h i 2 2 HE ðx; yÞ ¼ Uðx; yÞG k2E þ k ð1 BÞk2E þ k : 2
and
n–E
X
h i ðo=ot þ kE v Þ2 v 2 r2 ðHE GUÞ Z Z ½v HE ðx; y; tÞ roHE ðx; y; tÞ=ot ¼ ðBk2E =2pÞðo=ot þ kE v Þ ¼ Bðo=ot þ kE v ÞðkE v o=otÞHE ¼ B k2E v 2 o2 ot 2 HE :
We treat similarly all but the longest excitatory fibers, taking the activity as constant over their range for the purposes of the integration. We use the subscript E to denote the longest group, and break it out of the sum for special treatment.
The expression on the left contains the major features of a damped wave equation, while the 1=r factor causes the integral on the right to converge rapidly. In it we approximate HE by its value at (x, y, t) plus its first derivative in each variable. The linear terms in the spatial derivatives vanish under the symmetric integration, leaving
S0
ðqn k2n =2pÞ 0
S0
ð5Þ
expfkE rg dr dh
N X n¼1
Z Z
0
expfkE rg dx dy :
2 I kI =2
ð8Þ
Short wavelengths ðk k2 Þ are seen to have gain G. Long wavelengths are amplified by a larger factor, and stability under uniform excitation requires that B < 1. The initial portion of the readiness potential of interest here is distributed across the scalp with an effective wavelength approximately equal to the characteristic length of the long excitatory fibers. (Both are on the order of a hemispheric radius.) Therefore a moderate increase in gain is to be expected in the simulation during the priming phase. In fact, however, this probably understates the importance of the effect in real brains. By taking all excitatory fibers except those in the longest group to be very short in Eq. (4), we lost the ability to trace corresponding increases in gain as wavelengths become long with respect to the characteristic lengths of other excitatory groups. If an input is applied suddenly, the response will require some time to adapt. The transients have the form exp{st} where s satisfies the equation
J.W. Hartwell / Mathematical Biosciences 222 (2009) 86–91 2
ð1 þ BÞs2 þ 2kE v s þ ð1 BÞk2E v 2 þ v 2 k ¼ 0; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2 s ¼ kE v =ð1 þ BÞ 1 i k ð1 þ BÞ=k2E B2 :
ð9Þ
For short wavelengths the transients are oscillatory. They decay with time constant pffiffiffiffiffiffiffiffiffiffiffiffis ¼ ð1 þ BÞ=ðkE v Þ, and oscillate at radian frequency x ¼ kv = 1 þ B. Using a characteristic excitatory fiber length of 10 cm and a propagation velocity of 800 cm/s gives a time constant of about 20 ms. Long wavelengths are over-damped and build to their final values without oscillation. In the limiting case of uniform excitation, the slower time constant is s ¼ ð1 þ BÞ=½kE v ð1 BÞ, which again requires that B < 1. In our simulation we used B ¼ 0:95. Since the gain G is purely multiplicative, applied equally to all spatial distributions, we used the limiting value G ¼ 1. This may serve to emphasize the importance of long wavelengths in the initial priming.
References [1] W.G. Walter, R. Cooper, V.J. Aldridge, W.C. McCallum, A.L. Winter, Contingent negative variation: an electrical sign of sensorimotor association and expectancy in the human brain, Nature 203 (1964) 380.
91
[2] A. Givens, Signals of preparation: the contingent negative variation, in: E.D. Bigler (Ed.), Neuroimaging I: Basic Science, Springer, New York, 1996. ISBN: 0306452286. [3] H.H. Kornhuber, L. Deecke, Hirnpotentialänderungen bei Willkürbewegungen und passiven Bewegungen des Menschen: Bereitschaftspotential und reafferente Potentiale, Pflügers Arch 284 (1965) 1. [4] A summary of modern and historic work on the readiness potential, including graphs of its time course at various electrode sites is to be found in the Wikipedia. URL: http://en.wikipedia.org/wiki/Bereitschaftspotential. [5] S. Makeig, K. Gramann, T. Jung, T. Sejnowski, H. Poizner, Linking brain, mind and behavior, Int. J. Psychophysiol. 73 (2009) 95–100. [6] P.L. Nunez, Neocortical Dynamics and Human EEG Rhythms, Oxford University, New York, 1995. ISBN: 0-19-505728-7. [7] P.R. Kennedy, R.A. Bakay, Restoration of neural output from a paralyzed patient by a direct brain connection, Neuroreport 1 (9) (1998) 1707. [8] L.R. Hochberg, M.D. Serruya, G.M. Friehs, J.A. Mukand, M. Saleh, A.H. Caplan, A. Branner, D. Chen, R.D. Penn, J.P. Donoghue, Neuronal ensemble control of prosthetic devices by a human with tetraplegia, Nature 442 (2006) 164. [9] R.D. Katznelson, Normal modes of the brain: neuroanatomical basis and a physiological theoretical model, in: P.L. Nunez (Ed.), The Neurophysics of EEG, Oxford University, New York, 1981. ISBN: 0195027965. [10] P.L. Nunez, Generation of human EEG by a combination of long and short range neocortical interactions, Brain Topogr. 1 (1989) 199–215. [11] P.L. Nunez, R. Srinivasan, Electric Fields of the Brian: The Neurophysics of EEG, second ed., Oxford University, New York, 2006. ISBN: 019505038X. [12] V.K. Jirsa, H. Haken, A derivation of a macroscopic field theory of the brain from the quasi-microscopic neural dynamics, Physica D 99 (1997) 503–526. [13] P.L. Nunez, personal communication. See also Ref. [11], p. 498ff.
Contents lists available at ScienceDirect
Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs
Cortical transients preceding voluntary movement J.W. Hartwell * JW Hartwell & Associates, 3001 Hartwell Pond Drive, Hillsborough, NC 27278, USA
a r t i c l e
i n f o
Article history: Received 21 May 2009 Received in revised form 4 September 2009 Accepted 17 September 2009 Available online 25 September 2009 Keywords: Cortical dynamics Cortical transients Voluntary movement Readiness potential Contingent negative variation
a b s t r a c t The process of initiating a voluntary muscular movement evidently involves a focusing of diffuse brain activity onto a highly specific location in the primary motor cortex. Even the very simple stereotypic movements used to study the ‘contingent negative variation’ and the ‘readiness potential’ begin with EEG indicative of widely distributed brain activity. In natural settings the involvement of diffuse cortical networks is undoubtedly even more important. Eventually, however, activity must coalesce onto specific neurons for the intended movement to ensue. Here we examine that focusing process from a mathematical point of view. Using a digital simulation, we solve the global equations for cortical dynamics and model the flow from diffuse onset to localized spike. From this perspective the interplay between global and local effects is seen as a necessary consequence of a basic cortical architecture which supports wave propagation. Watching the process evolve over time allows us to estimate some characteristic amplitudes and delays. Ó 2009 Elsevier Inc. All rights reserved.
1. Introduction The process of initiating a voluntary muscular movement evidently involves a focusing of diffuse brain activity onto a highly specific location in the primary motor cortex. Even the very simple stereotypic movements used to study the ‘contingent negative variation’ (CNV) [1,2] and the ‘readiness potential’ [3,4], begin with EEG apparently indicative of widely distributed brain dynamics. In natural settings the involvement of diffuse cortical networks is undoubtedly even more important [5]. Eventually, however, activity must coalesce onto specific neurons for the intended bodily movement to ensue. Here we examine that focusing process from a mathematical point of view. Using a digital simulation, we solve the global equations for cortical dynamics [6] and model the flow from diffuse onset to localized spike. From this perspective the interplay between global and local effects is seen as a necessary consequence of the basic cortical architecture which supports wave propagation. An understanding of brain activity preceding normal movements could have important practical implications. As prosthetic devices become more fully functional imitations of the limbs and organs they replace, the interfaces for their control take on added significance. To achieve smooth activation of these complex devices, engineers seek ‘intuitive’ control mechanisms, mimicking the natural flow insofar as they are able and interfacing to the body’s neural system as early as understanding and technology permit. ‘Mind reading’ might summarize the goal. If a user’s inten* Tel.: +1 919 732 7951. E-mail address: [email protected] 0025-5564/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2009.09.003
tions for fine motor movements could be reliably detected, then mechanical devices might become much more useful to him. And if that detection took place at the cortical level, then response times could be shortened and prosthetic devices used even by those with severe injuries to the peripheral nervous system. Advances in this technology might prove especially beneficial to socalled ‘locked-in’ patients suffering (mainly) from ALS and brainstem strokes. Already some success with brain-computer interfaces has been achieved, enabling patients to move a cursor on a computer screen [7] and to control an artificial hand [8]. Regrettably little is known about the specific cortical events that precede voluntary muscle movements. The ‘motor strip’ itself is well mapped, and it is clear that highly specific neurons in it must fire to initiate the contraction of particular skeletal muscles. The necessary cortical activity must have a distinctly localized character, focusing on certain cells and excluding others located nearby. But these same neurons are richly interconnected with others in a structure that supports complex wave activity and precludes their firing as isolated entities. Indeed, EEG patterns recorded preceding voluntary movements are characterized by very diffuse waves at the onset. The mathematical simulation reported here represents an attempt to understand and model how this diffuse beginning leads to the eventual sharply localized result. Watching the process evolve over time allows us to estimate some characteristic amplitudes and delays. Both the ‘contingent negative variation’ (CNV) [1,2] and the ‘readiness potential’ [3,4] are EEG waves marked by a negative voltage spread over a wide region at the top of the head. This voltage grows in magnitude over a period of about half a second prior to a muscular response. In the case of the CNV, larger voltages are
J.W. Hartwell / Mathematical Biosciences 222 (2009) 86–91
associated with more intense involvement in the task and with shorter reaction times [2]. Occasionally the CNV is sufficiently large as to be observed in the raw EEG, but its study usually requires time-locked averaging. The readiness potential is similarly small, but somewhat differently distributed. In particular, it is more pronounced over the hemisphere contra-lateral to the responding hand. Alan Givens writes, ‘About 200 milliseconds after the second stimulus, the CNV ends in a sharp positive-going (10 Hz) wave termed the CNV resolution’ [2]. The readiness potential appears to end less suddenly [4]. In any event, this marks the onset of muscle movement, and concomitant inputs to the adjoining somatosensory cortex contribute to the EEG from this point forward. To see how a sharply localized response can evolve in this richly interconnected system, we have examined the global equations for ‘neocortical dynamics’, used successfully by Nunez [6,10,11], Katznelson [9], Nunez and Srinivasan [11] and others to elucidate widespread cortical rhythms as observed in EEG. These equations are based on a columnar model of the cortex with inhibitory interconnections that are predominantly short-range and with excitatory connections that at their longest span a cortical hemisphere. The density of these long fibers decreases with distance, with a characteristic length of about 10 cm. Their finite conduction velocity gives rise to the propagation of waves in synaptic activity. Both Nunez [10] and Jirsa and Haken [12] have joined the global equations to the contributions of local networks in the generation of EEG. The current work stops short of this goal, considering only global relations. The equations are normally written in integral form, with the activity at a particular cortical column being determined by the subcortical and commissural inputs to that location plus an integral over connections from the surrounding cortex, taking activity at earlier times as required by the intervening distance and the conduction velocity. About 92% of the long axon nerve fibers arriving at a cortical column are from other locations in the same hemisphere [9]. Although the inhibitory connections are often associated with a substantial delay, we follow treatments by Nunez in which their effects, along with those of the short-range excitatory fibers, are considered as ‘instantaneous’ (i.e. as comparable to post-synaptic delays within a cortical column, which are also neglected). The dynamics of the system depend in part on a ratio of overall positive to negative feedback as captured in a single parameter, B, which controls the gain for inputs of relatively long wavelength. We have taken B ¼ 0:95 which preserves stability of the linearized equations. Larger values of B lead to more labile behavior and are considered more realistic by some workers, but they require constant subcortical inhibitory inputs to keep the system operating in a quasi-linear region [13]. In the work reported here, the equations have been transformed from their usual integral representation to an equivalent partial differential form better suited to digital simulation (see Appendix A). They are solved on a closed spherical surface with a circumference of 80 cm, corresponding approximately to a cortical hemisphere ‘inflated’ so as to smooth its fissures and undulations. The conduction velocity of the long fibers is taken to be 800 cm/s. The digital simulation used discrete cortical elements with centers separated by no more than 0.8 mm and time steps of 1.0 ms.
2. Initial priming Our simulation begins with a diffuse, low-level input held constant for half a second. The input has a broad maximum over the eventual focus point in one hemisphere and tapers off to zero on its opposite side. The resulting cortical activity takes on a similar pattern and reaches somewhat larger amplitude. The response rises
87
slowly with a time constant of about 100 ms and levels off with a relative gain of about 1.8. (The absolute gain is considerably greater, but we wish to demonstrate a benefit of priming with long wavelengths. As mentioned in Appendix A, we have reduced to unity a multiplicative constant that affects all wavelengths equally.) Fig. 1 shows the steady input (green) and increasing response (red), now nearing its final amplitude at the peak location. This response represents increased post-synaptic activity in pyramidal neurons across a large region. These cells are aligned perpendicular to the cortical surface, and their collective activity is reflected in a similarly widespread negative voltage on the scalp, corresponding to the initial phases of the CNV and readiness potential. 3. Trigger onset At any time after initial priming, motion can be triggered by changing the input over a disc surrounding the focus point. Here a small negative input (i.e. a reduction from prevailing ambient) is applied over a disc of radius 35 mm as shown in Fig. 2. Responses corresponding to inputs of short wavelength are oscillatory and under-damped. A small wave front begins to propagate in both directions from the discontinuity. The out-bound wave decreases in amplitude as it propagates both because of the overall preponderance of negative feedback and because the wave front spreads over an increasingly larger perimeter. The in-bound wave, however, grows in amplitude as it concentrates over a decreasing perimeter, this effect out-weighing that of the negative feedback. A disc with radius 35 mm is near the optimum for the ratio of positive to negative feedback used in the simulation ðB ¼ 0:95Þ. For larger discs the inbound wave initially loses amplitude, while smaller ones would benefit from a greater circumference of initial activity. 4. A big bounce The wave of reduced activity grows in amplitude as it propagates toward the disc center as shown in Figs. 3 and 4 (red1 traces) where it converges to a trough after 62 ms. The locally reduced activity is reflected in a momentarily positive scalp voltage. The system will rebound in a few milliseconds. All input was removed 5 ms after the trigger pulse was applied. (Green traces in Figs. 3 and 4 show zero.) The consequent second wave of opposite polarity is also converging and will reinforce the oscillatory maximum. 5. Focusing pulse As the effects of the trigger pulse and its removal converge, a small focusing pulse is applied to select the particular motor neurons necessary for the desired muscle movement. Being of very short spatial extent, the effect of this pulse is almost immediate. It combines with the waves initiated earlier to drive the activity to a large value at a highly specific location. This is shown in Fig. 5 where the overall pulse amplification is more than 13:1. The narrow width of the resulting activity places it beyond the validity of the global wave equations. A more detailed accounting would require considering interactions with the local motor network [10]. Although the trigger and focusing pulses are of low amplitude, their relative timing is crucial. For the 35 mm trigger disc used in our simulation, the focusing pulse must follow the trigger by 69 ms to achieve its greatest effect. This constitutes an appreciable fraction of the very fast reaction times often recorded in experiments in which the CNV is prominently seen. 1 For interpretation of the references to color in Fig. 4, the reader is referred to the web version of this paper.
88
J.W. Hartwell / Mathematical Biosciences 222 (2009) 86–91
Fig. 1. Simulation through 503 ms. Initial priming and trigger at focus. Green trace shows input, constant for half a second, with sudden trigger onset. Red trace shows gradually increasing cortical activity now nearing its final amplitude. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)
Fig. 2. Cortical slice at 503 ms. Green trace shows input with trigger extending over a disc of radius 35 mm. Red trace shows activity. Small waves begin to propagate away from the discontinuity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)
Fig. 3. Simulation at 532 ms. Inbound wave increases in amplitude; outbound wave decreases. All input was removed at 505 ms (green trace shows zero). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)
Fig. 4. Simulation at 562 ms. Inbound activity has converged to a deep trough. Input remains zero.
Fig. 5. Simulation at 569 ms. Extremely narrow pulse of activity follows a small focusing input. Green trace shows spatial distribution of input and red the resultant activity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)
89
J.W. Hartwell / Mathematical Biosciences 222 (2009) 86–91
Fig. 6. Recovery at focus. Response to trigger and focusing pulse decays with rapid oscillations. Initial priming disappears with a long time constant.
Fig. 7. Simulation at 1000 ms. Only a portion of the initial priming remains.
6. Recovery The duration of the focusing pulse is not critical. The spike in cortical activity is extremely narrow in time as well as in space, persisting for only a few milliseconds, less than the refraction time of the motor neurons at its focus. The dynamics during this short interval must in fact be governed more by the local motor network than by the global wave theory used here. Recovery from the trigger and focusing pulses again involves global wave propagation, however. The global recovery requires about 80 ms during which the decay shows oscillations at approximately 150 Hz (too high to be seen in normal EEG). See Fig. 6. This, together with the 69 ms lead time required by the trigger pulse, gives a duration of around 150 ms for the firing process. Because the system is operating in a quasi-linear mode, there is no corresponding time delay required between successive reenactments. A new sequence could be initiated at any point. The long wavelength initial priming disappears in about half a second after its input is removed, decaying with the same long time constant that governed its build up. The state after 1000 ms is shown in Fig. 7. 7. Summary The initial diffuse wave provides a platform of coordinated activity that allows more sharply defined inputs to be magnified in effect. The mechanism benefits from the high gain associated with static inputs of long wavelength. The transition from this global preparation to the highly discrete final result utilizes converging waves initiated by a small but more localized trigger pulse. The in-bound pattern may be visualized as similar to what would result from throwing a rock into a pond and then playing the movie of the resulting waves backwards. Finally, the precise focus is determined by local interactions, represented here by a small short pulse whose effect is greatly amplified by the waves preceding it. The complete sequence utilizes the dynamic character of the global cortical interconnections somewhat like a gymnast exploits wave motion in the elastic medium of a trampoline, using preparation and precise timing to achieve notable effect. 8. Some consequences The good news for prosthetics engineers in this simulation is that transients associated with the trigger for a voluntary movement be-
gin about 69 ms prior the firing of the specific motor neurons involved. This gives an opportunity to detect and respond. The bad news is that the waves begin small and distributed over an area that covers much of the motor strip. Most easily detected is the widespread excitatory priming activity which seems only to facilitate motor output generally. The trigger, which lends the first hint of specificity, consists in a broad disc which gives rise to a slowly collapsing annulus. The specific point of focus stands out only very late in the sequence. Sensors and circuitry to detect such a developing pattern might have to resemble the cortical wiring itself. Acknowledgments Thanks to Paul Nunez and the reviewers for detailed and helpful commentary. Appendix A Following Nunez’ ‘Solutions to the Equations for Global Neocortical Dynamics’ as presented in Appendix A to Ref. [6], we adopt similar notation and write the linearized (small signal) equations for excitatory and inhibitory activity at cortical location ðx; yÞ and at time t as
HE ðx; y; tÞ ¼ Uðx; y; tÞ þ
N X
ðq00 k200 =2pÞ
n¼1
Z Z S0
½Q E HE ðx0 ; y0 ; t r=v Þ 0
0
0
0
Q I HI ðx0 ; y0 ; t r=v Þ expfkn rg dx dy ; Z Z HI ðx; y; tÞ ¼ ðqI k2I =2pÞ ½Q E HE ðx0 ; y0 ; t r=v Þ S0
Q I HI ðx0 ; y0 ; t r=v Þ expfkI rg dx dy ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ðx x0 Þ2 þ ðy y0 Þ2 : ð1Þ Here Uðx; y; tÞ denotes excitatory inputs arriving from the opposite hemisphere or from sub-cortical structures at the same time and place. The integrals sum the contributions of cortico-cortical fibers from the entire surface of one hemisphere, weighting each region by the density of fibers from it (which decrease exponentially with distance), and using activity from the earlier time t r=v to account for the finite neural propagation velocity v. The coefficients Q E and Q I give the sensitivity for small increases or decreases in excitatory and inhibitory inputs, respectively. kI is the reciprocal
90
J.W. Hartwell / Mathematical Biosciences 222 (2009) 86–91
of the characteristic length of the inhibitory axons. The excitatory fibers are parceled into N groups, with relative densities qn and characteristic lengths k1 n . Still following Nunez, we note the marked disparity in scale between inhibitory axons and the longest excitatory fibers, the latter having a characteristic length on the order of 10 cm, while the former are much more local, being generally confined to the 1 mm size of a cortical macro-column. We consequently treat the activity as sensibly constant over their range in carrying out the inhibitory integral. Using cylindrical coordinates centered at ðx; yÞ leads to an easy calculation.
tial differential form. This allows its evolution in time to be calculated from a much more limited history than would be required by a direct attack on the integral formulation. To that end we form the surface LaPlacian, noting the necessary involvement of derivatives with respect to time when taking the partials with respect to either x or y.
h i ðo=ot þ kE v Þ2 v 2 r2 ðHE GUÞ Z Z ðv =rÞHE ðx0 ; y0 ; t r=v Þ ¼ Bk2E =2p ðo=ot þ kE v Þ S0
0
HI ðx; y; tÞ ¼ ½Q E HE ðx; y; tÞ Q I HI ðx; y; tÞðq pÞ Z Z expfkI rgr dr dh ¼ qI ½Q E HE ðx; y; tÞ Q I HI ðx; y; tÞ;
ð2Þ
S0
HI ðx; y; tÞ ¼ Q E qI HE ðx; y; tÞ=½1 þ Q I qI : This fixes the relationship between excitatory and inhibitory activity, and allows the former to be determined from a single equation:
Z Z N X HE ðx; y; tÞ ¼ Uðx; y; tÞ þ ðqn k2n =2pÞ ½Q E HE ðx0 ; y0 ; t r=v Þ S0
n¼1
0
0
Q I HI ðx0 ; y0 ; t r=v Þ expfkn rg dx dy ¼ Uðx; y; tÞ þ Q E =½1 þ Q I qI
0
HE ðx0 ; y0 ; t r=v Þ expfkn rg dx dy :
ð3Þ
HE ðx; y; tÞ ¼ Uðx; y; tÞ þ Q E =ð1 þ Q I qI Þ " X qn HE ðx; y; tÞ þ ðqE k2E =2pÞ n–E
0 0 HE ðx0 ; y0 ; t r=v Þ expfkE rg dx dy ; S0 X qn HE ðx; y; tÞ½1 þ Q I qI Q E Z Z
n–E
ð4Þ
¼ ð1 þ Q I qI ÞUðx; y; tÞ þ ðQ E qE k2E =2pÞ Z Z 0 0 HE ðx0 ; y0 ; t r=v Þ expfkE rg dx dy ; S0
HE ðx; y; tÞ ¼ GUðx; y; tÞ þ ðBk2E =2pÞ Z Z 0 0 HE ðx0 ; y0 ; t r=v Þ expfkE rg dx dy ; S0
where
," G ð1 þ Q I qI Þ
1 þ Q I qI Q E
," B Q E qE
1 þ Q I qI Q E
X
#
qn
#
qn :
n–E
G will soon be seen to be the gain for static inputs of short wavelength. B is a dimensionless control parameter identified by Nunez. It reflects the ratio of positive to negative feedback, and governs the stability of the system. Its value is thought to change depending on the brain’s state, an understanding that is reinforced by comparing steady-state solutions to EEG spectra as recorded under varying levels of Halothane [2]. Nunez and Katznelson have solved this equation in several geometric configurations. We contemplate a digital simulation, so it is advantageous to render the remaining integral equation into a par-
ð6Þ
The partial differential equation is solved by setting D ¼ HEGU and grouping like terms.
h i ð1 þ BÞo2 =ot 2 þ 2kE v o=ot þ ð1 BÞk2E v 2 v 2 r2 D ¼ BG k2E v 2 o2 ot 2 U:
ð7Þ
We solved this equation digitally, estimating the Laplacian from the nearest neighbors and calculating the new value at each point on the surface from the current distribution and the next most recent value at that point. We used a spherical surface with a circumference of 80 cm, corresponding roughly to the extent of a cortical hemisphere when its convolutions have been expanded. Some intuition for the behavior to be expected can be gained by examining a few limiting cases. Suppose the input is spatially disrffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 tributed in a sinusoidal pattern with wave number k ¼ kx þ ky , and sufficient time is allowed for the cortical activity to come fully to equilibrium. At that point all derivatives with respect to time vanish, and we have
.h i 2 2 HE ðx; yÞ ¼ Uðx; yÞG k2E þ k ð1 BÞk2E þ k : 2
and
n–E
X
h i ðo=ot þ kE v Þ2 v 2 r2 ðHE GUÞ Z Z ½v HE ðx; y; tÞ roHE ðx; y; tÞ=ot ¼ ðBk2E =2pÞðo=ot þ kE v Þ ¼ Bðo=ot þ kE v ÞðkE v o=otÞHE ¼ B k2E v 2 o2 ot 2 HE :
We treat similarly all but the longest excitatory fibers, taking the activity as constant over their range for the purposes of the integration. We use the subscript E to denote the longest group, and break it out of the sum for special treatment.
The expression on the left contains the major features of a damped wave equation, while the 1=r factor causes the integral on the right to converge rapidly. In it we approximate HE by its value at (x, y, t) plus its first derivative in each variable. The linear terms in the spatial derivatives vanish under the symmetric integration, leaving
S0
ðqn k2n =2pÞ 0
S0
ð5Þ
expfkE rg dr dh
N X n¼1
Z Z
0
expfkE rg dx dy :
2 I kI =2
ð8Þ
Short wavelengths ðk k2 Þ are seen to have gain G. Long wavelengths are amplified by a larger factor, and stability under uniform excitation requires that B < 1. The initial portion of the readiness potential of interest here is distributed across the scalp with an effective wavelength approximately equal to the characteristic length of the long excitatory fibers. (Both are on the order of a hemispheric radius.) Therefore a moderate increase in gain is to be expected in the simulation during the priming phase. In fact, however, this probably understates the importance of the effect in real brains. By taking all excitatory fibers except those in the longest group to be very short in Eq. (4), we lost the ability to trace corresponding increases in gain as wavelengths become long with respect to the characteristic lengths of other excitatory groups. If an input is applied suddenly, the response will require some time to adapt. The transients have the form exp{st} where s satisfies the equation
J.W. Hartwell / Mathematical Biosciences 222 (2009) 86–91 2
ð1 þ BÞs2 þ 2kE v s þ ð1 BÞk2E v 2 þ v 2 k ¼ 0; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2 s ¼ kE v =ð1 þ BÞ 1 i k ð1 þ BÞ=k2E B2 :
ð9Þ
For short wavelengths the transients are oscillatory. They decay with time constant pffiffiffiffiffiffiffiffiffiffiffiffis ¼ ð1 þ BÞ=ðkE v Þ, and oscillate at radian frequency x ¼ kv = 1 þ B. Using a characteristic excitatory fiber length of 10 cm and a propagation velocity of 800 cm/s gives a time constant of about 20 ms. Long wavelengths are over-damped and build to their final values without oscillation. In the limiting case of uniform excitation, the slower time constant is s ¼ ð1 þ BÞ=½kE v ð1 BÞ, which again requires that B < 1. In our simulation we used B ¼ 0:95. Since the gain G is purely multiplicative, applied equally to all spatial distributions, we used the limiting value G ¼ 1. This may serve to emphasize the importance of long wavelengths in the initial priming.
References [1] W.G. Walter, R. Cooper, V.J. Aldridge, W.C. McCallum, A.L. Winter, Contingent negative variation: an electrical sign of sensorimotor association and expectancy in the human brain, Nature 203 (1964) 380.
91
[2] A. Givens, Signals of preparation: the contingent negative variation, in: E.D. Bigler (Ed.), Neuroimaging I: Basic Science, Springer, New York, 1996. ISBN: 0306452286. [3] H.H. Kornhuber, L. Deecke, Hirnpotentialänderungen bei Willkürbewegungen und passiven Bewegungen des Menschen: Bereitschaftspotential und reafferente Potentiale, Pflügers Arch 284 (1965) 1. [4] A summary of modern and historic work on the readiness potential, including graphs of its time course at various electrode sites is to be found in the Wikipedia. URL: http://en.wikipedia.org/wiki/Bereitschaftspotential. [5] S. Makeig, K. Gramann, T. Jung, T. Sejnowski, H. Poizner, Linking brain, mind and behavior, Int. J. Psychophysiol. 73 (2009) 95–100. [6] P.L. Nunez, Neocortical Dynamics and Human EEG Rhythms, Oxford University, New York, 1995. ISBN: 0-19-505728-7. [7] P.R. Kennedy, R.A. Bakay, Restoration of neural output from a paralyzed patient by a direct brain connection, Neuroreport 1 (9) (1998) 1707. [8] L.R. Hochberg, M.D. Serruya, G.M. Friehs, J.A. Mukand, M. Saleh, A.H. Caplan, A. Branner, D. Chen, R.D. Penn, J.P. Donoghue, Neuronal ensemble control of prosthetic devices by a human with tetraplegia, Nature 442 (2006) 164. [9] R.D. Katznelson, Normal modes of the brain: neuroanatomical basis and a physiological theoretical model, in: P.L. Nunez (Ed.), The Neurophysics of EEG, Oxford University, New York, 1981. ISBN: 0195027965. [10] P.L. Nunez, Generation of human EEG by a combination of long and short range neocortical interactions, Brain Topogr. 1 (1989) 199–215. [11] P.L. Nunez, R. Srinivasan, Electric Fields of the Brian: The Neurophysics of EEG, second ed., Oxford University, New York, 2006. ISBN: 019505038X. [12] V.K. Jirsa, H. Haken, A derivation of a macroscopic field theory of the brain from the quasi-microscopic neural dynamics, Physica D 99 (1997) 503–526. [13] P.L. Nunez, personal communication. See also Ref. [11], p. 498ff.