Recovery of Dynamics and Function in Spiking Neural Networks with Closed-Loop Control

There is a growing interest in developing novel brain stimulation methods to control disease-related aberrant neural activity and to address basic neuroscience questions. Conventional methods for manipulating brain activity rely on open-loop approaches that usually lead to excessive stimulation and, crucially, do not restore the original computations performed by the network. Thus, they are often accompanied by undesired side-effects. Here, we introduce delayed feedback control (DFC), a conceptually simple but effective method, to control pathological oscillations in spiking neural networks (SNNs). Using mathematical analysis and numerical simulations we show that DFC can restore a wide range of aberrant network dynamics either by suppressing or enhancing synchronous irregular activity. Importantly, DFC, besides steering the system back to a healthy state, also recovers the computations performed by the underlying network. Finally, using our theory we identify the role of single neuron and synapse properties in determining the stability of the closed-loop system.


Introduction
Open-loop brain stimulation has emerged as a common tool to restore aberrant neuronal activity. The most successful example is the application of high-frequency deep-brain-stimulation (DBS) used to ameliorate motor symptoms in Parkinson's disease (PD). However, even in this case the stimulation induces side-effects such as gait imbalance, cognitive impairment, speech impairment, depression etc. [1]. The main cause of these side-effects is likely to be the constant stimulation, but additional explanations are plausible, e.g. the inability of open-loop stimulation to recover the original computations carried out by the impaired brain area. Thus, there is a clear need for more sophisticated brain stimulation schemes [2][3][4].
Moreover, to exploit the full potential of external brain stimulation as a research and therapeutic tool it is important to obtain theoretical insights that can guide the design of novel stimulation protocols. The goal for these new stimulation methods should ideally be twofold: to alter the dynamical state of the brain activity in a desired manner and to recover the computations performed by the network. Here, we demonstrate that DFC, a method with its origin in chaos control [5,6], can achieve these objectives.
To show that DFC is effective in altering the global activity state, we focus on its ability to switch the network state between synchronous-irregular (SI) oscillatory and asynchronousirregular (AI) non-oscillatory activity. This choice is motivated by the fact that several brain diseases are manifested as a transformation of the AI state to persistent SI oscillations, e.g. in PD [7] and in certain forms of epilepsy [8], or as the inability of the network to generate transient SI activity, e.g. in schizophrenia [9]. To demonstrate that DFC facilitates the recovery of certain types of computations, we also illustrate how a network under DFC can effectively process and route rate as well as temporally coded signals. Thus, DFC not only steers the system to a more physiological activity regime, but it also recovers to a considerable degree the coding abilities of the network as they were present before the onset of the pathology.
Previous theoretical models of closed-loop stimulation are not suitable to study the control of SI oscillations because the dynamics that arise in networks of phase oscillators [10][11][12], in networks of Hodgkin-Huxley neurons [13,14] and in Wilson-Cowan type firing rates models [15] are qualitatively different from the SI oscillations [16,17]. In addition, the physiologically plausible SI oscillations are known to be robust to both noise and heterogeneities [18][19][20] and, therefore, require a more differentiated control approach.
Our control strategy is applicable to any network with arbitrary connectivity that undergoes a Hopf bifurcation and it is useful but not critical to know the parameters a priori. In fact, we show how adaptive tuning methods can be used to estimate the control parameters if the precise network parameters are not known. The DFC based stimulation method proposed here can be in principle applied in all animal models where the use of optogenetic tools allows for direct modulation of the membrane potential. In human patients, where optogenetics is currently not an option, manipulation of the membrane potential can be achieved indirectly via subthreshold electrical stimulation. Importantly, the method is not restricted to stimulation of deep structures, but could be applied non-invasively to modulate activity in cortical layers as well [21]. Finally, the theoretical insights we provide into the mechanisms of feedback control in SNNs could also explain the recent success of event-driven stimulation schemes [22][23][24].
synchronous-regular (SR) oscillations arise when the mean input to the individual neurons exceeds their spiking threshold, resulting in high firing rate and high frequency regular oscillations [16,17]. By contrast, the synchronous-irregular (SI) oscillations arise because of strong synaptic coupling and increased variance of the total input to the neurons. Importantly, the emergence of the SI oscillations is accompanied by a change in the network transfer function and its ability to represent stimulus-related activity. Persistent SI oscillations often are signature of brain diseases, e.g. in PD [7] and epilepsy [8]. The altered network transfer function and the robustness of the oscillations to noise and neuronal heterogeneities pose a serious challenge for stimulation-based therapeutic approaches. In the following we show that DFC is able to both quench SI oscillations and to recover the original network transfer function.
In order to provide an explanatory account of the exact mechanism by which DFC operates we provide specific examples of neuronal networks with known parameters that can generate SI activity. To this end, we used mean-field theory to describe the dynamics of homogeneously and randomly connected recurrent SNNs. The mathematical analysis was corroborated with appropriate numerical simulations of sparsely connected SNNs. We considered networks with only inhibitory and both excitatory and inhibitory neurons.

Control of SI activity in I-I networks
While our goal is to reveal the mechanisms by which DFC controls SI activity in excitatoryinhibitory SNNs, it is more instructive to first demonstrate the concept in a simple, purely inhibitory SNN. In the AI state the population average of the firing rate is constant in time r(t) = r 0 . Therefore, the mean recurrent input that each neuron receives is also constant: where C is the average in-degree, J/C the synaptic coupling strength and s(t) the postsynaptic current. In such a network, the emergence of SI oscillations can be investigated by analyzing the stability of the network firing rate in the AI state [18,19]. A small perturbation in the steady-state firing rate where e λt is an eigenmode of the network dynamics with complex eigenvalue λ, leads to a perturbation in the recurrent input where J is the synaptic coupling strength and S is the synaptic response function. In a recurrent network both perturbations have to be consistent, that iŝ where R is the neuron response function. This results in the self-consistency equation: In a purely inhibitory network J is negative, but here the negative sign of J has been absorbed in the phase S(λ). We can then compute the eigenvalue spectrum, that is the roots λ that satisfy Eq (1). When the eigenvalues have a positive real part, the AI state is unstable and the SNN settles in the SI state. Note that due to the synaptic delays the spectrum is infinite. However, in time-delay systems of the retarded type that we are considering here, the total number of unstable eigenvalues is always finite [25]. Increasing J shifts the spectrum towards more positive values on the real axis. For a critical value J cr a complex pair of eigenvalues crosses the imaginary axis and the system becomes unstable through a supercritical Hopf bifurcation [18] (Fig 1). In the following we consider a SNN in which J > J cr , thus resulting in the emergence of SI oscillations. We aim at designing a controller that can alter the global activity state from SI to AI by placing the unstable eigenvalues back to the left half-plane (Fig 1e).
For the implementation of our DFC controller we made two assumptions. First, we chose the instantaneous population rate to be the output state of the system. This state needs to be continuously monitored by the controller. Based on this state, the controller evaluates an appropriate control signal. The second assumption is that the control signal can be applied via current injection directly to the somata of the neurons. Note that the control signal is identical for all neurons in the network. Thus, the synchrony in the network activity was not decorrelated because each neuron received different input. To include the contribution of DFC the self-consistency Eq (1) needs to be modified: where K is the control gain, d c the control delay and M the control kernel. The roots of the above equation (see methods) yield the range of parameters K, d c that move the unstable eigenvalues back to the left-half plane (Fig 1e), which results in a switch of activity from SI to AI.
To verify the analytical solution we simulated an Erdös-Rényi type inhibitory recurrent network of N = 10,000 sparsely connected leaky-integrate-and-fire (LIF) neurons. The neurons were connected with a connection probability = 0.1. Switching on the DFC with parameters estimated from Eq (2), almost immediately results in suppression of oscillations and in a network state that resembles the AI regime (Fig 2a-2d). The suppression of stochastic oscillations is evident both in the spiking activity of single neurons (Fig 2a) and in the population activity of the network (Fig 2b). The spike count variability and the irregularity of single neuron interspike intervals, estimated by the Fano Factor (FF) and the coefficient of variation (CV) respectively, confirm that under DFC the firing of individual neurons in the network follows Poisson statistics (AI: FF = 1.04, CV = 1.01, DFC: FF = 1.02, CV = 0.99). Moreover, the oscillation index that captures the degree of oscillatory activity (see methods) is in both conditions comparable (AI: P T = 1.47, DFC: P T = 1.45) and significantly smaller than in the SI state (P T = 3). The change in the network spiking activity is also observed in the subthreshold membrane potential of individual neurons (Fig 2c and 2d).

Control of SI activity in E-I networks
Next we demonstrate the applicability of DFC in changing the SI state in recurrent networks of excitatory and inhibitory neurons. To this end we simulated a SNN composed of 8000 excitatory and 2000 inhibitory neurons with Erdös-Rényi type connectivity and a connection probability of = 0.1. As in the I-I network, we used the mean-field approach to derive appropriate values for the parameters to attain an SI state. The self-consistency equation for the coupled EI-network under DFC is given by where J ij , d ij is the synaptic coupling strength and delay from population j to population i and R E (R I ) is the neuron response function of excitatory (inhibitory) neurons. Note that we ignored the recurrent couplings within E and I populations, because we were interested only in the oscillations created by the EI-loop. We implemented DFC by recording the activity of neurons in the inhibitory population while stimulating excitatory neurons. Here again K is the control gain, d c the control delay and M the control kernel. Switching on the controller yielded a near instantaneous transition in the network activity from SI to AI (Fig 2e-2h). In this case the original physiological state that we wanted to recover was characterized by slightly less In a coupled network with more than one population additional possibilities for recording and stimulating neurons exist. For instance, we could both record and stimulate the excitatory population (see below "stability and robustness of control domains"). Our results, however, do not depend on the exact identity of the recorded and stimulated neurons. The oscillation frequency in the SI state of the network was in the beta range, f % 25-30Hz, which is characteristic for PD [7]. This suggests that if PD-associated beta oscillations are caused by a strong coupling between STN and GPe [26,27], then our DFC approach could be used to suppress these beta band oscillations. The results presented here are general and the same approach can be applied to suppress oscillations in other frequency bands as well as long as the oscillations underlie a Hopf bifurcation.

Stability and robustness of control domains
To determine the range of values that led to stable control we fixed the control kernel M, using a box function of width 1 ms, and parametrized the system by the control gain K and delay d c .
For each pair of values we simulated the SNN and computed the oscillation index (Fig 3a and  3c). The (K, d c )-plane shows that a stable control domain exists at 7 ms. That is, an effective control delay of d c,eff = 7 ms yields the maximum stability for the resulting AI state. The semianalytical results from mean-field theory (red contour lines derived from Eqs 2 and 3 for the I-I network and E-I network, respectively), are in good agreement with the numerical simulations (blue shades). The only discrepancy occurs when the difference between synaptic and control coupling is small. In such a scenario it is more difficult to maintain constant rates of the stimulated population and the system may become effectively excitatory leading to rate instabilities [28]. Moreover, fluctuations in the mean input that are ignored in our mean-field approach could also become more important. The analysis of these fluctuations is beyond the scope of this work and will be addressed in a future study.

Differential control
Despite the fact that one stable control domain exists, a compensation mechanism to maintain constant firing rates is required to achieve stable control. In real-life applications a detailed fine-tuning may not always be possible. Therefore, we modified our control protocol and introduced an additional delay term d c2 , thus, effectively feeding into the controller the difference between two time-delayed versions of the population activity. For such differential DFC scheme the control signal is given by (see methods): Differential control has been previously used to control unstable periodic orbits [5,29] and to suppress synchrony in networks with discrete-time neuron models [30]. With the differential control we accounted for the fact that recording neural activity and injecting a control current into the neurons introduces a finite time-delay. Therefore, we used a small but non-zero value for d c2 , i.e. d c2 = 1 ms, which is close to the overall closed-loop delay introduced by current technologies [31,32]. A crucial advantage of differential DFC is that no additional rate compensation is required, because the mean contribution of the control signal vanishes Moving in the control parameter space (K, d c ), therefore, did not affect the firing rates of the neurons. This was reflected in the near perfect overlap of theoretical predictions and numerical simulations of the SNN (Fig 3b). In addition, differential DFC introduced two positive effects on the stability of the control domains: (i) The first control domain was expanded, which amounts to an increase in the robustness in the parameter variation. That is, small deviations from the estimated values of the gain and the delay would not be critical for the stability of the AI state achieved by differential control. (ii) A new stable control domain appeared at t = 23 ms. Thus, with differential control there is an increase of the range of parameters that lead to stability. One stable control domain appears at d = 7 ms. The results from mean field theory (red contour) correctly predict the location and shape of the domain. The discrepancy for K > 200 mV is explained by a deviation of the firing rates between numerical simulations and mean-field theory. (b) Differential DFC in the inhibitory network. The first stable control domain around 7 ms is enlarged. An additional small, but stable control domain appears around d = 23 ms. Moving in the state space does not affect the firing rates and therefore no rate compensation is required. The numerical and analytical results (red contour) are in perfect agreement. (c) Direct DFC in an E-I network. The excitatory population is being stimulated while the activity of the inhibitory neurons is recorded. A stable control domain appears around d = 7 ms which reflects the effective delay from the inhibitory to the excitatory population. (d) Same as (c) but here the excitatory population is both recorded and stimulated at the same time. The theoretical analysis yields a small but stable control domain around d = 15 ms. The location of the domains reflects the larger E-I-E loop (see text). In the numerical simulations this control domain does not arise, because fluctuations, which are ignored by the mean field-approach, DFC also enhanced the robustness of the system to external disturbances, e.g. undesired signals at the controller output, measurement noise etc. This becomes evident when we consider the distance B cr of the complex eigenvalues λ i from the imaginary axis for the main stable control domain at t = 7 ms. A more robust closed-loop system is reflected in higher values of B cr . Differential and direct control yielded B cr diff ¼ maxðReðl i ÞÞ ¼ À257 and B cr direct ¼ maxðReðl i ÞÞ ¼ À224, respectively, clearly revealing a more robust system with differential DFC.
Both direct and differential control were effective in a E-I network as well. The location of the stable control domains depended on the exact implementation (Fig 3c-3e). When the activity of the inhibitory population was monitored while stimulating the excitatory population, the main stable control domain appeared at t = 7 ms (Fig 3c). This location is identical with the purely inhibitory network and reflects the overall delay of the I-E path (I-I loop) in the E-I (I-I) network. Indeed for both the I-E path and I-I loop the effective delay is d eff = 7 ms (see methods). By contrast, when the excitatory population was both recorded and stimulated then the location of the domains shifted to around t = 14 ms reflecting the larger overall delay in the E-I-E loop. (Fig 3d and 3e). Note that in this case the stable control domain for direct control was smaller. The reason is that the size of the stable control domains shrinks for larger delays.

DFC control vs noise injection
In both the I-I and E-I SNNs we applied an identical control signal to all stimulated neurons. That is, we did not disrupt oscillations and decorrelated network activity by injecting different currents to each of the neurons. This is in contrast with a widespread assumption that common input always tends to increase correlations in neural activity [33]. The results from the application of DFC reveal that common input can both increase or decrease correlations in SNNs. It is the timing and the amplitude of the common input that determines the direction in which correlations are affected.
It is important to point out that injection of a control signal is not equivalent to the application of additive noise to the system. To demonstrate this we simulated an I-I network and injected Gaussian noise with the same mean and variance as the control signal to all neurons. This stimulation approach failed to suppress SI oscillations (Fig 4a and 4b) indicating that the temporal structure of the control signal is crucial for successful control. Increasing further the noise intensity, e.g. by a factor of ten, eventually resulted in desynchronization of the activity and in quenching of oscillations (Fig 4d and 4e). However, with such strong external noise the network dynamics is predominantly influenced by the input rather than the recurrent activity. This condition is disastrous from a computational point of view, because any information processing taking place within the stimulated brain region would be severely impaired. To illustrate this we recorded the subthreshold dynamics of ten randomly selected neurons in the network (Fig 4f). The huge fluctuations in the membrane potential under the influence of strong external noise are rather pathological. By contrast, the fluctuations in the case of DFC are comparable to those in the physiological AI regime. quickly destabilize the system. (e) Same as (d) but differential DFC is used. The stable control domain is enlarged showing that differential DFC yields more robust control in coupled E-I networks as well. Here only the first delay d c1 is varied and the second delay is kept constant to d c2 = 1 ms. In all panels regions of stable activity are indicated by shades of blue color (dark shades:higher stability). The dashed lines denote that total network coupling J. doi:10.1371/journal.pcbi.1004720.g003

Recovery of network function
The detrimental effect of strong external noise became even more apparent when we studied the response of the network to incoming stimuli. We examined two scenarios. First, we tested how a series of incoming pulse packets composed of randomly distributed spikes are processed by the SNN. We evaluated the network response by the area under the curve (AUC, see methods) for each of the following network states: AI, SI quenched DFC and SI quenched by noise stimulation. A high AUC value reflects better separability of two conditions. It is evident that the AUC in the AI state and in the DFC condition is close to unity indicating that both conditions are comparable in terms of stimulus separability (Fig 5a and 5b). By contrast, when the SI oscillations were quenched by the injection of strong external noise the AUC dropped significantly. That is, DFC, in contrast to strong noise stimulation, does not impair the ability of the network to detect incoming stimuli. Next we tested how DFC affects temporal aspects of the network response. To this end, we provided external correlated inputs to all stimulated neurons and measured the spike train similarity in the network response. We computed the spike distance D that captures the timeresolved degree of synchrony between individual spike-trains ( [34], see methods). Again DFC did not impair the temporal processing as indicated by a clear separation of the two clusters during baseline D B and stimulation D S (Fig 5c). For external noise, however, the two distributions of values strongly overlapped, showing that aspects of temporal processing as measured by pairwise synchrony are clearly compromised when the SI state is disrupted by open-loop noise injection. These results suggest that processing of incoming signals either locally or by downstream areas is feasible in a DFC scheme, but not in an external noise scheme.

Mechanism of DFC
The above two results clearly demonstrate that DFC has multiple advantages compared to the open-loop noisy stimulation. DFC does not only suppress SI activity steering the network to an AI regime, it also facilitates the recovery of the network's ability to process stimulus related information. From its design it is evident that DFC effectively counteracts the increase in coupling strength, which is one of the main causes for the emergence of SI activity. Indeed, the goal of the DFC design was to move the poles of the system at, or close, to their original positions. Ideally, the stimulation kernel M would match the synaptic kernel S with d c = d and the amplitude of the control gain K would be tuned to match the pathological increase of the coupling strength ΔJ. If this were the case, DFC would completely eliminate the effects on the mean recurrent input. This is evident if we consider the modulation to a perturbation in the average input to a neuron That is, under DFC the effects of ΔJ are not visible in the perturbed current term. In practical applications, of course, a perfect match between the control parameters (K,d c ,M) with the synaptic values is not feasible, because the exact shape of the synaptic kernels are not known a priori and have to be estimated. Nevertheless, within a certain reasonable range of parameters (see also section "stable control domains"), DFC still places the eigenvalues close to the initial position they had before the onset of pathology. Therefore, as we showed above, aspects of both rate and temporal coding that the network may be performing are recovered.

Effects of neuronal and synaptic response function
The understanding of the exact mechanisms by which DFC suppressed SI activity allowed us to precisely investigate how the neuron and synapse response function R and S respectively influence the stability of the closed-loop system. To this end, we again used a mean-field approximation, which explicitly incorporates the expressions for R and S. In general, the neuron response R depends on the specific neuron model as well as on the external input. Here, we did not change the neuron model, but altered the external Gaussian white noise input by using different values for the mean and variance (μ, σ 2 ). We then assessed the stability of the system. It is apparent that for a given pair of coupling and control parameters (J, d) and (K, d c ), respectively, the system becomes unstable as we move in the two dimensional parameter-space spanned by the mean and variance (Fig 6a). For meaningful comparison we used (μ, σ 2 )-combinations that yield constant rates. In the ideal case where M(λ) = S(λ) and d c = d Eq (2) becomes: where G s is the slope (or the static gain) of the 'f-I curve' at the operating point and R n the normalized neuron response (see methods). The critical effective coupling is then given by L cr ðlÞ ¼ ðJ À KÞðlÞ ¼ 1 G s jR n ðlÞ Á SðlÞj Lighter shades of red correspond to higher gain G s . The gain changes significantly even along the constant firing rate lines (blue lines, 10-60 sp/s). (c) The most stable control is achieved if the difference Δ d = (d+2τ r )−(d c +b/2) between effective delays of control and synaptic kernels is minimized. The blue dashed line corresponds to d = d c where the optimal kernel width is b = 4τ r . It predicts correctly the stable regime for the range where atan (ωτ r ) % ωτ r (see Methods). In our simulations the effective control delay is d c,eff = 7 ms, which is very close to the optimal value (star). In panels (a) and (c) the stable and unstable regimes are marked by red and white colors respectively. doi:10.1371/journal.pcbi.1004720.g006

Closed-Loop Control of Spiking Neural Networks
As we move along the constant output firing rate lines both G s and |R n (λ)| increase (S1b and S1b Fig) leading to a decrease of L cr . The changes in G s are significantly larger than those in | R n (λ)|, implying that the static gain is the dominant factor that affects stability. The changes in |S(λ)| are negligible (S1c Fig). This is expected, because the frequency range we are interested in is much smaller than the cut-off frequency of the synaptic filter ω < ω 3db . Thus, when the system operates in a dynamic regime in which single neuron responses have a higher gain the control domains shrink and the range of K values that stabilizes the system decreases.
Next, we investigated the interaction between the synaptic S(λ) and the control kernel M(λ). The amplitude responses for different kernels do not vary significantly (S2 Fig). Therefore, the important factor that influences stability is the phase difference or, alternatively, the difference Δ d between the effective delays of the synaptic d eff and the coupling kernel d c,eff . An optimal result is achieved if this difference vanishes (see methods) i.e. when This point is illustrated for the case where d c = d+1ms (Fig 6c). These results show that DFC does not strongly depend on the shape but rather on the effective delay of the kernel.

DFC induced SI activity
Interestingly, the same control strategy can be used to induce or enhance rather than to suppress oscillations. Choosing appropriate control parameters to increase the effective coupling, i.e. selecting K to have the same sign as J (see methods), results in SI activity (Fig 7). This may be helpful for the treatment of symptoms in several pathological conditions that are characterized by impaired oscillations, e.g. gamma power decrease in schizophrenia [35]. Thus, DFC is a generic control approach and the control parameters (K,d c ,M) can be tuned to quench or to enhance oscillatory activity, depending on the nature of aberrant activity.

Discussion
Open-loop stimulation has been the main non-pharmacological approach to control the symptoms in a wide range of pathological conditions. It has been successful to some degree, but it often introduces clinical side-effects [36]. Moreover, it inherits the drawbacks from open-loop systems: (i) The stimulation profile is predetermined and is not adjusted to the clinically observed short-term fluctuations in the patients' symptoms [37]. In addition, stimulation is continuously applied even though it may not be always necessary. (ii) The stimulation does not adapt to long-term changes of the system, e.g. structural alterations due to the progression of the disease. (iii) The operating point cannot be altered to deal with perturbations, caused, for instance, by a drift of the electrode lead [38]. (iv) External disturbances due to transient undesired signals are not being suppressed. In contrast to open-loop control, closed-loop control can by design deal with all these situations. For this reason there is a growing interest in investigating feedback-control both experimentally [2,[22][23][24] and theoretically [39]. The goal of the experimental work has been to demonstrate that closed-loop control is indeed effective, whereas the theoretical studies aimed at providing a deeper understanding of the underlying conditions and mechanisms.
Here, we provide a theory for DFC, a conceptually simple but powerful form of control [5,6], applied to the control of stochastic SI oscillations in SNN. These oscillations are generic, they occur in many brain areas and in multiple conditions [40] and they emerge via a supercritical Hopf bifurcation [16]. Therefore, the control objective was specific: to counteract this bifurcation. We provide a mean-field approximation to estimate the DFC parameters and confirm the analytical predictions in numerical simulations in purely inhibitory and in coupled excitatory-inhibitory SNNs. DFC is effective in networks with arbitrary connectivity as long as the temporal instabilities, that is uniform oscillations, dominate the dynamics.
We used two control approaches, direct and differential control, and demonstrated that both schemes are effective in suppressing oscillations. Consistent with previous findings [5,30], our results reveal that differential control has two main advantages over direct control. First, the control domains are enlarged, which renders the selection of control parameters an easier task. Larger control domains imply increased robustness of the system both to perturbations in the parameters and to disturbances. This means that neither small deviations from the nominal values of K, d c 1 , d c 2 nor external signals compromise its stability. Second, in differential control the stimulation signal vanishes which translates to decreased power consumption. In clinical settings this is a highly desirable property and is, in fact, a basic requirement of any neuroprosthetic device.
The key advantage of the approach we presented here is that the system under control is being steered back towards its primary operating point (Fig 8). That is, DFC effectively decreases the synaptic coupling strength and, therefore, it counteracts the causes that originally led to the instability. This is obviously true only for the first-order statistics, because DFC does not counteract changes in the variance of the input that a random neuron in the network receives. Nevertheless, this is sufficient for the network to recover basic processing abilities both for rate and temporal coding schemes. Alternative approaches that rely on increased external noise are able to suppress oscillations [41], but they do not allow the network to perform any meaningful computations. We think that a similar explanation is valid also for the traditional open-loop DBS. The exact mechanisms of this type of DBS are still debated [42], but one of the reasons for the induced side-effects may be compromised information processing.
DFC suppresses oscillations in SNN not by decorrelating individual neurons, but rather by applying a common signal to all neurons that counteracts the mean input they receive from the network. Besides the many advantages described earlier, the utility of DFC lies in the fact that Alternative approaches, e.g. noise injection may suppress oscillations, but they drive the system to a dynamical regime characterized by a different transfer function T 2 (ω) where physiological computations are impaired. doi:10.1371/journal.pcbi.1004720.g008

Closed-Loop Control of Spiking Neural Networks
it is a very general control strategy, which does not depend qualitatively on lower level properties such as the specific coupling kernels of the connections. Neither does it depend qualitatively on the exact neuronal type (S3 Fig). Thus, DFC can in principle deal with more complex scenarios, e.g. heterogeneities in network and in neuronal properties. Nonetheless, the exact shapes of the neural and synaptic response functions do affect the system quantitatively and do modify the stability landscape. Thus, it is essential to have a good understanding of their precise contribution.
The results we presented provide a clear picture about how exactly the static gain of the neuron affects the size of the stable control domains. We also showed that the width rather than the shape of the control kernel affects the stability boundary. Further, we explained how the coupling strength and delay influence the overall stability landscape. Besides their theoretical value, these insights have a direct implication for the design of neuroprosthetic devices and are, therefore, of immediate practical and clinical relevance.
It is also important to address certain limitations in our approach. First, we assumed that we can directly affect the neuron membrane potential. This is currently only possible in animal models, where optogenetic tools are applied not only for superthreshold excitation/inhibition but also for subthreshold modulation [32]. In humans only electrical stimulation is realistic at the moment, therefore incorporation of volume conduction models [43] to describe the effects of electrical fields on individual neurons would be required [44]. Second, we used the population average of single neuron firing as our observable. Again, in more realistic settings the population rate has to be reconstructed from the available recordings, e.g. from multi-unit activity. This could be done, for instance, via the use of Kalman filters [45]. Third, our theoretical analysis was based on a mean-field approximation that ignores fluctuations in the input. Analytical and numerical results were largely in a good agreement, but additional work is necessary to specifically deal with the fluctuations in the activity. Last, we had access to the relevant parameters required for the tuning of the controller. This was intended, because our goal was not to simply demonstrate that the method is effective in large networks of spiking neurons, but also to provide an explanation in terms of mean-field theory of the precise mechanism. In real applications these parameters have to be estimated online from the recorded activity. The delay could be inferred from the frequency of the oscillatory activity. Inference of the coupling strength is less straightforward, but may still be feasible. Alternatively, once the delay is estimated, methods of adaptive tuning could be used to retrieve also the optimal control gain (S4 Fig). Tuning the controller is in general a difficult problem, even for open-loop DBS, and additional research in this direction is required.

Differences from previous work
Few studies have addressed the problem of suppressing oscillations in neural activity (see [39] for a detailed review). They are based (i) on population dynamics [46] [15], (ii) on detailed single neuron descriptions [13], (iii) on simplified but computationally efficient models (e.g. Rulkov maps [47]), (iv) or on combinations thereof [14]. These approaches have their merits, but they come with limitations: (i) the parameters cannot be directly mapped to experimental measurable quantities (ii) it is not clear if the results scale to large networks of neurons.
The approach that we presented here is a trade-off between biophysical realism and analytical tractability. We used the LIF model, which captures single-neuron dynamics to a sufficient degree, while at the same time allows computationally efficient simulations of large networks. We applied DFC that was originally proposed in the context of chaotic systems as a method to control unstable periodic orbits [5]. DFC has been also used to control dynamics in extended media [48] and has been applied in different contexts as well [49]. It was later used to control coherence [50] and to suppress synchronous activity in networks in which the neurons themselves act as oscillators [10][11][12].
Here, we did not use simplified population dynamics or phase oscillators. Instead, we used spiking neurons that fire irregularly and are nevertheless able to generate oscillations. We used a single proportional control term to provide a proof of principle of the method and to be able to delineate the control domains semi-analytically. Alternative approaches, e.g. linear PI/PID control [4] or non-linear control schemes [11] are also possible and may improve performance, however, their theoretical analysis is less straightforward. We also inserted realistic descriptions of synaptic dynamics and, therefore, were able to explicitly study their contribution to stability. This allowed us to design an appropriate control kernel, which resulted in increased control domains. In addition, by using a mean-field theory that explicitly incorporates the synaptic and neuronal response functions we could study their contribution in a systematic way. The neuronal response function enabled us to investigate the influence of external and recurrent inputs and to relate them to experimentally measurable quantities. Indeed, as we showed above, the statistics of the mean field for activity states with very similar firing rate profiles may be significantly different affecting stability. Therefore, feasible measurements of the population activity can be directly used to characterize the operating point of the network and to fine-tune the control parameters to achieve the desired results.
Finally, the results presented here provide us with an understanding of the recent success of event-triggered control strategies. Event-triggered control can be placed between open-loop and continuous closed-loop control (e.g. DFC). Open-loop provides constant magnitude stimulation independent of the ongoing activity. Event-triggered approaches provide also constant magnitude stimulation, but only if a certain event occurs, for instance the power of beta oscillations crosses a certain threshold. Thus, the overall stimulation time and, therefore, the undesired stimulation side-effects are reduced. In DFC the stimulation side-effects are likely to be further reduced, because the stimulation amplitude is continuously adjusted to the ongoing activity. That is, no excessive stimulation is applied. This is particularly true for differential DFC, in which the stimulation amplitude vanishes over time.

Conclusions
We used DFC, a relatively simple form of control that includes only a proportional gain term, because it is still possible to analytically study the stability of the closed-loop control system. More sophisticated control strategies could further increase the performance of the system. They come, however, at the price of increasing the number of control parameters that have to be estimated and of increasing complexity precluding a formal proof of stability. Our approach spans multiple levels of analysis of neuronal dynamics, enabling an understanding of how the control stimulus interacts with both low-level synaptic and high-level properties of the population activity to influence stability. At the same time the complexity of the controller is low enough to be of practical relevance. Thus, here we have provided a general conceptual framework for future studies that address both theoretical and practical aspects of closed-loop control in neuronal systems.

Numerical simulations
We simulate networks of N LIF neurons randomly connected with a probability of = 0.1. Thus each neuron receives on average C = N connections from other neurons in the network. For the purely inhibitory network we use N = N I and for the coupled excitatory-inhibitory case N = N E +N I . The subthreshold dynamics of a neuron i in the network is given by where R m is membrane resistance, τ m is the membrane time constant and v rest is the resting potential. The recurrent input term describes the total synaptic current arriving at the soma due to presynaptic spikes. c ij are elements of the binary connectivity matrix. Each presynaptic spike causes a stereotypical postsynaptic current s(t) modeled as an α-function [51] sðtÞ where τ s is the synaptic time constant and H(t) the Heaviside function. The double sum in Eq 6 runs over all firing times t k j of all presynaptic neurons connected to neuron i. For all connections in the network we use the same synaptic coupling strength J ij = J/ C, where C is the average in-degree and d ij = d the transmission delay. The external input contains a mean term μ and a fluctuating term resulting from the Gaussian white noise η i (t) that is uncorrelated from neuron to neuron with <η i (t)> = 0 and <η i (t)η i (t 0 )> = δ(t−t 0 ).

Asynchronous state
In the stable asynchronous state the population average of the firing rate is constant in time, r(t) = r 0 . The mean recurrent input that each neuron receives is therefore also constant and given by similarly the variance of the recurrent input is We study the stability of the asynchronous state following a linear perturbation approach [18]. A small oscillatory modulation of the stationary firing rate r(t) = r 0 +r 1 e −λt with v 1 ( 1 and λ = x+jω where ω is the modulation frequency leads to corresponding oscillation of the synaptic current The firing rate in response to an oscillatory input is given by ðy t ; lÞ À @U @y ðy r ; lÞ Uðy t ; lÞ À Uðy r ; lÞ The function U is given in terms of combinations of hypergeometric functions In a recurrent network the modulation of the firing rate and the modulation of the synaptic input must be consistent. Combining Eqs (9) and (10)  SðlÞ ¼ e Á t s ð1 þ l Á t s Þ 2 describe the neuronal and synaptic response functions respectively. The negative sign of J is absorbed in the phase of S(λ).
The critical coupling values at which modes have marginal stability with frequency ω i can then simply be computed by Stability. The solutions λ of the self-consistency Eq (11) determine the stability of the system. If for all solutions the real part is negative, Re{λ} <0, then the system is stable otherwise it is unstable. The stability border, λ = jω, is characterized by the occurrence of a supercritical Hopf bifurcation. At this point the population activity will be oscillatory with frequency ω.

Delayed feedback control
In the simulations we implement DFC by recording and stimulating all neurons in the network. The subthreshold dynamics of a neuron i with DFC is given by where I C (t) is the control input. Note that I C (t) is identical for all neurons in the network given by where v(t) is the instantaneous population activity at time t and ? denotes the convolution operation ðf ? gÞðtÞ ¼ R 1 À1 f ðt À tÞ gðtÞ dt: We used as control kernel m(t) a box function mðtÞ ¼ Hðt À aÞ À Hðt À bÞ where H(t) is the Heaviside function Thus the control input I C (t) was updated in steps of 1ms. Direct control. In the case of direct DFC a modification of the self-consistency Eq (11) yields where M(λ) describes the control kernel in the frequency domain and the negative sign captures the fact that the control stimulus counteracts the effects of the synaptic coupling J.
The box function has response characteristics given by where a = 0, b is the width of the kernel and d c is the control delay. Differential control. In differential DFC where the control input I C (t) is a function of the difference between two time-delayed versions of the population activity v(t). It is given by In this case the self-consistency Eq (11) is modified to give Rate compensation. In all simulations we adjust the mean μ and variance σ of the external input to the neurons such that the firing rates are approximately equal for all conditions, that is where v 0 , v J , v K are the firing rates of the uncoupled, coupled and network under DFC respectively. In this way we can exclude any effects due to changes in the firing rates.
Stability analysis. The eigenvalues λ of the self-consistency Eqs (14) and (16) determine the stability of the system. We compute for both direct and differential control the real part of the rightmost eigenvalue Re{λ 1 } that determines stability. We use the (I K , d c )-parameter pair with I K 2 [0, 300] mV and d c 2 [0, 30] ms. The second delay term in differential control was in both cases d c2 = 1 ms.
DFC induced SI activity. If the control gain K has the same sign as the synaptic coupling J and the control delay is chosen to be close to the synaptic delay, d c ' d, then the effective coupling in the network increases resulting in SI activity (Fig 7). In this case the self-consistency equation is given by Static gain. In the frequency domain the neuron response function R(f) is simply the Fourier Transform of the impulse response h(t), i.e. Rðf Þ ¼ F ½h 0 ðtÞðf Þ. The impulse response h(t) can be separated in two parts where h n (t) is the normalized impulse response such that R 1 0 h n ðtÞdt ¼ Rðf ¼ 0Þ ¼ 1 and G s = @ μ r is a constant term that we denote as the static gain of the response. It corresponds to the slope of the 'f-I' curve at the operating point and captures the 'susceptibility' of the rate to small changes in the mean μ.
We can rewrite the self-consistency Eq 4 as ðJ À KÞ Á RðlÞ Á SðlÞ Á e Àld ¼ 1 ðJ À KÞ Á G s Á R n ðlÞ Á SðlÞ Á e Àld ¼ 1 ðJ À KÞ Á G s Á jR n ðlÞ Á SðlÞje i ¼ 1 where we have separated the complex expression R n (λ) Á S(λ) in an amplitude and phase part. Splitting this complex equation in two real ones we get and ðJ À KÞ Á G s jR n ðlÞ Á SðlÞj ¼ 1 ) L cr ðlÞ ¼ ðJ À KÞðlÞ ¼ 1 G s jR n ðlÞ Á SðlÞj The first equation gives us the critical frequencies ω cr for which the modes exhibit marginal stability, i.e. λ = iω cr . From the second equation we can compute the corresponding effective critical coupling L cr (ω cr ).
Control kernel dependence. To study how stability depends on the control kernel we write again the self-consistency equation ðJ þ DJÞ Á RðlÞ Á SðlÞ Á e Àld À K Á RðlÞ Á MðlÞ Á e Àld c ¼ À1 ððJ þ DJÞ Á SðlÞ Á e Àld À K Á MðlÞ Á e Àld c Þ Á RðlÞ ¼ À1 where ΔJ is a pathological increase in the synaptic coupling that the controller needs to counteract. Separating amplitude and phase we get ððJ þ DJÞ Á A s ðlÞe i S À K Á A M ðlÞe i M Þ Á RðlÞ ¼ À1 J Á e i S þ ðDJ Á e i S À K Á e i M Þ Á RðlÞ ¼ À1 where we used A S (iω)) % A M (iω) = 1, which is valid for the frequency range that we are interested in ω < 300 rad (S2 Fig). To counteract the increase ΔJ we need to minimize the effective coupling term For the control kernel we use a box function defined as whereas the synaptic kernel is the α-function SðioÞ ¼ e Á t r ð1 þ iot r Þ 2 e ÀioÁd and if we assume that atan(ω Á τ r )%ω Á τ r then Thus, the most optimal control is achieved when the effective delays of synaptic and control kernel d s,eff and d c,eff respectively are identical D d ¼ d s;eff À d c;eff ¼ 0 In our simulations we used d = 5 ms and τ r = 1ms, thus d s,eff = 7ms. For the controller we used d c = 6.5 ms and b = 1ms, thus d c,eff = 7ms.

Recovery of network function
Response to random pulse packets. We computed the response of the network to incoming stimuli that arrived in form of random Gaussian pulse packets. A pulse packet was composed of a predefined number of spikes n pp = 100 with normally distributed random displacements t k * N(μ = 0, σ pp = 10ms) from the center time t c of the pulse. It was fully defined by the tuple (t c ,n pp , σ pp ). In total, ten pulse-packets with center times t k = (200k+200) ms and 1 k 10 were applied to 100 randomly chosen neurons in an inhibitory network of size N = 1000. We computed the population response at time points t k + 10 ms and compared it with the population activity during baseline at time points t k + 100 ms. For this, we performed a receiver operating characteristic (ROC) analysis evaluating the true positive and false positive rate for various thresholds. We then computed the area under the ROC curve (AUC), which indicates how well the response can be distinguished from baseline activity. An AUC value of 1 means full separability of the two activity states, whereas an AUC value of 0.5 indicates full overlap of the activity sampled during the two different conditions. We computed the AUC values for three different scenarios: (i) physiological AI state (ii) pathological (oscillatory) state controlled with DFC (iii) pathological (oscillatory) state with noise. The results are shown in Fig 5a and 5b.
Response to common input. We defined a spike train of n ST = 500 equally spaced spikes in a window of T ST = 50 ms. Ten copies of exactly the same spike train with time onset t k = (200k + 200) was provided as input to 100 randomly chosen neurons in an inhibitory network of size N = 1000. Thus, in this scenario all stimulated neuron received identical input during the stimulation periods. However, in this case we were interested in the temporal aspects of the network response. To this end, we measured the synchrony between the spike trains of all neurons in the network using the SPIKE-distance metric [34]. The SPIKE-distance is a measure of (dis)-similarity which allows for a time-resolved analysis and can track instantaneous changes. We computed the multivariate SPIKE-distance S both during the 50ms of stimulation (S ST ) and also during 50 ms of baseline activity (S BL ). We then computed the temporal average for the stimulation SðtÞ Á dt and for the baseline D BL (t) = D ST (t+100). The results, again for the three scenarios described above (AI, DFC, noise) can be seen in Fig 5c. Oscillation index. We estimated the discrete power spectral density P(ω) of the population activity r(t) using the standard Fast Fourier Transform (FFT) method. We then computed the total power in the range [0, 250/π] rad P T ¼ P i Pðo i Þ and used log 10 P T as a descriptor of oscillation strength.