Closed-loop stimulation of a delayed neural fields model of parkinsonian STN-GPe network: a theoretical and computational study

Several disorders are related to pathological brain oscillations. In the case of Parkinson's disease, sustained low-frequency oscillations (especially in the β-band, 13–30 Hz) correlate with motor symptoms. It is still under debate whether these oscillations are the cause of parkinsonian motor symptoms. The development of techniques enabling selective disruption of these β-oscillations could contribute to the understanding of the underlying mechanisms, and could be exploited for treatments. A particularly appealing technique is Deep Brain Stimulation (DBS). With clinical electrical DBS, electrical currents are delivered at high frequency to a region made of potentially heterogeneous neurons (the subthalamic nucleus (STN) in the case of Parkinson's disease). Even more appealing is DBS with optogenetics, which is until now a preclinical method using both gene transfer and deep brain light delivery and enabling neuromodulation at the scale of one given neural network. In this work, we rely on delayed neural fields models of STN and the external Globus Pallidus (GPe) to develop, theoretically validate and test in silico a closed-loop stimulation strategy to disrupt these sustained oscillations with optogenetics. First, we rely on tools from control theory to provide theoretical conditions under which sustained oscillations can be attenuated by a closed-loop stimulation proportional to the measured activity of STN. Second, based on this theoretical framework, we show numerically that the proposed closed-loop stimulation efficiently attenuates sustained oscillations, even in the case when the photosensitization effectively affects only 50% of STN neurons. We also show through simulations that oscillations disruption can be achieved when the same light source is used for the whole STN population. We finally test the robustness of the proposed strategy to possible acquisition and processing delays, as well as parameters uncertainty.


Introduction
Pathological conditions in Parkinson's disease, Huntington's disease or epileptic seizures are often associated with abnormal oscillations in specific brain areas (Bragin et al., 2010;Devergnas and Wichmann, 2014). In particular, parkinsonian motor symptoms result from a dopamine depletion in the striatum and strong neural synchronization, bursting and sustained oscillations in low frequencies can be observed (Bergman et al., 1994;Brown et al., 2001;Jenkinson and Brown, 2011). This oscillatory activity is especially detectable in global signals such as local field potentials (LFPs) (Stein and Bar-Gad, 2013). The precise relation between these sustained oscillations and motor symptoms is still a matter of debate. Nonetheless, there is clear evidence that low-frequencies activity in the LFP of the subthalamic nucleus (STN) is correlated to parkinsonian motor deficit. More precisely, bradykinesia and rigidity are related to strong and sustained activity in β-frequencies (13-30 Hz) (Kühn et al., 2006;Hammond et al., 2007;Little et al., 2012), while low γ -frequencies (31-45 Hz) correlate with tremor severity (Beudel et al., 2015). Sustained β-activity is not reported in healthy conditions (Bar-Gad et al., 2003;Heimer et al., 2006;Schwab et al., 2013). Among the possible hypotheses for the mechanisms underlying this pathological synchronization onset, a possible pacemaker role of the network formed by the excitatory STN and the inhibitory external globus pallidus (GPe) has been proposed based on in vitro experiments (Plenz and Kital, 1999). With striatal origins and cortical STN patterning, this role of the STN-GPe as a pacemaker is one of the main hypotheses for this pathological oscillations onset (Bevan et al., 2002;Lang and Zadikoff, 2005). Experiments conducted on non-human primate models of Parkinson's disease confirm the importance of reciprocal STN-GPe connections on oscillatory activity onset (Nambu and Tachibana, 2014). Although not formally proven in vivo, this pacemaker effect of the STN-GPe network is supported by mathematical models, fitted to neurophysiological data, which effectively reproduce sustained oscillations in parkinsonian conditions Holgado et al., 2010;Merrison et al., 2013;Pasillas-Lépine, 2013). In order to test this hypothesis, there is a critical need for developing theoretical strategies of neuromodulation that could be implemented experimentally in order to control sustained oscillations in the STN.
Deep Brain Stimulation (DBS, Benabid et al., 1991) is a neuromodulation method that consists in chronically implanting electrodes in specific targeted brain areas and in applying electrical currents to alleviate neurological symptoms (Lozano and Lipsman, 2013). In the case of Parkinson's disease, three different targets (thalamus, globus pallidus pars internalis and subthalamic nucleus) have been identified to reduce motor symptoms such as tremor or/and rigidity, akinesia, the STN being the most prescribed. Although the electrical DBS's mechanisms of action are still a matter of debate, it has been suggested that DBS of the STN may affect β-synchrony in the targeted neural population (Eusebio et al., 2011). In most existing DBS treatments, the electrical signal is delivered in an open-loop manner, meaning regardless of the patient's state or brain activity.
As reviewed in Carron et al. (2013), several attempts have recently been made to optimize DBS by taking into account some measurements of the patient brain activity. These include adaptive and on-demand stimulation (Marceglia et al., 2007;Graupe et al., 2010;Rosin et al., 2011;Santaniello et al., 2011;Little et al., 2013), delayed and multi-site stimulation (Omel chenko et al., 2008;Batista et al., 2010;Pfister and Tass, 2010;Lysyansky et al., 2011;Tass et al., 2012), optimal control strategies (Feng et al., 2007), and closed-loop firing rate regulation (Wagenaar et al., 2005;Luo et al., 2009;Franci et al., 2011;Pasillas-Lépine et al., 2013). Another strategy, not included in that review, is Schiff (2010), which advocates for the use of Kalman filtering to enhance DBS strategies. The review Carron et al. (2013) underlines a clear gap between model-based closedloop DBS strategies, whose efficiency is most of the time assessed in silico, and the DBS signals actually tested in vivo.
Optogenetics is a promising tool to unravel pathological oscillations onset mechanisms and to explore closed-loop strategies to optimize DBS treatment. Using a viral vector, genes coding for neuronal membrane ionic channels, sensitive to precise wavelengths, are selectively transferred to precisely targeted neurons, which become photosensitive. Thus, simple pulses of intense light, e.g., through implanted optical fibers, can induce the response of the photosensitized neurons (Boyden et al., 2005;Pastrana, 2010). By contrast with electrostimulation, optical neuromodulation does not induce any stimulation artifacts on electrical recordings. It also allows sharpest targeting of specific neuronal subpopulations. In addition, the development of bi-color optogenetics stimulation and bi-switchable opsins such as Halo-ChR2 (Han and Boyden, 2007), now allows for simultaneous excitation and inhibition of the neural population in a controlled way, whereas electrical DBS may perform both excitation and inhibition on different compartments of the same neurons or on different neurons depending on their orientations relative to the electrical field (Ranck, 1975;McIntyre et al., 2004). Thus, combined with electrode measurements, optogenetics may offer unprecedented possibilities to achieve closed-loop stimulation of neural structures, in particular for counteracting pathological oscillations in animal models (Chaillet et al., 2014). First experimental attempts include those described by Paz et al. (2012), where closed-loop optical stimulation was used on rats in order to interrupt seizures after cortical injury.
A particular way to avoid the generation of sustained oscillations in a parkinsonian STN-GPe network is to artificially reduce the gain of the STN-GPe loop. This strategy has been developed in Pasillas-Lépine et al. (2013) 1 using tools from linear control theory. That theoretical paper demonstrates that, by delivering a photostimulation whose frequency is proportional to the measured unitary activity, sustained oscillations are impeded in the parkinsonian basal ganglia model proposed in Holgado et al. (2010). It further proposes a filtering strategy to cope with inherent acquisition and processing delays. In view of its simplicity and ease of technological implantation, this closed-loop DBS strategy is very appealing. Nonetheless, this encouraging result still suffers from important limitations. The main one stands in the employed model, which summarizes the activity of each neural population by a single averaged variable and thus fails at taking into account spatial heterogeneity within STN and GPe. The first aim of this theoretical paper is to propose a model of the parkinsonian STN-GPe network, able to take into account both temporal and spatial evolutions of the activity.
Another limitation of the study conducted in Pasillas-Lépine et al. (2013) is that it relies on a linearized model, and thus neglects the inherent nonlinear effects arising from the neurons f-I curves and the saturation of the stimulation signal. The works Faye and Faugeras (2010); Veltz and Faugeras (2011) provide a relevant mathematical framework to study the proposed spatiotemporal model analytically, despite the combined effects of delays and nonlinearities. The second objective of this paper is to propose a closed-loop stimulation strategy enabling disruption of sustained oscillations in this spatiotemporal model of the STN-GPe network, and to prove its efficiency both mathematically and numerically. This stimulation signal is designed to be realistically implementable in order to test the role of sustained oscillations in parkinsonian symptoms in future experiments.
The third objective of this paper is to validate both the model and the associated closed-loop stimulation signal based on in silico experiments that take into account specific experimental constraints that would occur in an in vivo testing. These constraints include parameter uncertainties, imperfect photosensitization or insufficient illumination of STN neurons, use of a single light source for the whole STN population, and unavoidable acquisition and processing delays.

Neural Fields Model of the STN-GPe Network
In parkinsonian conditions, pathological oscillations are commonly observed in the dorsolateral region of the STN (Zaidel et al., 2010). We focus on this specific area and on the GPe neurons that either project or receive projections from that region. From now on, we will simply refer to these regions as STN and GPe for simplicity. A finer study of the pathological oscillations onset in the parkinsonian STN-GPe network requires to include the evolution of their activity in both time and space. To this aim, we model the STN and GPe populations as delayed neural fields (Coombes, 2005;Bressloff, 2012). Indicating the STN and GPe populations with indexes 1 and 2 respectively, the STN-GPe network is thus modeled as coupled delayed neural fields: τ 2 ∂z 2 ∂t (r, t) = −z 2 (r, t) z 1 (r, t) and z 2 (r, t) represent the neuronal activity of the neurons at position r and at time t for STN and GPe respectively. denotes a compact set of R p , p ∈ {1, 2, 3}, representing the physical support of the network neuronal populations. τ 1 and τ 2 are the time constant STN and GPe activities. w ij (r, r ′ ) describes the synaptic strength between location r ′ in population j and location r in population i. d 1 (r, r ′ ) and d 2 (r, r ′ ) represent the axonal time delays between a pre-synaptic neuron at position r ′ and a post-synaptic neuron at position r. S 1 and S 2 are the activation functions of STN and GPe respectively. I ext 1 (r, t) and I ext 2 (r, t) account for the external inputs of STN and GPe. In this work, we consider that the STN neurons expresses the gene that codes for photoactivatable proteins and receive sufficient illumination from the stimulation device. In other words, the STN activity can be modulated with light stimulation. I ext 1 (r, t) thus decomposes into a stimulation signal u(r, t) that we can regulate and the inputs I 1 (r, t) from other cerebral structures, especially cortical inputs. The GPe is not directly affected by the stimulation and hence sums up to inputs I 2 (r, t) from other brain areas, specifically from the striatum. The overall setup is depicted by Figure 1.

Closed-loop Stimulation Signal
As will be seen in Section 3, the STN-GPe model (Equation 1) may exhibit sustained oscillations even when the exogenous inputs I ext 1 and I ext 2 are constant in time. These oscillations result from the combined effects of transmission delays and too strong synaptic weights between STN and GPe. As originally FIGURE 1 | STN-GPe network. This figure illustrates the physical arrangement of STN and GPe. We assume that the zones involved in pathological oscillations in STN and GPe measure 2.5 mm each (Zaidel et al., 2010); they are indicated in black clouds. The light-green box illustrates a controller that reads data from STN and applies a photostimulation to the same neural population. Elements of this figure have been found in Wikipedia and they are licensed according to CC (http://creativecommons.org/licenses/ by-sa/3.0). proposed in Holgado et al. (2010), these increased synaptic gains in pathological conditions may result from dopamine depletion. In order to attenuate these pathological oscillations, we rely on a proportional stimulation law, in the spirit of Pasillas-Lépine et al. (2013). Proportional control constitutes the basics of control theory; it consists in applying an input which is proportional to the difference between the measured state and a targeted reference: where k c is a positive control gain that can be adjusted at will, and z ref (r) represents a prescribed rate for the STN population. The stimulation signal is thus assumed to act as an external rate on the STN. It may act as an either inhibitory or excitatory input, depending on the sign of u(r, t). This is in agreement with the recent development of bi-color optogenetics stimulation and bi-switchable opsins such as Halo-ChR2 (Han and Boyden, 2007), which enables both excitation and inhibition of the same neuronal population. The rationale behind this proportional closed-loop stimulation is rather intuitive: when the measured STN activity z 1 (r, t) overpasses the prescribed rate z ref (r), the stimulation tends to decrease the STN activity, translating in a negative sign of the stimulation signal u(r, t). On the contrary, when the measured STN activity is too low compared to the target rate, the stimulation excites the STN (u(r, t) comes with a positive sign). In both cases, the amplitude of the delivered stimulation signal is proportional to the difference between the STN measured activity and the prescribed one. Note that the implementation of this closed-loop stimulation signal requires the real-time measurement of the STN activity. No measurements on the GPe or other brain areas are needed.
Plugging the closed-loop stimulation signal (Equation 2) into the spatiotemporal model (Equation 1) leads to the following dynamics: α(r) ≥ 0 is a position-dependent gain accounting for the strength of the stimulation influence at position r. It typically decreases as a function of the distance to the stimulation device, due to light absorption in brain tissues (Deng et al., 2014). Its value may also be affected by the quality of photosensitization of the targeted neurons. In the case when photosensitization is not effective or light stimulation is not received, it holds that α(r) = 0.

Tools to Assess Performance and Robustness
In order to numerically investigate the effects of the proposed closed-loop DBS method described above, we propose five in silico experimental protocols.
• Protocol A tests the ability of the proposed spatiotemporal model to generate sustained oscillations in the β band. This test is conducted both on a nominal set of parameters deduced from data available in the literature (see Section 2.4) and on variations around these nominal values. • Protocol B assesses the efficiency of the proposed closed-loop stimulation signal to disrupt these sustained oscillations. • Protocol C studies the performance of the closed-loop stimulation signals in case of poor STN photosensitization. • Protocol D tests the implementability of the closed-loop stimulation policy when a single light source is available for the whole STN population. • Protocol E investigates how acquisition and processing delays may affect the performance of the closed-loop stimulation signal.

Parameters Selection
For simplicity and clarity of exposition, we consider only 1dimensional spatial distributions: is picked as the interval [0, 15]mm, as the STN and GPe are both contained in zone of length 15mm in adult human brain. The STN activity lies in the subregion 1 = [0, 2.5]mm, while the GPe stands in the subregion 2 = [12.5, 15] mm. All these values are in accordance with adult human brain physiology (Mai et al., 1997).
An advanced identification of parkinsonian activity of STN and GPe rate dynamics has been performed in Holgado et al. (2010) based on experimental data available in the literature. We partly rely on that work to identify some parameters of the model. More precisely, in that reference, the time constants τ 1 and τ 2 were taken as the membrane time constants of the neurons involved and can be taken as 6 ms and 14 ms according to e.g., Paz et al. (2005) and Kita and Kitai (1991) respectively. The activation functions S 1 and S 2 were fitted to experimental data available in Kita and Kitai (1991) and Deister et al. (2009) among others. They were taken as sigmoidal functions: where m 1 = 300 spk/s, b 1 = 17 spk/s, m 2 = 400 spk/s and b 2 = 75 spk/s. It can easily be checked that the maximum slopes of these functions are ℓ 1 = ℓ 2 = 1. The remaining parameters cannot be readily taken from Holgado et al. (2010) as they are space-dependent, while Holgado et al. (2010) relied on an averaged model. The axonal delays d i (r, r ′ ) are derived based on the physical distance between r and r ′ and the velocity c i of spike propagation along the axons under concern: We picked c 1 = 2.5 m/s and c 2 = 1.4 m/s, in accordance with Kita et al. (1983). The inputs I 1 (r, t) and I 2 (r, t) encompass activities from other brain areas, mainly from cortex for STN and from striatum for GPe. Although they were assumed constant in time in our theoretical analysis (see Section 3.1), we have taken them as white noises with variance 0.05 centered at 27 spk/s for cortex and 2 spk/s for striatum (the values have been taken from Holgado et al., 2010). Also in line with Holgado et al. (2010), these input rates were multiplied by uniform synaptic weights of 12.5 (cortical input to STN) and 110 (striatal input to GPe). The addition of a noise to these constant inputs makes the cortical and striatal influences more physiologically plausible and in turn allows to test the robustness of the proposed closed-loop stimulation strategy. The striatal influence on GPe being mostly inhibitory (Kita, 2007), I 2 (r, t) was taken with a minus sign. The kernel w 11 accounting for synaptic coupling within STN was taken as zero since there is no evidence of self-excitatory STN connections (Marani et al., 2008). The other kernels, modeling synaptic weights between the populations (w 12 and w 21 ) or within the GPe (w 22 ), were taken as Gaussian functions as represented by Figure 2 (left panel). Synaptic strengths w 12 and w 21 follow the principle of short-range excitation/inhibition. This implies that if a neuron of STN project to a neuron of GPe, the former excites the neighbor neurons in GPe as well. More precisely, w 12 and w 21 were taken as 2 where µ 1 = 1.25mm and µ 2 = 13.25mm indicate the centers of STN and GPe respectively, and the Gaussian functions g ij are defined as The variance σ 12 and σ 21 were taken identical (0.03 mm 2 ) based on the numerical and computational studies Park et al., 2011), which indicate that STN and GPe neurons are connected to a small portion of GPe and STN neurons respectively. The variance σ 22 within GPe is taken smaller (σ 22 = 0.015) as GPe lateral projections are more local (François et al., 1984). The amplitudes K ij of the synaptic weights distribution were picked in order for the mean activity to fit the simulations presented in Holgado et al. (2010): K 12 = 30, K 21 = 38, and K 22 = 2.55. Note that the amplitude of w 22 is much lower than the other two, as the number of lateral connections within GPe 2 All these functions are taken as zero outside their domain.
is small (Sadek et al., 2005). w 12 and w 22 were taken negative as the projections from GPe neurons are mostly GABAergic (Sato et al., 2000a); on the other hand w 21 was picked positive, as STN projections are mostly glutamatergic (Sato et al., 2000b). Note that, although this has little influence on the simulations result, w 22 is not a pure Gaussian function: it is zero when r = r ′ . This was done in order to provide little synaptic coupling with GPe neurons in the immediate vicinity as little correlations have been observed between neighboring neurons in the GPe ( Bar-Gad et al., 2003). This parameters choice lead to a L 2 -norm of lateral GPe connections ℓ 2 w 22 = 1 × 0.95 < 1, thus fulfilling the requirement (Equation 4) of our theoretical result. This rather small value can be interpreted as the absence of endogenous oscillations within GPe (see Section 3.1 for details). In line with e.g., Plenz and Kital (1999); Nambu and Tachibana (2014), the observed pathological oscillations do not find their roots within GPe, but result from interactions between neuronal populations. Finally, the function α(r) accounting for the impact of light stimulation was taken as a Gaussian function (zero outside 1 ) of amplitude 1, centered at the middle of the STN population (1.25mm), and variance 1.25 mm 2 : see right panel of Figure 2. This choice corresponds to a photo-sensitization of all STN neurons and takes into account light absorption in the medium (as the central STN neurons are more impacted by stimulation than neurons in the periphery).
The parameters used in each in silico protocols are summarized in Table 1.

Simulation Details
In order to simulate the delayed neural fields described by Equation (3), we normalized the spatial domain to sim = [0, 1]. All synaptic strengths and delays are computed according to this normalized interval. We discretized sim using m spatial nodes representing unitary subregions of the populations: m/6 for the STN, m/6 for the GPe and the rest for silent neurons modeling the physical distance between STN and GPe. The temporal integration of the system has been done by the Forward Euler method. All the simulation parameters are gathered in Table 2.
We run all the simulations on a Dell Desktop Machine equipped with an Intel i5 processor and 8 GB memory. The average running time of a simulation was 0.7 s. All the simulations are written in Python (Numpy, Scipy, Matplotlib). The source code is freely distributed under the GPL License and can be found on-line at Github 3 .

Provable Disruption of Pathological Oscillations
Our main theoretical result states that the delayed neural fields (Equation 3) can always be stabilized by the proportional closedloop stimulation (Equation 2), provided that the GPe does not exhibit endogenous sustained oscillations. More precisely, we have the following result, whose proof is provided in Appendix.  Mai et al., 1997). Right panel shows a realization of function α(r), the black curve indicates a STN photosensitization of 100%, the red line corresponds to a random photosensitization of 50%.  Consider any admissible target pattern z ref (r). Let the activation functions S 1 and S 2 be bounded, increasing, and with maximum slope ℓ 1 and ℓ 2 respectively. Assume that the lateral synaptic weights in the GPe satisfy and that the photosensitization is efficient over all STN, meaning that α(r) > 0 for all r ∈ . Assume further that the inputs I 1 (r, t) and I 2 (r, t) from other brain areas are constant in time. Then, for any feedback gain k c sufficiently large, the neural fields under closed-loop stimulation (Equation 2) is stable and, from almost all admissible initial conditions 4 , the STN and GPe activities exponentially converge to constant patterns.

Condition (Equation 4
) requires a sufficiently small L 2 -norm for the lateral synaptic weights w 22 of the GPe. It has a physical interpretation: it imposes that the internal dynamics of the GPe should not exhibit instability (see e.g., Faye and Faugeras, 2010). Thus, proportional closed-loop stimulation can compensate for oscillations generated either within the STN (large w 11 ) or due to 4 We call admissible any initial condition that lies within the ranges of S 1 and S 2 .  too strong interactions between the STN and GPe (large w 12 or w 21 ), but might be insufficient to tackle endogenous oscillations within the GPe (large w 22 ).

Protocol A: Sensitivity Analysis
We start by assessing the capability of the proposed model (Equation 3) to robustly produce sustained oscillations with the parameters choice described in Section 2.4. Tight parameters identification is usually incompatible with experimental constraints: the robustness of the model is therefore a crucial necessity. To this aim, we systematically modified the amplitude of the synaptic weights distribution (K 12 , K 21 , and K 22 ) and the axonal transmission velocities (c 1 and c 2 ) in ranges of plus or minus 35% their nominal values. Ten values in this range were tested for each parameter, thus corresponding to a total of 10 5 simulations. No stimulation was applied in this protocol (k c = 0). For each simulation, we estimated the main harmonic of the STN and GPe activity, as the maximum peak of the Fourier transform of the signals after transients. The result is summarized by a histogram that illustrates the different oscillations frequencies that can be generated by the model: see Figure 3.
This histogram suggests reasonable robustness to parameter uncertainty as, in all simulations, the main oscillations harmonic either lied within 13-25 Hz or was zero. In other words, in these parameters ranges, oscillations either take place in the β-band or do not take place at all.
In order to assess more tightly the role of each parameters in the sustained oscillations onset, we estimated their frequency by varying only two parameters in a wider range, while leaving all other parameters to their nominal values. It appeared that the amplitude K 22 of lateral GPe synaptic weights may result in sustained oscillations or constant activity, but has little effect on the oscillations frequencies (curves not reported). On the contrary, the synaptic weights between STN and GPe (K 12 and K 22 ) and the axonal velocities (c 1 and c 2 ) do play a crucial role in the frequency of sustained oscillations. Figure 4 shows the frequencies as a function of K 12 and K 21 (left panel) and as a function of c 1 and c 2 (right panel).
It can be observed that each of the considered parameters can be significantly modified without compromising the β nature of sustained oscillations. Nonetheless, too low synaptic weights or two high axonal velocities (meaning lower transmission delays between STN and GPe) may compromise the very oscillations onset. This is in agreement with the analyses FIGURE 3 | Effects of parameters uncertainty on the frequency of sustained oscillations. This histogram illustrates the occurrences of different oscillations frequencies when the parameters K 12 , K 21 , K 22 , c 1 , and c 2 evolve in a range of ±35% their nominal values (see Tables 1, 2). In all 10 5 simulations, when oscillations take place, they lie within the β-band.
conducted on an averaged rate model in Holgado et al. (2010) and Pasillas-Lépine (2013), which clearly showed that sustained oscillations results from the combined effects of two factors: too strong synaptic coupling between STN and GPe resulting from dopamine depletion and axonal propagation delays. This is also in agreement with the basics of control theory, which predicts that instability is favored by high coupling and delays.

Protocol B: Sustained Oscillations Disruption
As seen in the above section, the nominal parameters lead to sustained β-oscillations (approximately 19 Hz) in both STN and GPe. In order to counteract these pathological oscillations, we apply the proposed closed-loop stimulation (Equation 2) from time t = 0.5 s. The results are presented in Figure 5, which clearly shows oscillations disruption. Once stimulation is ON, the remaining oscillations mostly result from the cortical and the striatal stochastic inputs (I 1 (r, t) and I 2 (r, t)).

Protocol C: Poor STN Photosensitization
Depending on the efficacy of transgenesis linked to the viral vector itself, its envelope, its promoter, and to the targeted neurons, it may happen that photosensitization does not affect all STN neurons. It is important to check whether the closedloop stimulation signal still succeeds in disrupting pathological oscillations. To that end, we randomly put to zero 50% of the values of α(r), corresponding to a photosensitization of only half of the STN population (see Figure 2, right panel), while keeping the same value of the feedback gain (k c = 2). The results are presented in Figure 6. They indicate that the proposed closed-loop stimulation signal is robust to poor photosensitization, as the oscillations are still reduced although less efficiently than when all STN cells are successfully photosensitized.
In order to further investigate the link between the level of photosensitization of STN neurons and the efficacy of the proposed closed-loop stimulation strategy, we plot the maximum amplitude of observed β-oscillations as a function of the α(r) degeneracy (meaning the percentage of STN cells that do not respond to light stimulation). The results are presented in Figure 7, for various values of the percentage of STN neurons that do not respond to light stimulation (0, 25, 50, 75, and 100%) and of the stimulation gain k c (2, 6, and 12).
Not surprisingly, the efficacy of the stimulation diminishes as α(r) degeneracy increases. While a degeneracy of 50% has an impact on oscillations disruption for k c = 2 (amplitudes of approximatively 30 sp/s), it can efficiently be tackled by using a larger k c (k c = 6 leads to low-amplitude remaining oscillations).

Protocol D: Single Light Source
The size of the STN impedes the use of several light sources. Although holographic techniques could be envisioned (see Section 4 for a discussion about this), it is crucial to test the validity of the proposed closed-loop stimulation signal when only one light source is available for the whole STN population.  To this end, we applied the following closed-loop stimulation signal: This closed-loop stimulation signal can be seen as a spatial average of the original one proposed in Equation (2) over the whole STN population. Like Equation (2), this averaged stimulation signal relies on measurement of STN activity only. Nonetheless, it owns the advantage of being no longer dependent on position: it is taken proportional to the mean STN activity, and is the same for all STN neurons (modulo the function α(r)).
In this case, a larger feedback gain is required to efficiently disrupt oscillations (k c = 6.5). The results are presented in Figure 8, and show that proper β-oscillations disruption is achieved. We did not succeed in showing this analytically yet. Surprisingly enough, the required amplitude of stimulation is smaller when a single light source is used (compare the right panels of Figures 5, 8).

Protocol E: Effect of Acquisition and Processing Delays
The experimental implementation of the proposed closed-loop stimulation signal in an animal model of Parkinson's disease would require a real-time estimation of the activity in the STN based on the measurements provided by the recording electrodes. This signal acquisition and its processing require some time, which is likely to induce a significant delay in the applied stimulation signal. In order to test the robustness of the proposed strategy to acquisition and processing delays, we consider the following closed-loop stimulation signal: FIGURE 6 | Poor STN photosensitization. 50% of the function α(r) was randomly silenced. Closed-loop stimulation is applied from time t = 0.5 s with gain k c = 2 and target reference z ref (r) = 0. Left: red and black curves illustrate time evolution of the mean activities of GPe and STN, respectively. It is apparent that the attenuation is not as successful as in Figure 5, but oscillations are still attenuated. Right: spatiotemporal evolution of the activity in STN (left) and GPe (right).
FIGURE 7 | Sustained oscillations amplitude vs. photosensitization degeneracy. We measure the efficiency of the proposed closed-loop stimulation through the amplitude of the remaining oscillations. We consider an increasing degeneracy of α(r), corresponding to an increasing number of STN neurons that do not respond to photostimulation. The three curves correspond to different values of the feedback gain k c (2, 6, and 12). The darkest curve corresponds to the highest value and the lighter one to the lowest value. The red dashed line indicates the oscillations maximum amplitude when no stimulation is applied.
where T indicates the time needed to acquire and process the STN activity. The results are presented in Figure 9 for a feedback gain k c = 2.
When T is too large (T ≥ 10 ms), β-oscillations are no longer disrupted (Figure 9, left panel) and can even be amplified by the stimulation, as already observed in Pasillas-Lépine et al. (2013) in an averaged rate model. Nevertheless, if T is reasonably small, oscillations disruption remains satisfactorily achieved (Figure 9, right panel, T = 5 ms). This indicates that strong attention should be paid on the code optimization for data processing in view of experimental implementation. Filtering strategies, such as the one presented in Pasillas-Lépine et al. (2013) and a careful choice of the target reference z ref (r) could probably increase robustness to acquisition and processing delays (not implemented here).
In order to assess more finely the impact of acquisition and processing delays T and their link with the feedback gain k c , we tested 10 different delays (1, 3, 5, 7, 8, 9, 10, 13, 15, and 20 ms). For each delay, we measured the oscillations amplitude after transients for different values of k c (2, 6, and 12), as we did for protocol D. The results are reported in Figure 10.
For small gains (k c = 2), the model is quite robust to acquisition and processing delays: only delays above 9 ms significantly degrade the performance of the stimulation. This robustness to delays decreases as the feedback gain increases. Depending on the delay value, it may happen that stimulation is even counter-productive, in the sense that sustained oscillations get enhanced rather than disrupted.

Discussion
In this theoretical work, we have proposed a model of the STN-GPe network able to take into account both temporal and spatial evolutions of electrophysiological activity, and a closed-loop stimulation strategy to suppress sustained oscillations related to Parkinson's disease and test their causality in future experiments with animal models.  oscillations disruption is lost for too large delays (T = 10 ms). Right: β-oscillations disruption is still efficiently achieved for an acquisition and processing delay T = 5 ms.
Our model is based on delayed neural fields to account for spatiotemporal dynamics. The need to consider spatial heterogeneity in STN and GPe arises from the facts that: (i) the efficiency of DBS requires a fine positioning of the stimulation electrodes (a millimetric scale precision is required, Anheim et al., 2008;Luján et al., 2009), thus indicating a strong positiondependency of the populations sensitivity; (ii) connectivity organization between neurons among a population or between populations cannot be taken into account in averaged models, although some studies underline its crucial role in oscillations onset (Schwab et al., 2013); (iii) the complex functioning of the structures involved, assimilated in Montgomery (2007) to a hierarchical interconnection of multiple nested, loosely coupled oscillators, requires spatially distributed models; (iv) the use of multi-plot electrodes with the most recent technologies permits a spatial resolution below 50 µ m, which cannot be exploited in averaged models; and (v) human recordings indicate a possible role of spatial correlation of STN synchronization in the parkinsonian motor symptoms (Cagnan et al., 2015). We chose to model this spatiotemporal dynamics of the STN-GPe network by neural fields (Coombes, 2005;Bressloff, 2012). Neural fields were originally introduced by Amari (1977) to model cortical layers. They have been extensively used in the study of many different cortical phenomena such as propagating waves (Cremers and Herz, 2002;Bressloff and Webber, 2012;Meijer and Coombes, 2014), synaptic depression (Bressloff and FIGURE 10 | Oscillations amplitude vs. acquisition and processing delays. We measure the efficiency of the closed-loop stimulation (Equation 6), which includes processing and acquisition delays T, through the amplitude of the remaining oscillations after transients. We consider increasing values of T (1, 3, 5, 7, 8, 9, 10, 13, 15, and 20 ms). The three curves correspond to different values of the feedback gain k c (2, 6, and 12). The darkest line corresponds to the highest value and the lighter to the lowest value. The red dashed line indicates the oscillations maximum amplitude when no stimulation is applied.
Kilpatrick, 2011), and neural oscillations (Hutt et al., 2008). As proposed in Holgado et al. (2010) and confirmed in this work, transmission delays between STN and GPe are likely to play a key role in pathological oscillations onset and are thus included in the model. Neural fields affected by axonal propagation and synaptic delays have been studied in Faye and Faugeras (2010); Veltz and Faugeras (2011), which provide the mathematical framework this study relies on. They have been used in (Modolo et al., 2010) for closed-loop stimulation, using a term-cancellation that requires precise knowledge of the dynamics involved, including synaptic weights distributions.
The proposed spatiotemporal STN-GPe model reproduces sustained β-oscillations in a robust way, meaning even for relatively large uncertainty on parameters. We have proposed a closed-loop stimulation policy that provably disrupts sustained oscillations in this spatiotemporal model of the STN-GPe network. Inspired by Pasillas-Lépine et al. (2013), this closedloop stimulation signal is taken proportional to the STN activity. Assuming that STN neurons have successfully undergone photosensitization, we have formally shown that, provided that no endogenous oscillations take place within the GPe, stabilization of the whole network can always be achieved by this closed-loop stimulation if the feedback gain is taken sufficiently large. This result shows that any oscillations generated within the STN or resulting from a pacemaker effect between STN and GPe populations can be tackled in this model. It is worth stressing that the assumptions of this theoretical result do not rely on the precise knowledge of the parameters and functions involved in the model. No precise information on the activation functions are required: an estimate of their maximum slopes is enough. Similarly, the synaptic weight distributions do not need to be known precisely: all we need to know is whether the intensity of the lateral GPe connections is strong or weak. Also, no knowledge is required on the time constants or on the delays. Finally, no precise tuning of the feedback gain is required: it is enough to pick it sufficiently large. These features make our closed-loop stimulation signal intrinsically robust to parameter uncertainties and promising from an experimental testing perspective. The requirement that lateral synaptic gains of the GPe be low enough is in agreement with the experimental findings of Bar-Gad et al. (2003), where functional correlations between GPe neurons have been reported as either weak or non-existent.
The spatial nature of the model allows to conduct more precise in silico experiments exploiting the geometry of STN and GPe and the synaptic distributions of neurons within and between these structures. Our numerical experiments revealed that the model is robust to poor STN photosensitization and that, in case of low photosensitization, a larger feedback gain can be used in order to successfully preserve oscillations disruption.
From an experimental testing in animal models of Parkinson's disease perspective, we stress that the use of optogenetics neuromodulation is crucial since it enables simultaneous neuro-excitation, neuro-inhibition and neural recordings, with biotechnological tools already validated experimentally (Han and Boyden, 2007). Moreover, these numerical findings point toward a necessary trade-off between robustness to poor photosensitization, which is favored by the use of large feedback gains, and robustness to acquisition and processing delays, which is better with small feedback gains. Simulations revealed that the proposed strategy is little robust to acquisition and processing delays if a large feedback gain is used. This issue needs to be carefully taken into account when performing experimentations in animal models of Parkinson's disease. Careful attention should be paid to the optimization of the acquisition algorithms in order to reduce the round-trip-time as much as possible. Future work should be conducted to make the proposed closed-loop stimulation more robust to these delays. This may take the form of a low-pass filter, as in Pasillas-Lépine et al. (2013), or of more advanced strategies inspired from control theory (including predictive features).
For experimental translation, the ability to deliver a spacedependent light for photostimulation is also crucial. Light stimulation is usually delivered through an optical fiber connected to a laser or directly by a light-emitting diode, so that all neurons expressing opsins and receiving photons are simultaneously photo-stimulated, receiving a number of photons depending simply on the output of the fiber, the distance to the fiber and the physical properties of the biological tissue interacting with photons. Patterned illumination constitutes a promising way to illuminate precisely shaped regions of interest, with neighboring neurons receiving potentially different patterns of illumination in intensity and in time (Bovetti and Fellin, 2015). Patterned illumination can be performed by phase modulation of light: liquid crystal spatial light modulators (LC-SLMs) can modulate two-photon light according to arbitrary shapes, thus allowing simultaneous illumination of multiple neurons with different phase maps over time in order to probe neuronal networks with complex spatial and temporal patterns of neuronal activation. For instance, optogenetic neuromodulation with millisecond temporal precision and cellular regulation has been performed in retinal ganglion cells expressing channel-rhodopsins (Reutsky-Gefen et al., 2013). Although we have shown through simulations that proportional stimulation remains efficient even when a single light source is used for the entire controlled population, holographic optogenetics thus constitutes a promising technique to provide precise time and space resolution of the closed-loop stimulation signal.
Future work should deepen the dynamics induced by light stimulation in optogenetics experiments. The model used in this work assumes that the impact of photostimulation is similar to that of a fictitious afferent neural population whose rate and excitatory/inhibitory features can be regulated in real time. In particular, the combined effects of intensity and frequency of light stimulation are not captured by the present model. Tighter modeling could be developed in the future based on existing works such as Nikolic et al. (2009).
Synaptic plasticity and homeostasis are not taken into account in the employed model. Recent studies evoke the role of synaptic plasticity in the long-term effects of DBS (Tass et al., 2012;van Hartevelt et al., 2014). Neural fields models including learning mechanisms have been studied in Galtier et al. (2012); Detorakis and Rougier (2014); Fix (2014) and we believe that the results presented here could constitute a relevant framework to deepen the analysis by including plasticity or homeostasis mechanisms. Furthermore, the analysis conducted here focuses on the STN-GPe interaction. The theoretical result presented here can be readily adapted to more than two neuronal populations by modifying the synaptic coupling distributions. Consequently, our framework is well fitted to assess the possible role of other brain structures in the pathological oscillations onset.
Finally, as stressed in the introduction, the precise mechanisms of parkinsonian oscillations in deep brain structures is still a matter of debate (Bevan et al., 2002;Lang and Zadikoff, 2005). This paper focuses on one possible mechanism, namely the STN-GPe network acting as a pacemaker. This hypothesis is supported by in vitro experiments reported in Plenz and Kital (1999) using mature organotypic cortex-striatum-STN-GPe cultures. To the best of our knowledge, no in vivo experiments has yet confirmed or infirmed this hypothesis, although Nambu and Tachibana (2014) underlines the role played by STN-GPe reciprocal connections. We believe that the theoretical insights provided by the model presented here and its associated closed-loop stimulation constitute interesting bases for further experimental investigations in animal models of Parkinson's disease in order to confront this hypothesis with other possible mechanisms, including cortical patterning of STN (Magill et al., 2001) and striatal origins of oscillations (McCarthy et al., 2011).