Using a virtual cortical module implementing a neural field model to modulate brain rhythms in Parkinson’s disease

We propose a new method for selective modulation of cortical rhythms based on neural field theory, in which the activity of a cortical area is extensively monitored using a two-dimensional microelectrode array. The example of Parkinson’s disease illustrates the proposed method, in which a neural field model is assumed to accurately describe experimentally recorded activity. In addition, we propose a new closed-loop stimulation signal that is both space- and time- dependent. This method is especially designed to specifically modulate a targeted brain rhythm, without interfering with other rhythms. A new class of neuroprosthetic devices is also proposed, in which the multielectrode array is seen as an artificial neural network interacting with biological tissue. Such a bio-inspired approach may provide a solution to optimize interactions between the stimulation device and the cortex aiming to attenuate or augment specific cortical rhythms. The next step will be to validate this new approach experimentally in patients with Parkinson’s disease.

Introduction

The development of neuroprosthetic devices designed to restore brain function is currently becoming the focus of intense research (see for example the special issue of Expert Rev. Med Devices, in 2007; Neuroscience and Biobehavioral Reviews in 2008, focusing on Deep Brain Stimulation, DBS; and special issue on brain–computer interfaces of Neural Networks in 2009). In DBS, a stimulating device composed of a single macroelectrode (on each side of the brain) and a stimulator is used to alleviate a variety of neurological symptoms. Electrodes are typically implanted in a subcortical structure, e.g., the subthalamic nucleus in patients with Parkinson’s disease (PD). Each electrode is then connected to a stimulator (similar to a cardiac pacemaker) implanted in the chest wall of the patient. In this class of devices, once stimulation parameters have been fixed, the same stimulation signal is applied continuously as long as the stimulator is turned on. A single location of the motor network (i.e., the cortico-basal ganglia-thalamo-cortical loop) is stimulated, independently of the spatial distribution of network activity. Current DBS limitations include the need to replace batteries regularly, and in some cases the presence of side effects such as depression, speech difficulties, and weight gain (Guehl et al., 2007 ). Despite these limitations, DBS has become a standard procedure for the relief of PD symptoms (currently more than 40,000 patients receive DBS in the world) and to a lesser extent for other disorders such as pain (1500–2000 patients), non-parkinsonian tremor (500–1000), and dystonia (300–500) (Pereira et al., 2007 ). In addition DBS is used to treat cluster headaches, epilepsy, obsessive compulsive disorder, depression, and Tourette syndrome. The main structures stimulated belong to the basal ganglia, thalamus, internal capsule, zona incerta, and hypothalamus. Recently, a superficial structure, i.e., the cortex has been added as a new target (Canavero et al., 2002 ; Pagni et al., 2003 ). The physiological mechanisms underlying the improvement of symptoms by stimulating different neural structures are still being actively investigated. Briefly, it is admitted that DBS affects the cortico-basal ganglia-thalamo-cortical network, which is composed of a series of bidirectionally interconnected structures. DBS likely induces functional changes in the whole network as a consequence of a complex cascade of dynamic orthodromic and antidromic effects (Priori and Lefaucheur, 2007 ; Modolo and Beuter, 2009a ).

According to Pereira et al., (2007) the next generation of stimulators is currently being developed. Ideally, these stimulators should take into consideration the scale of the motor network (neither microscopic: at the level of the individual neuron; nor macroscopic: at the level of behavior; but rather mesoscopic: at the level of interacting neural populations). Furthermore, these new stimulators should incorporate the dynamical characteristics of the network such as feedback loops, time delays or synchronization. Finally, according to the principles of individualized medicine stating that therapy should be tailored to each patient’s physiology and needs, these new devices should be adaptable to each patient’s specific needs and focus on customized therapeutic stimulation patterns. This last aspect implies taking into account both spatial and temporal properties of the patient’s cortical activity.

The first aim of this paper is to propose a stimulation pattern that is adaptive both in space and time. Such a stimulation pattern can be derived from the endogenous dynamics of cortical tissue and applied to the motor cortex via a matrix of electrodes. In this paper, we consider PD and focus on the alleviation of tremor, one cardinal symptom of this disease also characterized by rigidity and bradykinesia. Normal physiological tremor is defined as an involuntary oscillation of the extremities at about 8–12 Hz (Elble and Koller, 1990 ). Parkinsonian tremor is slower (4–6 Hz) than normal physiological tremor (Rodriguez-Oroz et al., 2001 ) and its amplitude is augmented and tends to fluctuate over time. It is correlated with frequencies recorded in various cortical and subcortical structures at about 10 Hz (Timmermann et al., 2003 , 2007 ; Timmermann and Fink, 2009 ). Thus, the second aim of this paper is to propose a neuroprosthetic device used for selective modulation of cortical rhythms and based on an extension of the method proposed in the first part of the paper.

The proposed new class of neuroprosthetic devices is an extension of the method proposed in the main part of the paper for selective modulation of cortical rhythms. It is a natural extension of the adaptive feedback idea to see the multielectrode array as an artificial neural network interacting with biological tissue. Such a bio-inspired approach may provide a mean to fine-tune interactions between the stimulation device and the cortex in order to attenuate parkinsonian aspects of cortical activity patterns that cannot be identified by a simple frequency band.

Materials and Methods

The principle of closed-loop stimulation for Parkinson’s disease is not a new idea. Indeed, it has already been studied by different research teams, using methods from various backgrounds (control theory, phase oscillator models). One of the first authors that put forward a model of closed-loop stimulation for Parkinson’s disease is Tass (1999) , who developed with his group methods designed to counteract pathological neuronal dynamics using a mean-field recording (i.e., a mean measure of neuronal activity with a macroelectrode for instance) as information to trigger and compute the feedback. Among these techniques, we mention “non-linear delayed feedback”, where the measured mean-field Z is delayed in time, a non-linear function is applied to this value, and the resultant signal is sent back to the network (Popovych et al., 2008 ). This method has been validated using networks of coupled phase oscillators, and shows that this technique is effective in reducing the abnormal amount of synchronization in the network. Another method developed by this group aims to modulate synaptic plasticity in the neuronal target, in order to “unlearn” pathological activity to the network (Tass and Majtanik, 2006 ). Using a model of phase oscillators with coupling weights modulated by a spike-timing dependent plasticity (STDP) rule, they show that pathological activity is suppressed even after the stimulation is stopped. This approach is promising, since it opens the way to therapeutic applications in which only a transient stimulation is applied, and lasting clinical benefits can be obtained. However, this control is applied with a single electrode, so the neuronal target is not optimally controlled in space. Furthermore, the authors did not investigate the potential perturbation of “physiological” oscillations by the feedback signal. Finally, these methods are still to be tested experimentally.

Another relevant work is the modeling work by Santaniello et al. (2008) , who propose a closed-loop feedback method for DBS, based on control theory (model-based minimum variance law). Interestingly, even if this method is tested on a network of bursting neurons, the method does not require any particular knowledge on the detailed physiological characteristics of the neuronal target, which is convenient. One drawback is that this method appears to modify significantly normal physiological rhythms of neural activity (see Figure 7 in Santaniello et al., 2008 ).

An interesting approach was proposed by Feng et al. (2007) , who used a genetic algorithm (GA) to obtain a model-free stimulation waveform aiming to suppress pathological bursting in the subthalamic nucleus (STN). They applied their method to the well-known model of Rubin and Terman (2004) , and showed that their GA was able to effectively restore a more physiological activity pattern to STN cells. Among effective waveforms, a stochastic train of pulses, in which the mean frequency was the same as classical DBS but where the period was computed stochastically, suppressed abnormal bursting activity effectively. The proposal to use a GA to optimize DBS is undoubtedly interesting, but the effectiveness of a stochastic high-frequency train of pulses appears to contradict experimental results in patients with PD. Indeed, in a clinical paper (Birdno et al., 2008 ), it was shown that the more stochastic the stimulation pattern was, the worse was the clinical improvement. This inconsistency could be caused by possible limitations of the Rubin and Terman model. Another drawback of all of the above proposed methods is that, like the existing DBS procedure, it requires access to deep brain structures, whereas cortical stimulation would not. In this paper, we propose a different approach, in that we aim at a minimal perturbation of physiological brain oscillations, but also at a control both in space and time of neuronal activity. These two constraints represent important advances for the development of clinical applications.

Continuous neural field theory provides a realistic representation of coarse-scale dynamics of cortical tissue in which non-local interactions occur with time delays and a variety of external inputs (Wilson and Cowan, 1973 ; Amari, 1977 ). The basis of neural field models comes from the idea that cortical tissue can be considered as a two-dimensional sheet of tissue (with negligible thickness), and that this tissue can be seen as a continuum rather than discrete elements (Wilson and Cowan, 1973 ; Amari, 1977 ). This is justified by the high density of neurons in cortical tissue (on the order of 100,000 neurons per square millimeter of cortex). The “elements” or “basic units” in question are neural masses, i.e., aggregates of neurons containing excitatory and inhibitory neurons. The neural mass may be a convenient way to model the cortical macrocolumn, a motif of neural pattern repeated throughout the cortex (Horton and Adams, 2005 ). Indeed, since cortical tissue appears composed of synaptically coupled cortical columns, a neural field is a continuous entity composed of coupled neural masses with a connectivity function W(d) depending on the distance, d, between two neural masses. Summation is performed over the neural field, including the spike propagation delay, τ, depending on the conduction speed, c, of the fiber and d. The general form of the function W(d) is a sum of two exponential functions, as explained in Atay and Hutt (2005) and in the Appendix, and depending on the choice of parameters (spatial range of fibers and strength of excitatory/inhibitory connections), the connectivity function may give rise to four different connectivity patterns: (1) locally excitatory, (2) locally inhibitory, (3) locally excitatory and laterally inhibitory, and (4) locally inhibitory and laterally excitatory. It is assumed that connectivity is invariant by translation, that is to say, the connectivity function is the same at each spatial position. Thus, the value of the neural field reflects the mean depolarization of cell membranes. Overall, neural fields appear especially appropriate to evaluate new stimulation patterns aiming to achieve selective effects on cortical activity in predetermined frequency bands.

Using a neural field model as a quantitative framework to describe cortical dynamics in a physiological or pathological situation appears reasonable, since neural field models have proven their descriptive, and more importantly their predictive value in the understanding of cortical dynamics. Indeed, neural field models have been used to understand the formation of visual patterns during hallucinations (Ermentrout and Cowan, 1979 ), have predicted the existence of spiral waves in disinhibited cortex (Huang et al., 2004 ) and of traveling waves (Bressloff, 2001 ; Rubino et al., 2006) . Such a successful description and prediction of phenomena observed experimentally suggests that neural field models can provide a reasonable description of cortical dynamics as measured by a microelectrode array (MEA).

Let us consider a two-dimensional sheet of cortical tissue [which is a good approximation since the thickness of the cortex, about 2.5 mm, is small in relation to its surface area which is about 2,000 cm²; or even compared to the surface area of the primary motor cortex (M1) considered below, which is approximately 135 cm² (Gredal et al., 2000 )]. This cortical tissue is subject to an activity termed as “pathological”, i.e., an activity with defined frequencies and amplitude associated with the presence of symptoms. Of course, normal “physiological” activity is also present. Our assumption that it is possible to separate unambiguously “physiological” and “pathological” activity of brain networks needs to be taken with caution. However, in the particular case of Parkinson’s disease, a significant cortico-muscular coupling has been demonstrated and rest tremor appears to be generated in an extensive cerebral oscillatory network involving cortical and subcortical structures. These structures show extensive coupling among each other with activity at around 8–12 Hz which corresponds to double the tremor frequency (Schnitzler et al., 2006 ). Thus, it is possible to identify with confidence the “pathological” components of neural activity. However, for other neurological disorders, “pathological” activity could be more difficult to identify, although in the case of abnormal postural tremor caused by hepatic encephalopathy, similar changes in the thalamo-cortical network have been observed by Schnitzler et al. (2006) . Overall these results suggest that changes in coupling at different frequencies between structures of this network may represent a common pathophysiological mechanism of tremor disorders. In addition, in healthy humans, the mu rhythm (centered on 10 Hz, ranging from 7 to 14 Hz also termed alpha rhythm elsewhere in the cortex) is noticeable at rest, and disappears before the onset of movement and reappears after the termination of movement (Hatospoulos, 2009 ). One putative role of this rhythm that has been proposed is visuo-motor integration, since this rhythm is modulated either by a motor action or observation of a motor action (Pineda, 2005 ). However, it appears that this rhythm is exaggerated in Parkinson’s disease as measured by MEG (Timmermann et al., 2003 ). Consequently, even if we present results aiming at suppressing as much as possible this mu rhythm, a more realistic view would be to aim not to completely suppress the 10 Hz (in the mu range) rhythm, but rather to decrease it to physiological levels. Attenuating the amplitude of the 10 Hz component by a factor N is straightforward from the formalism presented here, as mentioned below.

In the following, the sheet of cortical tissue is termed the cortical target. Let ψ(r,t) be the neural field value at position r and at time t. We assume that the measured cortical activity is accurately described by the following neural field equation:

where ψ(r, t) is the “potential” of the neural field at position r at time t. More precisely, the potential of the neural field is an “effective” potential, i.e., the deviation of the measured potential, V, with respect to the rest potential, V₀, ψ(r,t) = V(r,t) − V₀. In Eq. 1, k₁ is the time constant of synapses (chosen as 1 ms) with first-order kinetics, γ is the synaptic coupling strength, is the distance-dependent connectivity function, S is a sigmoid function converting the potential into a firing rate, is the propagation speed of spikes along axons and is the sum of external inputs. Finally, is the spatial domain of the cortical target.

Then, we assume that a MEA, part of a neuroprosthetic device, monitors the activity of the two-dimensional cortical target in space and time. We have restricted our focus to modeling the stimulation of the primary motor cortex because it would allow a treatment for tremor in Parkinson’s disease that is much less invasive than current DBS, as well as making adaptive (“on demand”) stimulation more feasible. However, our method is not inherently restricted to the case of cortical targets, but is also applicable, in principle, to stimulation of deeper structures. A crucial aspect of the solution is that the MEA can be used not only to record, but also to stimulate the cortical target. Thus, the MEA provides a mapping of ψ(r,t), the field evolution over time throughout the cortical surface, since the MEA records neuronal activity at the sub-millimeter scale, which is also the scale of description of neural field theory. Figure 1 shows a schematic representation of MEA recording the state of the cortical target.

FIGURE 1

Figure 1. Schematic representation of the cortical target and the MEA. The MEA provides a spatial mapping of neural activity in the cortical target compatible with the spatial scale of description of the neural field model (mesoscopic scale). Consequently, it is possible to exploit the neural field theory to process multi-site recordings provided by the MEA, and to combine these recordings into multiple stimulation signals aiming to modulate neuronal activity.

It should be noted that, for clinical applications, neural field equation parameters will have to be adjusted to fit experimental data. Below, we introduce a closed-loop stimulation term α(r,t) which is an unknown function of the potential ψ(r,t) throughout the spatial domain Ω of the network, but also a function of past values of the potential due to finite conduction delays. The closed-loop stimulation term, or control term, plays the role of a control input to the cortical target. One of the goals is to selectively attenuate, or augment, a predetermined frequency band (Modolo et al., 2010 ). In the following, the case of attenuation is presented, but the case of augmentation is also briefly mentioned afterwards.

Let [ω′,ω″] be the “pathological” frequency band, where ω = 2πf and f is the frequency. For instance, one may consider ω′ = 2π × 9 rad/s and ω″ = 2π × 11 rad/s that includes the 10 Hz rhythm correlated with rest tremor in PD (Timmermann et al., 2003 ), as illustrated in Figure 2 .

FIGURE 2

Figure 2. Localization, frequency analysis and coherence of primary motor cortex (M1) activity with tremor. Magnetoencephalographic (MEG) and electromyographic (EMG) data acquired in a PD patient with right hand tremor. The arrow shows where the MEG signal is recorded. The power spectrum of M1 MEG recording presents a strong peak centered around 10 Hz (middle). Cerebro-muscular (CM) coherence between M1 activity as measured by MEG and EMG of the Extensor Digitorum Communis (EDC) has a maximum around 10 Hz (Fig. modified from Timmermann et al., 2003 , with permission).

Let us write the neural field, ψ(r,t), under the form of a wave packet (Modolo et al., 2010 ) expressed as the summation of “physiological” and “pathological” frequencies:

where is the contribution of undesirable frequencies in intrinsic neural field dynamics. Here, is the Fourier Transform of ψ(r,t) in space and time. Furthermore, let I(r,t) = κ(r,t) + ζ(r,t), where is the “parkinsonian” contribution of the input, i.e., the frequency component causing the “pathological” frequency band centered on 10 Hz, be the input to the cortical target. By including this decomposition of ψ(r,t) and I(r,t) in the main neural field Eq. 1, and by linearizing the sigmoid function, one obtains (Modolo et al., 2010 ):

where the term represents the total input to cells in the “pathological” frequency band. By doing so, we assume that neural field dynamics occur mainly in the linear range. Using perturbation theory, it is possible to improve this approximation by including non-linear contributions, but this is not covered in this paper.

To suppress a “pathological” frequency band (one may also attenuate by a factor N, which is straightforward to derive from the following), the closed-loop stimulation can be expressed as:

which controls neural activity in the “pathological” band only where and when this “pathological” frequency component exceeds a predetermined threshold. When this stimulation term is added to the neural field Eq. 3 the field evolves according to:

In summary, at each position, r, of the cortical target, the recorded membrane potential of each individual microelectrode of the MEA is used in the neural field Eq. 1 to compute the total input of the network, i.e., the sum of the local synaptic activity within the cortical target (corresponding to the integral term of the right-hand side of the neural field equation) and of the inputs from distant brain areas. This estimated total input of the network for the last 100 ms is then filtered using a 9–11 Hz bandpass filter (centered on the 10 Hz target frequency). Then, we use this filtered input to compute the closed-loop stimulation signal, according to Eq. 4. In physical terms, the closed-loop stimulation signal prevents rhythmic depolarization of cell membranes by filtering their synaptic inputs in the “pathological” frequency band. Indeed, the closed-loop stimulation is computed by recording the mean membrane potential at a given point of the network, and using this to control network activity, i.e., by attenuating the “pathological” frequency band.

As illustrated in Figure 1 , the MEA provides the values of the “pathological” frequency component (as detailed below in the Results section) after a simple signal processing consisting in amplifying (since measured signals are small, on the order of the μV), filtering (to keep only the part of the signal corresponding to undesirable frequencies) and analyzing frequency (to compute the power spectrum), since the neural field is mapped in space. Overall, this stimulation signal affects not only individual dynamics of neural masses, but also their coupling. Furthermore, we have shown in a previous study that, in a limit case of our stimulation term (linear delayed feedback, proposed by Rosenblum and Pikovsky, 2004 ), the propagation of cortical waves at predetermined frequencies was strongly disrupted by the stimulation, depending on the delay of the feedback loop (Modolo et al., 2009b ). Figure 3 summarizes the principle of our closed-loop stimulation method.

FIGURE 3

Figure 3. Principle of the proposed closed-loop stimulation method. Local Field Potentials (LFPs) are recorded by micro-electrodes of the high-density (on the order of 25 electrodes per square millimeter, similar to the Utah Array) MEA. These signals are independently amplified, filtered and analyzed to detect the presence of undesirable frequency components above a predetermined threshold. Stimulation is provided to the cortical target only at the spatial position(s) where the threshold is crossed.

In this previous study (Modolo et al., 2010 ), we performed adaptive control of a network of spiking neurons (900 neurons), showing its effectiveness at selectively attenuating bursting activity at a frequency of 10 Hz. In the current work, we consider a closed-loop feedback aiming to control pathological activity in time, as before, but also in space, over a spatially distributed field of neural masses. Such spatio-temporal control aiming at a narrow frequency band is an advance in the control of pathological neural activity and is relevant for clinical use, since cortical activity is not necessarily homogeneous in space in a “real” brain. The use of a neural field model instead of the network of spiking neurons previously developed is consistent with the temporal and spatial resolution capabilities of the MEA (millisecond and sub-millimetric, respectively). We illustrate below our method of spatio-temporal control of network dynamics by closed-loop stimulation in a large cortical network, with a time- and space-dependent stimulation signal.

Results

We used the two-dimensional neural field Eq. 1 to describe the spatio-temporal dynamics of a network of cortical neurons (Eq. 1) of 1 cm² in size (1 cm × 1 cm), with a spatial resolution of 0.1 mm in x and y directions. The cortical network was simulated during 1 s (500 ms without closed-loop stimulation, 500 ms with closed-loop stimulation) with a 1-ms time step. 10 Hz activity was induced in the neural field by applying an external input, simulating 10 Hz input from the premotor cortex that has been observed in patients with PD (Timmermann et al., 2003 ). In order to test if our closed-loop stimulation method affects inputs from other brain areas, we also added two control inputs in the gamma range (50 and 85 Hz respectively), possibly simulating auditory inputs (Ermolaeva and Borgest, 1980 ; Crone et al., 2001 ).

To compute the closed-loop stimulation term, we evaluated the undesirable components from the neural field equation:

where is the band-pass filtered component (not containing the “pathological” signal contribution) of . Thus, using multi-site measurements of membrane potential (provided by the MEA), and by fitting the neural field equation parameters to cortical activity measurements (estimating the parameters γ, c and the parameters of functions W(r − r′) and S[V]), it is possible to use Eq. (6, left term) to evaluate the total input (sum of synaptic and external inputs) to neurons at each point of space and time. In the following, we assume that these parameters (connectivity function, synaptic coupling, synaptic time constant) are known. We extracted the last 100 ms window of this estimated input signal at a sampling rate of 1,000 Hz. High- and low-pass filters (with cut-off frequencies at 9 and 11 Hz, respectively) were then applied to extract the numerical values of the signal’s “pathological” component. Figure 4 illustrates the activity of the 1 cm² neural field during 1 s of stimulation, where 10 Hz incoming activity was implemented to simulate the “pathological” state (500 ms where closed-loop stimulation is “off” then 500 ms “on”).

FIGURE 4

Figure 4. Network activity (2D mapped to 1D) before (“control off”) and during (“control on”) closed-loop stimulation. “Pathological” 10 Hz oscillations of the recorded membrane potential (0–500 ms) are durably attenuated by the closed-loop stimulation signal (switched “on” between 500 and 1000 ms) controlling network activity in space and time. When the stimulation signal is present, the network is still active in other (“physiological”) frequency bands, i.e., network activity are not completely suppressed by the stimulation signal.

The closed-loop stimulation signal efficiently attenuates the 10 Hz oscillation that was the target frequency. Inputs at 50 and 85 Hz (with amplitude of 0.1 mV for both inputs) were added to the cortical network (to evaluate potential perturbation by the closed-loop stimulus to other rhythms present in the signal), and the spectrum of neural activity was computed for two periods: without and with stimulation (as illustrated in Figure 5 ).

FIGURE 5

Figure 5. Power spectrum of neural activity before (blue) and during (red) stimulation. Before closed-loop stimulation, the dominant frequency is 10 Hz. During closed-loop stimulation specifically targeting the 10 Hz frequency, the power spectrum is strongly decreased at 10 Hz (and also at its sub-harmonics) whereas other frequency bands appear minimally affected. Note: (1) = strong attenuation of the 10 Hz peak, (2) and (3) = slight perturbation of the 50 and 85 Hz peaks respectively. This means that the closed-loop stimulation signal minimally interferes with the processing of inputs to the network at non-“pathological” frequencies.

As expected, the power spectrum of neural activity in Figure 5 shows a strong peak at 10 Hz in the absence of closed-loop stimulation. Conversely, during closed-loop stimulation, this 10 Hz peak (and also sub-harmonics) is strongly attenuated (roughly by a factor of 10). However, the amplitude of other frequency components is much less affected by the stimulation. Indeed, the response of the cortical target to the 50 Hz input is only slightly affected by the stimulation (7.9% increase in peak area). The 85-Hz input response appears more affected (14.5% increase in peak area). One possible explanation is that the linear approximation of the sigmoid function in the neural field equation is reasonable, but that the introduction of supplementary higher-order terms to approximate the sigmoid function may decrease this perturbation of other rhythms. This likely explains why the 10 Hz frequency component is strongly attenuated instead of completely suppressed. The power spectrum presented in Figure 5 shows that the cortical target’s response to 50 Hz inputs is comparable to the closed-loop stimulation signal. Furthermore, the cortical target’s response to 85 Hz inputs is moderately affected (however, approximately six times less than the target frequency to attenuate, that is decreased by about 90%, see Figure 5 ). This illustrates the selectivity of our cortical stimulation method, which minimally interferes with “physiological” brain rhythms, assuming that “physiological” rhythms do not include 10 Hz components, which appears to be the case for neuronal activity in M1. An important clinical consequence of minimizing the interference with other brain rhythms is the potential reduction of undesirable side effects associated with such perturbations (see the Introduction section).

With this method, we map the cortical target with a two-dimensional MEA and apply feedback in space and time to the cortical target. The feedback term is computed according to Eq. 4. Another possibility is to consider the MEA as a dynamical system behaving according to the same equations as the cortical target itself, i.e., Eq. 1. In the following, we explore this second option, with the aim to develop an innovative neuroprosthetic device, which we call a virtual cortical module (by “virtual”, we mean that both the module and the cortical area targeted are governed by the same equations, resulting in a neuroprosthetic device mimicking the activity of an artificial tissue interfaced with a biological cortical tissue. This artificial tissue, that assists the biological tissue to function correctly, comes into play only when the presence of pathological rhythms is detected by this same artificial device).

Discussion

Virtual Cortical Module, a Perspective for a New Neuroprosthetic Device

In the first section of this paper, we presented a method for selective correction of brain rhythms by a closed-loop system. Now let us consider the most general case of closed-loop stimulation, where an arbitrary stimulation signal is sent at position r at time t. This feedback term F(r,t) may be written in the form:

where Φ(r,t) is the effective potential of an external neural field interacting with the field ψ(r,t) such that Φ(r,t) depends on ψ(r,t), since the feedback may depend on potential values throughout the network, and thus on delayed potential values due to the finite conduction time along axons. We assume that the field Φ(r,t) also obeys to Eq. 1, but where the input function, instead of I(r,t), is a function of ψ(r,t), the actual cortical field activity (this function may be, for example, a convolution over space with a connectivity function, similarly to that in a classical neural field equation) as an input instead of I(r,t). The function , combining potential values into an appropriate stimulation value, is unknown. One possible form for this function is:

Equation 8 means that the feedback received by the cortical target, at each location, has the form of an interaction with another cortical network Φ(r,t), not necessarily having the same activity. In our theoretical framework (see Introduction), we suppose that a MEA is appropriate to map the activity of the cortical target, and we propose a stimulation term especially designed to attenuate predetermined frequency bands. It is possible to generalize this, and to consider the following: the MEA may be seen as an “artificial” neural field Φ(r,t) in two dimensions, with each microelectrode emulating the activity of a neural mass. The connectivity function, time constants, firing threshold, and external inputs of this artificial neural field are completely configurable. The potential of this artificial neural field may be, via the sigmoid function, converted to a firing rate, i.e., a stimulation frequency. Indeed, the dynamics of both the cortical target and the stimulation device would be governed by the same equations and the same physical principles. In practice, this new class of devices may be implanted (eventually through the dura) on the cortical target using a MEA (as done in Kellis et al., 2009 ), the MEA having the possibility to be flexible to adapt to target topology (see Figure 6 ), in order to control local cortical rhythms in a predetermined manner (attenuation or augmentation of rhythms). The application we are interested in is to control the 10-Hz frequency component in the motor cortex of patients with PD. Such a virtual cortical module could potentially be used to normalize, augment, or restore brain functions, since it is optimally interfaced with brain tissue.

FIGURE 6

Figure 6. Block diagram of the proposed virtual cortical module. The MEA is positioned on the cortex, and communicates via an interface INT with the chipset CPU. Recorded signals are amplified by the AMPL module, and then analyzed to detect the presence of “pathological” rhythms. If the undesirable rhythm is present, the TRIG module triggers the solving of neural field equations with the NFE module. The computed potential values are stocked in the BUFF module, to be later recalled when delayed values of the potential are needed. The potential values are then converted by the module CONV into a firing rate, i.e., a stimulation frequency. This stimulation frequency is transmitted to the STIM module that will trigger the stimulation to the cortical target. An internal clock SYNC synchronizes the various modules (Fig. from Beuter and Modolo, patent EP 09305432.8, 2009).

The parameters of the virtual cortical module should be adjusted to produce the desired effect, i.e., to attenuate, enhance, or reshape neural activity in the cortical target. Such adjustment may be achieved by either analytical or numerical study of two-population networks, where one population represents the cortical target, and the other the virtual cortical module. Interestingly, the technology required to elaborate this new class of devices is already available. Indeed, several groups have developed neuromorphic circuits, i.e., electronic devices emulating the behavior of spiking neurons, or of small spiking neuron networks (see Zou et al., 2006 for an example of circuits emulating Hodgkin–Huxley neurons). Dynamical equations are “naturally” solved by analog circuits performing mathematical operations (addition, differentiation…). These analog chipsets allow real-time simulation of spiking neurons, but their large size (several square millimeters for a chipset simulating a single neuron, Saighi et al., 2005 ) makes their use impossible today for therapeutic or augmentation purposes. Conversely, virtual cortical modules, while implementing only coarse-scale neural activity, represent a realistic and practical way to implement in brain tissue artificial neural networks with plausible dynamical behaviors. To this end, analog circuits solving the delayed integro-differential neural field equation may accurately emulate cortical circuits, to interact with cortical tissue with predetermined objectives.

Future Directions

Until now, DBS has been surprisingly effective despite the fact that, in large part, its underlying physiological processes are still unknown. However, this technique also has its drawbacks and side effects, and cannot be proposed to a significant proportion of PD patients. Furthermore, the parameters are adjusted periodically and empirically by the clinician based on clinical outcomes, which may be seen as a way to “close the loop” once in a while. In this context, the natural evolution of brain stimulation devices is to propose a tailored therapy for individual patients with PD. The method proposed here, based on neural field theory from the theoretical point of view, and on the use of a dual function high-density MEA from the technical point of view, illustrates that such a tailored therapy might be achieved with high-density MEA sequentially recording and stimulating appropriately the motor cortex of PD patients in order to control the rhythm associated with “pathological” tremor. Furthermore, in principle, this method should not be limited to PD since the principle of using a MEA to target a brain area in order to normalize its activity is general enough to be applied clinically to other neurological disorders. Depending on the disease, and thus on the cortical area targeted, it is conceivable to adjust the size of the MEA to the size of the target, for example. Also, with the proposed closed-loop stimulation method, it is possible not only to decrease neural activity in predetermined frequency bands, but also to increase it. This latter possibility may be desirable in neurological disorders decreasing neural activity in well-identified frequency bands with functional consequences. Another consideration regarding the clinical application of our method is the possible damage induced by the insertion of the MEA. This issue has been partly addressed by House et al. (2006) that studied cortical tissue after having implemented a MEA. Even if the insertion of the MEA left visible traces in cortical tissue, the authors did not find any significant hemorrhage or any clinical complication. This work confirms that the device implementing our proposed method should not induce damage or clinical complications.

To unambiguously distinguish between normal and “pathological” activity might be in general a difficult task. However, PD is unusual in that the “pathological” activity does tend to be concentrated in a well-defined frequency band, though there is, of course, some variation in this frequency both within and between subjects. Furthermore, for the particular case of Parkinson’s disease, additional confirmation of the appropriate frequency might be obtained by adding a monitoring signal in order to detect the onset, presence, and frequency of tremor. One may imagine that a small accelerometer worn in a finger ring might communicate with the neuroprosthetic device, in order to confirm that control stimulation is needed. By doing so, one may avoid “false positives”, i.e., the device delivering stimulation when it is not needed. This could also prevent, or at least minimize, the interference of the stimulation signal with “physiological” frequency bands.

Furthermore, one challenge is that there might be differences in the pathological frequency and its power from patient to patient. Consequently, it is not possible to propose a unique setting (target pathological frequency and power threshold) for the triggering of the closed-loop stimulation, since it may interfere with physiological brain oscillations. For instance, the mu rhythm, present in the primary motor cortex at rest, is centered at 10 Hz in the frequency band aimed by the closed-loop feedback, and may be affected if the settings are not appropriate. This might have adverse effects on other functional modalities since the mu rhythm is speculated to play an important role in visuo-motor integration (Pineda, 2005 ). Moreover, since desynchronization of the mu rhythm is associated with motor activity, it remains to be seen experimentally how modulation of the parkinsonian “pathological” frequency component affects the motor circuitry of patients. It might be added that the exact role and origin of the mu rhythm are still uncertain (Hari et al., 1997 ; Hari, 2006 ; Caetano et al., 2007 ). Thus, a potential area for future work would be to study the effects of brain stimulation in the context of parkinsonian tremor on the mu rhythm and its consequences on motor and/or sensory activities of patients. Another example is the complex interplay between low-frequency and high-frequency oscillations in the motor cortex during movement (Miller et al., 2007 ). This suggests that, depending on the motor activity, oscillations in the pathological frequency band should not always be suppressed. Thus, it may be required to measure the baseline for the various frequency components of brain activity in various states of physical movement, in order to adapt the method as effectively as possible, avoiding the suppression of needed rhythms.

Possible extensions of the model presented here may include a contribution of axonal activation by electrical stimulation. Indeed, DBS seems to activate axons at the same frequency as the DBS current (McIntyre et al., 2004 ). One possible way to achieve this would be to consider the cortical tissue as a superposition of two fields: one field of somas and one of fibers. Depending on stimulus duration (the main factor for the activation of neural elements, indeed, somas and axons have very different chronaxies, see Ranck, 1975 ), one may be able to add the electrical stimulus to the soma field or the fiber field, or even a combination of both, using an appropriate transfer function. Another possible extension is to convert the calculated signal to a charge-balanced signal, to be used in clinical practice. Indeed, the stimulation signal derived from the neural field model is continuous, and stimulation signals used in DBS systems are discrete, under the form of a train of charge-balanced pulses. Another possibility would be to use a multi-compartment neuron model that accounts for differential activation of somatic and axonal parts of neurons, such as used by McIntyre et al. (2004) . Finally, even if we are currently focused on Parkinsonian tremor, we need to address the question of the attenuation of other parkinsonian motor symptoms in order to propose a complete and effective method. This point will be addressed in forthcoming papers. Also, the proposed method does not constitute an end in itself, and definitely needs to be tested experimentally to provide a proof of concept.

In a broader context, we propose that devices based on recording/stimulating, flexible, high-density MEA represent an appealing possibility to modulate brain rhythms in space and time by emulating the activity of biological cortical circuits. Virtual cortical modules may provide “replacement parts” for the cortex, or even represent a possibility of augmentation. While the example presented in this paper is related to tremor control in PD, other possible applications may include elimination of phantom limb pain after amputation, elimination of neuropathic pain and other neurological disorders caused by cortico-thalamo-cortical dysrhythmia. Development of these devices will require integration of signal processing, electronic conception, and a combination of biophysical and more abstract mathematical models of cortical tissue.

Appendix

Numerical Implementation of the Neural Field Model

To solve the neural field Eq. 1 numerically, time and space are discretized so that the cortical sheet is represented as a 2D grid of N = N_x × N_y neural masses, where N_x and N_y denote neural masses along the X and Y co-ordinates, respectively. Let W(d) be the spatial connectivity between the nth and the mth neural mass separated by a distance d, expressed as (Atay and Hutt, 2005 ):

where a_e and a_i denote excitatory and inhibitory synaptic weights, respectively; and r = σ_e/σ_i (Atay and Hutt, 2005 ) the ratio of spatial ranges for excitatory and inhibitory connections, respectively. For instance, if r < 1 and r < a_e/a_i, the connectivity is locally excitatory and laterally inhibitory. Parameter values have previously been published (Atay and Hutt, 2005 ), i.e., a_e = 60, a_i = 55, and r = 0.5. The spatial distance between two locations M = (x_m, y_m) and N = (x_n, y_n) of the discretized cortical sheet was computed as the Euclidian distance . Considering the finite conduction speed c of action potentials along axons (here, c = 100 mm/s), the resulting delay in signal propagation across a distance d is τ = d/c. The effect of the activity of the mth neural mass at some time t − τ on the firing rate σ_n(t) of the nth neural mass at time t is denoted by σ_n,m(t − τ) and is defined as:

where υ_max = 100 Hz is the maximum firing rate of a single neuron, V₀ the threshold for spike initiation, V_m(t − τ) the mean voltage of the mth neural mass at time t − τ, and λ the steepness parameter for the sigmoid function defining the firing rate.

Then, let V_n(t) be the mean potential voltage of the nth neural mass at time t. We provided an initial stimulus κ_n(t) to the nth neural mass for each n at time t, to induce the 10-Hz frequency component experimentally measured in M1 of patients with PD and coherent with muscular activity at the periphery. The neural field evolves over time according to, in discretized form:

where χ_n,m(t) is defined as:

and where ζ_n(t) accounts for uncorrelated white spatio-temporal noise of mean value 0.05 mV, added to simulate background synaptic activity.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by the European Network of Excellence BioSim (contract No. SHB-CT-2004-005137) to Anne Beuter and by a Post-Doctoral Research Fellowship from the Government of Canada (program of the Department of Foreign Affairs and International Trade 2009-2010) to Julien Modolo.

References

Amari, S. (1977). Dynamics of pattern formation in lateral-inhibition type neural fields. Biol. Cybern. 27, 77–87.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Atay, F. M., and Hutt, A. (2005). Stability and bifurcations in neural fields with finite propagation speed and general connectivity. SIAM J. Appl. Math. 65, 644–666.

CrossRef Full Text

Birdno, M. J., Kuncel, A. M., Dorval, A. D., Turner, D. A., and Grill, W. M. (2008). Tremor varies as a function of the temporal regularity of deep brain stmulation. Neuroreport 19, 599–602.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bressloff, P. C. (2001). Traveling fronts and wave propagation failure in an inhomogeneous neural network. Physica D 155, 83–100.

CrossRef Full Text

Caetano, G., Jousmaki, V., and Hari, R. (2007). Actor’s and observer’s primary motor cortices stabilize similarly after seen or heard motor actions. Proc. Natl. Acad. Sci. U.S.A. 104, 9058–9062.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Canavero, S., Paolotti, R., Bonicalzi, V., Castellano, G., Greco-Crasto, S., Rizzo, L., Davini, O., Zenga, F., and Ragazzi, P. (2002). Extradural motor cortex stimulation for advanced Parkinson disease. Report of two cases. J. Neurosurg. 97, 1208–1211.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Crone, N. E., Boatman, D., Gordon, B., and Hao, L. (2001). Induced electrocorticographic gamma activity during auditory perception. Clin. Neurophysiol. 112, 565–582.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Elble, R. J., and Koller, W. C. (1990). Tremor. Baltimore: Johns Hopkins University Press.

Ermentrout, G. B., and Cowan, J. D. (1979). A mathematical theory of visual hallucination patterns. Biol. Cybern. 34, 137–150.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ermolaeva, V. Y., and Borgest A. N. (1980). Intercortical connections of the auditory areas with the motor area. Neurosci. Behav. Physiol. 10, 210–215.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Feng, X. J., Greenwald, B., Rabitz, H., Shea-Brown, E., and Kosut, R. (2007). Toward closed-loop optimization of deep brain stimulation for Parkinson’s disease: concepts and lessons from a computational model. J. Neural Eng. 4, 14–21.

CrossRef Full Text

Gredal, O., Pakkenberg, H., Karlsborg, M., and Pakkenberg, B. (2000). Unchanged total number of neurons in motor cortex and neocortex in amyotrophic lateral sclerosis: a stereological study. J. Neurosci Methods 95, 171–176.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Guehl, D., Edwards, R., Cuny, E., Burbaud, P., Rougier, A., Modolo, J., and Beuter, A. (2007). Statistical determination of optimal STN stimulation site in Parkinson’s disease. J. Neurosurg. 106, 101–110.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hari, R. (2006). Action–perception connection and the cortical mu rhythm. Prog. Brain Res. 159, 253–260.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hari, R., Salmelin, R., Makela, J. P., Salenius, S., and Helle, M. (1997). Magnetoencephalographic cortical rhythms. Int. J. Psychophysiol. 26, 51–62.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hatospoulos, N. G. (2009). “Rhythms in motor processing: functional implications for motor behavior,” in Short course organized by the Society for Neuroscience “Rhythms of the Neocortex: Where Do They Come From and What Are They Good For? Short course organised by the Society For Neurosciencequotidns 39th Annual Meeting, Chicago, IL.

Horton, J. C., and Adams, D. L. (2005). The cortical column: a structure without a function. Philos. Trans. R. Soc. Lond. B Biol. Sci. 360, 837–862.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

House, P. A., MacDonald, J. D., Tresco, P. A., and Normann, R. A. (2006). Acute microelectrode array implantation into human neocortex: preliminary technique and histological considerations. Neurosurg. Focus 20, E4.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Huang, X., Troy, W. C., Yang, Q., Ma, H., Laing, C. R., Schiff, S. J., and Wu, J. Y. (2004). Spiral waves in disinhibited mammalian neocortex. J. Neurosci. 24, 9897–9902.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kellis, S. S., House, P. A., Thomson, K. E., Brown, R., and Greger, B. (2009). Human neocortical electrical activity recorded on nonpenetrating microwire arrays: applicability for neuroprostheses. Neurosurg. Focus 27, E9.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

McIntyre, C. C., Savasta, M., Kerkerian-Le Goff, L., and Vitek, J. L. (2004). Uncovering the mechanism(s) of action of deep brain stimulation: activation, inhibition, or both. Clin. Neurophysiol. 115, 1239–1248.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Miller, K. J., Leuthardt, E. C., Schalk, G., Rao, R. P. N., Anderson, N. R., Moran, D. W., Miller, J. W., and Ojemann, J. G. (2007). Spectral changes in cortical surface potentials during motor movement. J. Neurosci. 27, 2424–2432.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Modolo, J., and Beuter, A. (2009a). Linking brain dynamics, neural mechanisms and deep brain stimulation in Parkinson’s disease. Med. Eng. Phys. 31, 615–623.

CrossRef Full Text

Modolo, J., Campagnaud, J., and Beuter, A. (2009b). Effect of complex delayed feedback in a neural field model. Conference Abstract: Bernstein Conference on Computational Neuroscience. Front. Comput. Neurosci. doi: 10.3389/conf.neuro.10.2009.14.113

CrossRef Full Text

Modolo, J., Edwards, R., Campagnaud, J., Bhattacharya, B., and Beuter, A. (2010). Past, present and future of brain stimulation. Math. Model. Nat. Phenom. 5, 185–207.

CrossRef Full Text

Pagni, C. A., Zeme, S., and Zenga, F. (2003). Further experience with extradural motor cortex stimulation for treatment of advanced Parkinson’s disease. Report of 3 new cases. J. Neurosurg. Sci. 47, 189–193.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pereira, E. A. C., Green, A. L., Nandi, N., Aziz, T. Z. (2007). Deep brain stimulation: indications and evidence. Expert Rev. Med. Devices 4, 591–603.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pineda, J. A. (2005). The functional significance of mu rhythms: translating “seeing” and “hearing” into “doing”. Brain Res. Rev. 50, 57–68.

CrossRef Full Text

Popovych, O. V., Hauptmann, C., and Tass, P. A. (2008). Impact of nonlinear delayed feedback on synchronized oscillators. J. Biol. Phys. 34, 267–279.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Priori, A., and Lefaucheur, J. P. (2007). Chronic epidural motor cortical stimulation for movement disorders. Lancet Neurol. 6, 279–286.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ranck, J. B. Jr. (1975). Which elements are excited in electrical stimulation of mammalian central nervous system: a review. Brain Res. 98, 417–440.

Pubmed Abstract | Pubmed Full Text

Rodriguez-Oroz, M. C., Rodriguez, M., Guridi, J., Mewes, K., Chockmann, V., Vitek, J., DeLong, M. R., and Obeso, J. A. (2001). The subthalamic nucleus in Parkinson’s disease: somatotopic organization and physiological characteristics. Brain 124, 1777–1790.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rosenblum, M. G., and Pikovsky, A. S. (2004). Controlling synchronization in an ensemble of globally coupled oscillators. Phys. Rev. Lett. 92, 114102.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rubin, J. E., and Terman, D. (2004). High-frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model. J. Comput. Neurosci. 16, 211–235.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rubino, D., Robbins, K. A., and Hatsopoulos, N. G. (2006). Propagating waves mediate information transfer in the motor cortex. Nat. Neurosci. 9, 1549–1557.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Saighi, S., Tomas, J., Bornat, Y., and Renaud, S. (2005). “A conductance-based silicon neuron with dynamically tunable model parameters,” Proceedings of the 2nd International IEEE EMBS Conference on Neural Engineering, Arlington, Virginia.

CrossRef Full Text

Santaniello, S., Fiengo, G., and Glielmo, L. (2008). “Adaptive feedback control in deep brain stimulation: a simulation study,” In Proceedings of the 17th congress of the International Federation of Automatic Control, Seoul, Korea, July 6–11, 11624–11629.

Schnitzler, A., Timmermann, L., and Gross, J. (2006). Physiological and pathological oscillatory networks in the human motor system. J. Physiol. Paris 99, 3–7.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tass, P. A. (1999). Phase Resetting in Medicine and Biology. Stochastic Modelling and Data Analysis, Series: Springer Series in Synergetics. New York: Springer.

Tass, P. A., and Majtanik, M. (2006). Long-term anti-kindling effects of desynchronizing brain stimulation: a theoretical study. Biol. Cybern. 94, 58–66.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Timmermann, L., and Fink, G. R. (2009). Modulating pathological oscillatory activity in Parkinson’s disease: what is the rhythm? Exp. Neurol. 215, 209–211.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Timmermann, L., Florin, E., and Reck, C. (2007). Pathological cerebral oscillatory activity in Parkinson’s disease: a critical review on methods, data and hypotheses. Expert Rev. Med. Devices 4, 651–661.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Timmermann, L., Gross, J., Dirks, M., Volkmann, J., Freund, H. J., and Schnitzler, A. (2003). The cerebral oscillatory network of parkinsonian resting tremor. Brain 126, 199–212.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Zou, Q., Bornat, Y., Saighi, S., Thomas, J., Renaud, S., and Destexhe, A. (2006). Analog-digital simulations of full conductance-based networks of spiking neurons. Netw. Comput. Neural Syst. 17, 211–233.

CrossRef Full Text