Modulation of epileptic activity by deep brain stimulation: a model-based study of frequency-dependent effects

A number of studies showed that deep brain stimulation (DBS) can modulate the activity in the epileptic brain and that a decrease of seizures can be achieved in “responding” patients. In most of these studies, the choice of stimulation parameters is critical to obtain desired clinical effects. In particular, the stimulation frequency is a key parameter that is difficult to tune. A reason is that our knowledge about the frequency-dependant mechanisms according to which DBS indirectly impacts the dynamics of pathological neuronal systems located in the neocortex is still limited. We address this issue using both computational modeling and intracerebral EEG (iEEG) data. We developed a macroscopic (neural mass) model of the thalamocortical network. In line with already-existing models, it includes interconnected neocortical pyramidal cells and interneurons, thalamocortical cells and reticular neurons. The novelty was to introduce, in the thalamic compartment, the biophysical effects of direct stimulation. Regarding clinical data, we used a quite unique data set recorded in a patient (drug-resistant epilepsy) with a focal cortical dysplasia (FCD). In this patient, DBS strongly reduced the sustained epileptic activity of the FCD for low-frequency (LFS, < 2 Hz) and high-frequency stimulation (HFS, > 70 Hz) while intermediate-frequency stimulation (IFS, around 50 Hz) had no effect. Signal processing, clustering, and optimization techniques allowed us to identify the necessary conditions for reproducing, in the model, the observed frequency-dependent stimulation effects. Key elements which explain the suppression of epileptic activity in the FCD include: (a) feed-forward inhibition and synaptic short-term depression of thalamocortical connections at LFS, and (b) inhibition of the thalamic output at HFS. Conversely, modeling results indicate that IFS favors thalamic oscillations and entrains epileptic dynamics.


INTRODUCTION
Deep brain stimulation (DBS) for Parkinson's disease (PD) and other movement and psychiatric disorders-including dystonia, tremor, and depression-is clinically used today as a conventional therapeutic procedure for the alleviation of symptoms (Sillay and Starr, 2009). Since the early 90s, neurologists also attempted to apply DBS to other neurological disorders, typically to intractable epilepsies in order to suppress-or at least dramatically reducethe occurrence of seizures [see recent review in Boon et al. (2009)]. These studies followed early scientific evidence showing potentially beneficial effects of DBS on epileptic neural dynamics in animal models (Reimer et al., 1967;Hablitz, 1976) as well as in patients (Cooper et al., 1973;Davis et al., 1982; Weller, 1983). However, contrary to PD, the optimal "antiepileptic parameters" of DBS for reducing the frequency of seizures are much more variable among patients and the number of nonresponders to stimulation still perplexes scientists. Moreover, in responding patients, the fine tuning of stimulation parameters in a patient-specific manner remains indispensable for maximizing antiepileptic effects. On that account, many fundamental questions are frequently raised: where and when to stimulate, at which frequency, at which current intensity, and with which current waveform?
The answers to these questions remain bound to our current, and still limited, understanding of the mechanisms by which DBS modulates neuronal dynamics, whether normal or pathological. Today, the precise mechanisms of neuronal modulation by DBS remain elusive. In addition, these mechanisms are controversial as observed effects are sometimes opposite (McIntyre et al., 2004b). Among the many studies reported over the last decade, identified mechanisms regarding HFS include: local depolarization blockade by HFS (Beurrier et al., 2001), synaptic depression due to neurotransmitter depletion (Shen et al., 2003;Kim et al., 2012), synaptic inhibition (Filali et al., 2004), disruption of the thalamocortical network's dysrhythmia (McIntyre and Hahn, 2010;Kendall et al., 2011). As far as LFS is concerned, some studies described a transient synaptic depression that alters synaptic transmission (Jiang et al., 2003;Speechley et al., 2007). Finally, IFS is routinely used in the context of presurgical evaluation of patients with drug resistant epilepsy to map epileptogenic and functional brain areas. It has long been observed that this type of stimulation is prone to trigger epileptic afterdischarges (Goddard, 1967). This brief overview shows that the spectrum of involved mechanisms is very large and that distinct stimulation frequencies trigger distinct cellular/network processes. More precise insights into these processes will come with increased knowledge about both biophysical and neurophysiological effects of stimulation currents on underlying neuronal systems.
However, the access to cellular and network mechanisms induced by DBS is rather difficult in animal models of epilepsy and (almost) impossible in patients especially in large-scale systems like the thalamocortical loop. An alternative approach is the use of computational models based on physiological data to first reproduce and then explain changes in cerebral activity as a function of stimulation conditions (stimulation site, intensity, and frequency). This is precisely the objective of this study, with a special focus on the distinct effects of DBS frequency on cortical epileptic dynamics.
Our investigation combines computational modeling and clinical data. We explored stimulation effects in a lumped-parameter mesoscopic neural mass model of the thalamacortical loop, inspired from previously published models (Suffczynski et al., 2004;Lopes Da Silva, 2006;Roberts and Robinson, 2008;Crunelli et al., 2011).
Although these models are lumped representations of underlying neuronal systems, they offer a number of advantages in the context of this study. First, neural mass models include subpopulations of principal excitatory cells and inhibitory interneurons. Second, these models were shown to produce realistic activity as observed in LFPs or EEG under normal (Freeman, 1973;Lopes Da Silva et al., 1974) or epileptic conditions [review in Lytton (2008); Wendling (2008)]. Third, main parameters (mean membrane potential and firing rate) provide access to the investigation of several stimulation-induced (patho)physiological mechanisms. For instance, a neural mass model was successfully used in the context of direct low-intensity pulse stimulation in the hippocampus to explain the behavior of evoked responses during the transition to seizures (Suffczynski et al., 2008).
In particular, using this model, we analyzed the neurophysiological effects induced by direct thalamic stimulation on epileptic cortical dynamics at low frequency (LF, < 20 Hz), intermediate frequency (IF, and high frequency (HF,. Model parameters were tuned to reproduce a typical pathological oscillatory activity observed in a neocortical lesion (focal cortical dysplasia, or FCD) in a patient with drug-resistant epilepsy. Intracerebral EEG (iEEG) signals observed during thalamic stimulation (centromedian nucleus) of this patient revealed particularly pronounced frequency-dependent modulation of the FCD pathological activity. Therefore, this data set offered the unique opportunity to identify key model parameters for which such a frequency-dependent modulation could be reproduced and, subsequently to get insights regarding the mechanisms underlying the modulatory effects, in the FCD, of thalamic stimulation. Results revealed that LFS favors feed-forward inhibition and short-term depression at the cortical level and that HFS inhibits the thalamic activity, while IFS reinforces reticulothalamic oscillations thus entraining cortical pathological epileptic dynamics.

MATERIALS AND METHODS
In this section, we present (1) the neurophysiologically-relevant computational model that we developed to study thalamic DBS, (2) the real depth-EEG dataset used for model tuning and, (3) the signal processing methods used for characterizing real and simulated EEG signals.

MODEL OF THE THALAMOCORTICAL LOOP
In order to study the effects of thalamic DBS on cortical dynamics, we implemented a physiologically-plausible mesoscopic model of the thalamocortical loop. This model accounts for the average activity of both cortical and thalamic compartments which include various types of neuronal populations interacting via synaptic transmission. This modeling approach was first proposed in the early 70s (Wilson and Cowan, 1972) and further enriched in order to interpret electrophysiological recordings and study brain dynamics, in the olfactory (Freeman, 1973) and the thalamocortical (Lopes Da Silva et al., 1974) system, for instance, as well as the dynamics of cortical oscillations (Nunez, 1974). This approach was then developed by other research groups in the context of state changes in brain dynamics (Wright et al., 1985), visual evoked potentials (Jansen et al., 1993), dynamics of the human alpha rhythm (Stam et al., 1999) or pathophysiological mechanisms of ictal transitions in epilepsy (Wendling et al., 2000(Wendling et al., , 2002Suffczynski et al., 2001;Robinson et al., 2002;Liley and Bojak, 2005;Breakspear et al., 2006). Later, neural mass models were also used in studies dealing with the connectivity among cortical regions and the impact of model parameters on the power spectrum of EEG or MEG signals (Robinson et al., 1997;David and Friston, 2003;Zavaglia et al., 2006).

Model architecture
The model architecture was inspired from previously published models of the thalamocortical loop (Suffczynski et al., 2004;Lopes Da Silva, 2006;Roberts and Robinson, 2008;Crunelli et al., 2011). In a global view, the model was built of three interconnected compartments: a cortical compartment, a thalamic compartment, and a reticular compartment, in accordance with previously published models ( Figure 1A) and with anatomical data (Figure 1B). Each compartment includes one or several subpopulation(s) of neurons, either excitatory or inhibitory. Generally speaking, the input/output functions of a considered subpopulation are represented by two mathematical equations that were respectively named "pulse-to-wave" (input) and "wave-to-pulse" (output) by Walter Freeman (Freeman, 1992). The former is a linear transfer function that converts the presynaptic average density of afferent action potentials into an average postsynaptic membrane potential (PSP), either excitatory (EPSP) or inhibitory (IPSP). The output function is a static nonlinear function (sigmoid) that provides the average pulse density of action potentials fired by neurons depending on the sum of EPSPs and IPSPs at the input. This non-linear function accounts for threshold and saturation effects that take place at the somas and initial axonal segments of considered cells.
Formally, the input function is represented by a second order low-pass filter H(s) = W/(s+1/τ w ) 2 (where s is the Laplace variable). The impulse response of this filter is given by Parameters W and τ w are tuned such that h(t) approximates the shape of real excitatory (glutamatergic) or inhibitory (GABAergic) postsynaptic potentials (Lopes Da Silva et al., 1976). The quantity W.τ 2 w is the static gain of filter h. Lumped parameter τ w (expressed in s) is linked to the kinetics of synaptic currents. It determines both the rise time (t rise = τ w ) and the decay time (t decay = 3.146τ w ) of the second order filter impulse response h and it is usually adjusted with respect to the physiological rise and decay times of actual PSPs (Molaee-Ardekani et al., 2010). Given the time constantτ w , parameter W can be used to adjust the sensitivity of synapses (the maximal PSP amplitude is W.e −1 ). An alternative implementation of the h function was introduced in Bojak and Liley (2005) and is described in detail in Molaee-Ardekani et al. (2013). It is based on a bi-exponential pulse-to-wave function with two time constant parameters. This implementation allows for the separate adjustment of the rise and decay times of PSPs, and therefore a better approximation of actual PSPs in some circumstances. Besides, the output function , where 2e 0 is the maximum firing rate, v 0 is the postsynaptic potential corresponding to a firing rate of e 0 and r is the steepness of the sigmoid.

The cortical compartment
The cortical compartment was inspired from an existing model of the neocortex which proved its capability of generating both normal and epileptiform activity. Readers may refer to Molaee-Ardekani et al. (2010) for details. In brief, the cortical compartment integrates a subpopulation of pyramidal cells (P, W = A C , τ w = τ ac in Equation 1) and two inhibitory neuronal populations (I c 1 and I c 2 , Figure 1A) representing soma-and proximaldendrite targeting interneurons (GABA A, fast currents, W = G C , τ w = τ gc in Equation 1) and dendrite-targeting interneurons (GABA A, slow currents, W = B C , τ w = τ bc in Equation 1), respectively. Pyramidal collateral excitation was implemented as in Jansen et al. (1993).
In addition, these three cortical subpopulations receive excitatory input from the thalamic compartment. Therefore, feedforward inhibition (FFI) is represented in the model as the two subpopulations of interneurons project to the pyramidal subpopulation (see The Thalamic and Reticular Compartments paragraph below).

Short-term synaptic depression (STD)
STD is present in the neocortex (Boudreau and Ferster, 2005). It can be potentially involved in the context of direct stimulation of the thalamus as TC cells directly project to cortical pyramidal cells. Consequently, this mechanism was implemented at the interface of thalamic/cortical compartments. To our knowledge, an implementation of STD in neural mass models has not been proposed before.
In our model, we represented a modulatory effect of the amplitude of the average EPSP (parameter A C ') at the level of subpopulation P depending on the density of action potentials [d AP (t)] coming from the thalamic compartment. This modulatory effect was obtained by multiplying A C ' by a time-varying coefficient κ(t) ∈ [0.6, 1] where the function describing the evolution of κ(t) was derived from Chung et al. (2002). This study shows that: (i) cortical EPSPs drop by 40% under periodic lowfrequency intense thalamocortical (TC) cell firing and, (ii) this drop in cortical EPSP is directly linked to transient depression of thalamocortical monosynaptic projections to pyramidal neurons.
In line with these observations, STD was implemented as fol- AP is low-pass filtered (cutoff frequency = 10 Hz) to restrict the STD effect to LFS. Then, from each time t η at which the filtered signal d f AP (t) exceeds a firing rate equal to η, the κ(t) coefficient undertakes an exponential decay given by and where t − η is the time instant that just precedes t η . The decrease of κ(t) is limited to the time interval [t η + 0.45 s] and cannot exceed 40%, total. Parameters η and τ were set to 0.8 and 8 s, respectively.

The thalamic and reticular compartments
The thalamic compartment was limited to one population of excitatory neurons (known as glutamatergic thalamocortical -TC -cells) receiving glutamatergic EPSPs (W = A Th , τ w = τ aTh in Equation 1) from cortical pyramidal cells (P) and GABAergic IPSPs with slow (W = B Th , τ w = τ bTh in Equation 1) and fast (W = G Th , τ w = τ gTh in Equation 1) kinetics from the reticular compartment (RtN). Here, we increased the time constant (τ bTh ) with respect to τ bc to account for both GABA A, slow -and GABA Breceptor mediated currents in a single variable. TC cells directly target both cortical pyramidal cells and interneurons. The activation of these GABAergic interneurons subsequently promotes inhibition of pyramidal cells after a di-synaptic delay. Therefore, TC cells activation induces first an EPSP followed later on by an IPSP on cortical pyramidal cells, resulting in feed-forward inhibition (FFI). The RtN compartment comprised two inhibitory subpopulations, namely I RT 1 and I RT 2 which both receive excitatory input from the cortical (W = A Rt , τ w = τ aRt in Equation 1) and the thalamic (W = A Rt , τ w = τ aRt in Equation 1) compartments.

Simulation of stimulation effects
Stimulation currents induce a perturbation of the membrane potential of neurons. At cellular level, this effect can be accounted for by the "λE model", which is well grounded in the biophysics of compartment models (Rattay, 1998;McIntyre et al., 2004a;Manola et al., 2005Manola et al., , 2007) (see Miranda et al., 2009 for a review) and supported by in vitro experiments (Bikson et al., 2004;Frohlich and McCormick, 2010). This model V ≈ λ. E approximates the membrane potential variation V as a linear function of the electrical field E induced by stimulation ( λ representing the membrane space constant). In our neural mass model, the situation is less straightforward as space is not explicitly represented, conversely to detailed or mean-field models. However, within a certain range of intensity values, it has been shown that the membrane potential variation V is modified in a linear way with respect to the electrical field which is itself proportional to the stimulation intensity (Bikson et al., 2004). These considerations led us to also assume a linear variation for the mean membrane potential as a function of stimulation intensity, in stimulated sub-populations of neurons. In addition, stimulation was represented by a train of periodic monophasic depolarizing pulses. The pulse width was fixed to 1 ms (as in clinics). Pulses were low-pass filtered to account for the average time of repolarization (set to 4.8 ms) in stimulated sub-populations of cells. The resulting stimulation signal was added to the mean membrane potential of neuronal sub-populations included in the thalamic (TC) and reticular (I RT 1 and I RT 2 ) compartments of the proposed model. The depolarizing effect was weighted by three coefficients S TC , S Rt1 and S Rt2 ( Table 1) accounting for the possibly different stimulation impact at the thalamic and reticular level.

Model parameters, outputs, and implementation
Parameter values as well as physiological interpretation are provided in Table 1. Note that each synaptic connection in the model is weighted by a connectivity constant denoted by C SP1−SP2 where SP1 and SP2, respectively, denote the source and target subpopulations. In addition, two Gaussian noise inputs p P (t)∼ N(μ P , σ P ) and p TC (t)∼ N(μ TC , σ TC ) were used to represent nonspecific inputs on pyramidal and thalamocortical cell subpopulations. Finally, signals simulated at the level of pyramidal cells in the cortical compartment and at the level of TC cells in the thalamic compartment were chosen as model outputs. They correspond to the sum of PSPs at each compartment respectively. The temporal dynamics of these signals provide a good approximation of actual LFPs. The model was implemented in Simulink®, and all other complementary scripts were implemented in MATLAB®.

REAL DATA FOR MODEL TUNING
We used real clinical data to tune the model into a functioning mode which simulates pathological activity. The clinical data set was limited to a unique patient who underwent thalamic DBS during the presurgical intracerebral EEG exploration (iEEG performed with depth electrodes implanted under stereotaxic conditions) at the Epilepsy Surgery Unit, Rennes University Hospital. This particular patient was chosen for two main reasons: (1) the pronounced frequency-dependent stimulation effects observed during his preoperative diagnostic iEEG exploration at LF, IF and HF in addition to (2) the existence of an epileptogenic zone in a limited area of the premotor cortex (PMC).
In brief, this patient suffered from partial drug-resistant epilepsy since the age of two. MRI scans and EEG recordings pointed out the existence of a neuronal malformation known as FCD in the PMC at the origin of seizures. This type of cortical malformation is known for its epileptogenic features like neuronal hyperexcitability and hypersynchronization and its characteristic epileptiform discharges (continuous, rhythmic or semirhythmic spikes, and polyspikes) (Avoli et al., 2003;Palmini, 2010) as shown in Figure 2C. Based on various clinical studies reporting the modulation of epileptic cortical activity by the Parameters of the nonlinear sigmoid function (transforming the average membrane potential to an average density of action potentials)  stimulation of the CM nucleus (Velasco et al., 1995(Velasco et al., , 1997(Velasco et al., , 2000(Velasco et al., , 2001(Velasco et al., , 2007, it was decided by neurologists and neursosurgeons to implant a depth electrode in this nucleus, as potentially beneficial for the patient who gave his informed consent. During the presurgical exploration, the stimulation of the thalamic CM nucleus (CMN) induced frequency-dependent modulation of the pathologic activity of the FCD (Figure 2). Readers may refer to (Pasnicu et al., 2013) for detailed information. Interestingly, LFS (2 Hz, 4 mA) and HFS (70, 100, and 150 Hz, 0.8 mA) desynchronized the pathological activity of the FCD, while IFS (50 Hz, 0.8 mA) barely affected it. These segments of signals corresponding to either typical pathological activity or modulated activity (depending on stimulation conditions) were used to optimize the model parameters.

PROCESSING OF REAL AND SIMULATED SIGNALS
The use of signal processing techniques was necessary (i) to quantify the above-described effects of stimulation in real iEEG signals, and (ii) to define a feature-vector-based cost function for model parameter optimization. Figure 3A illustrates the feature extraction methodology. iEEG signals recorded in the FCD in absence of stimulation (LFPs FCD ) and under different stimulation conditions were decomposed using an orthogonal matching pursuit algorithm [matching pursuit toolkit-MPTK- (Krstulovic and Gribonval, 2006)]. First introduced in 1993 ( Mallat and Zhifeng, 1993), matching pursuit is signal processing algorithm used to decompose any time series into a linear sum of waveforms selected from a predefined dictionary based on a mother wavelet. To proceed, a proper multi-scalar dictionary of Gabor, Fourier, and Dirac atoms was first defined to account for real iEEG signal components (time-frequency atoms are waveforms well localized in both the time and the frequency domains). In line with (Krstulovic and Gribonval, 2006), the multi-scalar dictionary was formed by translation in time and amplitude/frequency modulation of atoms (defined as Gabor and Fourier functions in our case), over ten different user-defined time scales (i.e. the atom durations, ranging from 0.125 to 5 s). Then, the algorithm provided a table of time-frequency parameters associated to the detected atoms (atom type, central frequency, phase, scale, amplitude, position). Identified atoms were reconstructed using the extracted parameter table and their analytical expression. They were then associated to a given frequency band depending on their central frequency. These frequency bands corresponded to the classical EEG bands as defined in normal adults (δ 1 [0 -  (Figure 3A, blue). Finally, a 9D feature vector V F was defined from the normalized energy distribution in these frequency bands, itself computed as the sum of averaged (over time) atom energies relative to the total signal energy (Figure 3A, green).

MODEL OPTIMIZATION UNDER THE "NO STIMULATION" CONDITION
In order to simulate LFPs FCD , we optimized the excitation/inhibition ratio of the cortical compartment. Thus, the average EPSP/IPSP amplitude parameters of the cortical compartment {A C , B C , G C } were considered as free parameters while all other model parameters were set to fixed values ( Table 1). The optimization method is illustrated in Figure 3B. distance d(V F, real , V F, model ) when parameters A C , B C , and G C span pre-defined ranges of values according to a Brute-Force procedure.

RESULTS
In this section, results regarding the identification of cellular mechanisms underlying the modulation of cortical activity by thalamic DBS are reported. First, the model capability to reproduce signals similar to those recorded from the FCD in the patient was assessed, under two conditions (no stimulation and during stimulation). Three mechanisms contributing to frequency-dependant stimulation effects could be identified. Then, simulations were performed to analyze the marginal or joint contribution of these mechanisms at low, intermediate or high frequency stimulation.

SIMULATION OF LFPs FCD UNDER NO STIMULATION CONDITION
As a first step, we verified the ability of the model to generate signals that resemble those recorded from the FCD in the considered patient (LFPs FCD ). This procedure, described in sections Processing of Real and Simulated Signals and Model Optimization Under the "No Stimulation" Condition, led us to identify a minimal distance (Figures 4A-C) and thus an optimal parameter vector AC, BC , GC = {6, 14, 16.5} for which simulated signals under the no stimulation condition have similar features as compared with those of real signals ( Figure 4D).

SIMULATION OF LFPs FCD UNDER STIMULATION CONDITIONS
Actual LFPs FCD recorded at various stimulation frequencies (2, 50, 70, 100, and 150 Hz) were first characterized using the matching pursuit method described in section Processing of Real and Simulated Signals ( Figure 3A). Results are shown in Figure 5A where feature vectors of segments of LFPs FCD are represented in a 3D space where axes correspond to merged typical EEG frequency bands (δ 2 to θ 1 , θ 2 to β 1 , β 2 to γ). Results show that the distribution of points in the 3D frequency space is not random but clustered, indicating that the frequency content of LFPs FCD segments depends on the stimulation frequency. In addition, some clusters are very close. This is typically the case for i) the no stimulation (yellow) and the 50 Hz stimulation conditions (red) on the one hand, and ii) the 70 Hz (violet) and 150 Hz (cyan) stimulation conditions on the other hand.
To go beyond the qualitative clustering performed by visual inspection of 3D plots, a K-means clustering algorithm implemented in MATLAB and using a Mahalanobis distance was used to automatically detect the three types of stimulation effects. Initial centroids were randomly chosen. The optimal clustering that globally minimizes intra-cluster inertia is presented in Figure 5B. LFPs FCD segments were automatically classified into three subgroups. The first subgroup contains LFPs FCD segments obtained for low-frequency stimulation (LFS). The second subgroup gathers all segments recorded for high frequency stimulations (HFS, > 70 Hz). And finally, in the third subgroup, segments obtained under the no stimulation and the intermediate stimulation frequency (IFS, 50 Hz) conditions are merged together, suggesting that this stimulation frequency does not reduce the "epileptiform aspect" of the activity reflected in the LFP. Based on this characterization of local field potentials recorded in the FCD (LFPs FCD ), parameters S TC , S Rt1 and S Rt2 were manually tuned to lead the model to generate simulated signals which have spectral characteristics similar to those of actual LFPs FCD . Such a manual procedure was sufficient to reproduce stimulation effects observed in one patient. However, extending the study to a larger group of patients would have made imperative an automated parameter fitting procedure based on the spectral characteristics of real EEG signals as in Rowe et al. (2004). Figure 5B shows the projection of representative simulated LFPs FCD in the 3D frequency space ("M" triangles). As depicted, simulated signals obtained for LFS, IFS and HFS were close to corresponding clusters obtained from real signals for the exact same computation of feature vectors. Shown in Figure 5C, these representatives simulated LFPs FCD do not perfectly match actual signals. However, qualitatively similar bifurcations were observed in the model when the stimulation conditions are changed. Indeed, under the no stimulation (NS) and the IFS condition the model generates rhythmic slow oscillations (δ) with superimposed faster activity (β, γ), as observed in real data. For LFS and HFS conditions, strong modulation of this activity was also obtained in the model.

FIGURE 5 | Characterization and classification of real and simulated data. (A) 3-dimensional (3D) projection of feature vectors (V F, real )
corresponding to different stimulation conditions. This projection was obtained by summing some coordinates of initial 9D feature vectors to get 3D vectors [E(δ 2 )+E(θ 1 ), E(θ 2 )+E(α 1 )+E(α 2 )+E(β 1 ), E(β 2 )+E(γ)]/ [E(δ 1 )+ · · · + E(γ)]. Each vector was then represented by a point in the 3D space (δ 2 +θ 1 , θ 2 + α 1 + α 2 + β 1 , β 2 + γ).Three main classes can be visually identified. At LFS, in the model, the slow wave activity was strongly reduced but spike events occurred in the signals at the instant times of stimulation, mimicking, to some extent, comparable events also present in actual LFPs FCD . Finally, at HFS, slow oscillations (δ) were abolished in the model which generates quasi-normal background activity. This simulated activity was also comparable to real activity observed for HFS stimulation but disclosed less γ activity. Note that these are the effects which were quantified in Figure 5B. The qualitative optimization procedure of parameters S TC , S Rt1 , and S Rt2 was then complemented by an evaluation of parameter sensitivity aimed at studying the impact of random changes affecting the parameter vector = {A C , B C , G C , A Th , B Th , G Th , A Rt } on simulated signals. Parameter vector determines the excitability properties in the three model compartments. As shown in Figure 6, results show that the simulated signals obtained under the four stimulation conditions (NS, LFS, IFS, HFS) stay "quite robust" (in the sense that waveforms are conserved) when parameters stay in the range [ ± ζ. ] with 0 ≤ ζ ≤ 0.2.

MECHANISMS UNDERLYING FREQUENCY-DEPENDANT STIMULATION EFFECTS
Three main mechanisms implemented in the model are required to mimic actually observed effects of the CM nucleus stimulation. These mechanisms are the following: (i) the presence of feed-forward inhibition (FFI) at the level of thalamic projections to the FCD, (ii) the presence of short-term depression (STD) at the level of the thalamocortical glutamatergic synapses and, (iii) the depolarization of RtN inhibitory interneurons targeting TC cells. This result raises an additional question: to what extent the joint effect of these mechanisms is necessary to reproduce frequency-dependant stimulation effects (LFS, IFS, and HFS). In order to assess their individual contribution, we performed simulations where each mechanism was either present in-or removed from-the model (the model parameters remaining unchanged). Results are displayed in Figure 7. First, they confirmed that both FFI and STD mechanisms are jointly necessary in the model to suppress the epileptic activity in the FCD when LFS is being used since the withdrawal of either STD or FFI leads the model to generate epileptic activity at LFS. Second, results indicated that the RtN inhibitory interneurons targeting TC cells (both I RT 1 and I RT 2 subpopulations) must be affected (i.e., depolarized) by the stimulation to obtain a suppression of epileptic activity when HFS is being used, as observed in the patient. Third, and interestingly, an unexpected effect was observed at IFS when the depolarization of I RT 2 interneurons was removed from the model. Indeed, epileptic activity was abolished in this case, which is really unlikely to occur during actual stimulation as both subtypes of neurons are expected to be affected by the direct stimulation of the CM nucleus. These results were complemented by a deeper analysis of the thalamic output (i.e., the firing rate of TC cells) in response to stimulation at low, intermediate and high frequency. Results are provided in Figure 8. First, they showed that the thalamic output dramatically differs depending on the stimulation frequency ( Figure 8A). Under the no stimulation condition, the firing rate continuously oscillates around a certain value (referred to as , Figure 8A). At LFS, the firing rate was found to be lower, except at the stimulation times where it abruptly and transiently increased. At IFS, a balance was observed between time intervals for which the TC firing is above and below . Finally, at HFS, the output of TC cells was found to be very low, i.e., systematically under the threshold . From these observations, we could define (i) two time intervals, 1 and 2, for which the TC cells firing rate is either below ( 1) or above ( 2) and (ii) a "high to low firing" ratio (HtoLR) which provides an indication on the amount of time the TC cells spend firing (up state) relatively to the amount of time they do not fire (down state). Figure 8B provides the evolution of the HtoLR when the stimulation frequency is progressively changing from 0 to 150 Hz in the model. As depicted, these simulations indicated that three stimulation frequency ranges have dramatic effects on the firing of TC cells. First, from 0 to 20 Hz, the down state is predominant. Then, an abrupt jump was observed around 22 Hz indicating that beyond this value, the firing rate dramatically increased. Interestingly, from 55 to 65 Hz, a progressive decrease of the HtoLR was observed. Then, after 70 Hz, the ratio is equal to zero indicating that TC cells did not fire anymore. Finally, in order to relate the thalamic activity with the cortical activity, we plotted the phase portraits (TC cell firing vs. cortical LFP) as illustrated in Figure 8C. Results confirmed the visual inspection of signals simulated at the two sites. For the no stimulation (NS) and for the intermediate frequency stimulation (IFS) conditions, phase portraits were found to be quite similar. They indicated the presence of mixed slow/fast oscillations in both signals. For the low frequency stimulation (LFS) condition, oscillations in the simulated LFP in the FCD were reduced. They came along with short-duration, abrupt and rhythmic augmentations of the TC firing corresponding to stimulation pulses. Finally, for the high frequency stimulation (HFS) condition, oscillations in both types of activity stayed confined to small amplitude values.

DISCUSSION
We modeled the thalamocortical loop in order to investigate frequency-dependent effects of electrical stimulation performed in the thalamus and aimed at modulating the neocortical activity. We chose to elaborate our model at a mesoscopic level, i.e., intermediate between microsocopic and macroscopic.
Regarding the model architecture, we followed a similar approach to that used in previously proposed models of the thalamocortical loop (Robinson et al., 2002;Suffczynski et al., 2004;Breakspear et al., 2006;Roberts and Robinson, 2008;Marten et al., 2009;Crunelli et al., 2011). Our model includes three main compartments: cerebral cortex, reticular nucleus and thalamic relay. Subpopulations of neurons and interneurons located in these three structures interact via excitatory and/or inhibitory synaptic connections. The novelty with respect to aforementioned studies is threefold. First, we modified the cortical compartment in order to better approximate the temporal dynamics of epileptic signals recorded in the FCD. This modification consisted in the use of two types of interneurons (mediating GABAergic IPSPs with slow and fast kinetics on cortical principal cells), as reported in a previous study (Molaee-Ardekani et al., 2010). Second, our model accounts for the direct effects of electrical stimulation. At this stage, we used the V ≈ λ. E assumption according to which the perturbation of the mean membrane potential of neurons is a linear function of the electrical field magnitude induced by bipolar stimulation. This "λE" assumption was already used in neural mass models in the context of low-intensity direct hippocampal stimulation to anticipate seizures (Suffczynski et al., 2008) as well as in the analysis of the stimulus-response relationship of DBS in healthy animals (Adhikari et al., 2009). However, it is worth mentioning that in our model, the three subtypes of neurons (TC cells and both subpopulations of inhibitory neurons in the RtN) are depolarized by the stimulation, as suggested in (Molaee-Ardekani et al., 2013) and conversely to (Adhikari et al., 2009) where only principal cells are impacted. And third, our model includes two well-known mechanisms at the cortical level: feed-forward inhibition (FFI) and short-term depression (STD).
As in any modeling approach, our approach has some limitations. First, the chosen modeling level does not allow for analyzing sub-cellular mechanisms involved in stimulation-evoked changes. Similarly, it does not account for direct activation of axons by stimulation versus somatic inhibition (McIntyre et al., 2004b), nor for the mechanisms of orthodromic/antidromic propagation of action potentials due to stimulation (Degos et al., 2005;Hammond et al., 2007;Dorval et al., 2008). Second, a strong assumption in the type of model we used (neural mass) is related to the intrinsic synchronization among neurons included in a given sub-population. This assumption does not allow for representing either de-or weakly-synchronized firing patterns that may be observed during epileptic activity, in particular during high frequency oscillations that can be encountered in FCDs (Brázdil et al., 2010). Nevertheless, we could accurately reproduce the abnormal rhythms generated in the FCD suggesting that main pyramidal cells have a relatively synchronized activity in this epileptogenic tissue. Third, regarding plasticity-related mechanisms, we only implemented short-term effects (i.e., STD) and neglected long-term plastic changes that may be induced by DBS (Shukla et al., 2013).
Despite these limitations, we could identify a number of mesoscopic factors which could explain the frequency-dependent mechanisms of thalamic stimulation. The model was tuned using electrophysiological data recorded in a patient in whom the centromedian nucleus (CMN) stimulation was particularly efficient to reduce the epileptic activity of a FCD located in the premotor cortex, in a frequency-specific manner. The main findings are summarized in Figure 9.

"NO STIMULATION" (NS) CONDITION
In the model, under the NS condition, excitation among pyramidal cells had to be increased and inhibition had to be reduced in the cortical compartment for producing "pathological" oscillatory rhythms, as observed in the FCD. The thalamocortical loop was found to be responsible for these pathological dynamics, characteristic of FCDs. These findings are in line with histological studies showing that these typical oscillations are generated in altered brain tissue, where inhibition is partially deteriorated or dysfunctioning (Calcagnotto et al., 2005), and where excitation is heavily increased (Avoli et al., 2003). In addition to neuron alterations in the dysplastic tissue (Sisodiya et al., 2009), FCD keeps sufficient projections to-and input from-other brain structures to propagate pathological dynamics (Avoli et al., 2003). As mentioned, the presence of connections with subcortical structures was a necessary condition in the model for producing pathological oscillations resembling those actually recorded in the FCD ( Figure 9A).

LOW-FREQUENCY STIMULATION (LFS) CONDITION
For the low-frequency stimulation (LFS, f < 20 Hz) condition, two mechanisms were found to play a major role for the abortion of epileptic activity in the FCD: short-term depression (STD, i.e., decreased excitatory synaptic efficacy in thalamus-to-cortex connections) and feed-forward inhibition (FFI, i.e., excitation of inhibitory cortical interneurons by TC cells) (Figure 9B). STD was reported in previous studies concerning cortical adaptation to thalamic stimulation, and suggesting that electrical LFS of TC cell axons in vivo resulted in a 40% reduction in cortical EPSPs (Chung et al., 2002). In the same context, LFS trains in adult anaesthetized rats provoked transient long-term depression of thalamocortical synapses; this was measured by up to 40% drop in cortical EPSPs after LFS trains and under the effect of GABA antagonist (Speechley et al., 2007).
As mentioned above, the LFS effects could not be reproduced by the model without incorporating also FFI. Actually,

FIGURE 9 | Frequency-dependent mechanisms underlying DBS. (A)
Under the no stimulation (NS) condition, the thalamocortical loop is responsible for pathological oscillatory rhythms observed in the FCD. (B) For low-frequency stimulation (LFS), feed-forward inhibition (FFI, i.e., excitation of inhibitory cortical interneurons by TC cells) and short-term depression (STD, i.e., decreased excitatory synaptic efficacy in thalamus-to-cortex connections) was found to play a major role for the abortion of epileptic activity in the FCD. thalamocortical ascending fibers directly target pyramidal neurons as well as cortical GABAergic interneurons inducing EPSPs in both cell types (Pouille and Scanziani, 2001). In the model, while less efficient (STD) thalamic EPSPs arrive directly onto pyramidal neurons, IPSPs induced by thalamic stimulation also arrive on pyramidal neurons (FFI) lagging by 1-2 ms. This short latency between the onset of thalamocortical excitation and the onset of feed-forward inhibition presents a temporal "window of opportunity" for pyramidal cells to integrate excitatory and inhibitory inputs, thus keeping the transmembrane potential below firing threshold. In the literature, neuroanatomical and neurophysiological studies (Isaacson and Scanziani, 2011) showed the functional importance of FFI in regulating cortical dynamics by controlling cortical excitability (Gabernet et al., 2005). Our study suggests that LFS regulates cortical excitability by a dual mechanism of FFI and STD ( Figure 9B).

INTERMEDIATE-FREQUENCY STIMULATION (IFS) CONDITION
For the intermediate-frequency stimulation (IFS, 20 < f < 70 Hz) condition, results indicated that the thalamic output is reinforced (increase of TC cells firing) and leads to an increase of the average excitatory post-synaptic potential (EPSP) on cortical pyramidal cells ( Figure 9C). This effect corresponds to an increase of the spatiotemporal summation of unitary EPSPs. In this case, both the cortical excitability and the gain in the excitatory thalamocortical loop is increased, leading to "no anti-epileptic" effect. We did not find much studies using DBS stimulation in the intermediate frequency range of (20-60 Hz) in the context of epilepsy. Nevertheless, it is noteworthy that 50 Hz stimulation frequency is classically used during the presurgical evaluation of patient with intractable partial epilepsy in order to trigger seizures and delineate the epileptogenic zone (Talairach et al., 1974;Jayakar et al., 1992). The same frequency range is also known to provoke afterdischarges and was actually used in the kindling model of epilepsy (Goddard, 1967;Racine, 1972).

HIGH-FREQUENCY STIMULATION (HFS) CONDITION
Finally, for the high-frequency stimulation (HFS, f > 70 Hz) condition, the direct and sustained excitation of reticular nucleus (RtN) interneurons leads to strong inhibition of TC cells and thus to dramatic decrease of their firing rate. Despite the fact that TC neurons are also affected by stimulation, the response of reticular GABAergic neurons to stimulation and the higher efficiency of GABA-mediated currents ensure that IPSPs override EPSPs on TC cells. In this case, the reduced excitatory input to cortical pyramidal cells also leads to a suppression of epileptic activity (Figure 9D). This result corroborates reported stimulation studies where HFS (>100 Hz) was associated with significant decrease in epileptiform discharges in vitro, and reduction in seizure frequency in responding patients (Velasco et al., 2006;Fisher et al., 2010). This hypothesis is in line with recent findings suggesting that HFS of the globus pallidus (GPi) in dystonia patients decreased its firing by stimulation-evoked GABA release from afferent fibers and thereby the enhancement of inhibitory synaptic transmission by HFS (Liu et al., 2012). Similarly, HFS (100 Hz-130 Hz) of the STN neurons in vitro showed a suppression of the activity of the majority of neurons by the reinforcement of inhibitory responses (Filali et al., 2004). Other HFS studies also provided evidence on the inhibition of GPi output during HFS in human patients (Dostrovsky et al., 2000) as well as the disruption thalamocortical network's dysrhythmia (McIntyre and Hahn, 2010;Kendall et al., 2011).

CONCLUSION
In epilepsy research, it is well-admitted that there is, unfortunately, a lack of tangible results regarding the effects of electrical stimulation in the brain. Therefore, the very crucial issue of choosing the "optimal" stimulation parameters remains unsolved, whatever the stimulation procedure. Although computational models are always based on a number of simplifying assumptions, we think that they provide an efficient framework to (i) account for the many and essential factors that may intervene during stimulation procedures and (ii) analyze the links between these factors in a formal manner. This approach is particularly fruitful when models are well grounded in experimental/clinical data (Wendling et al., 2012). This is somehow a weak point of this study since we could make use of data sets recorded in one patient only. However, it should be mentioned that these very informative data sets stay relatively rare since many conditions have to be met (patient candidate to surgery, FCD, electrodes positioned in appropriate structures).
At this stage, the face value of the model is satisfactory. The next step is obviously to test the model predictions using animal models. Experiments can be undertaken in rodents with electrodes implanted in the cerebral cortex and in the thalamus. First, we could start with control animals to assess the modulation of cortical rhythms during/after direct thalamic stimulation at various frequencies and for controlled vigilance states (sleep, awake, resting, exploratory). In these controls, some drugs can be used to alter some parameters related to synaptic transmission (in a more or less specific manner) which have a correspondence in the model, on the other hand. Then, refined experimental models could be introduced to get closer to the epilepsy context including models of developmental dysplastic lesions [see review in Schwartzkroin and Wenzel (2012)]. Hopefully, this combined computational/experimental approach will help us to disclose some of the highly intricate effects of DBS either at local or at network level.