Controlling Synchronization of Spiking Neuronal Networks by Harnessing Synaptic Plasticity

Disrupting the pathological synchronous firing patterns of neurons with high frequency stimulation is a common treatment for Parkinsonian symptoms and epileptic seizures when pharmaceutical drugs fail. In this paper, our goal is to design a desynchronization strategy for large networks of spiking neurons such that the neuronal activity of the network remains in the desynchronized regime for a long period of time after the removal of the stimulation. We develop a novel “Forced Temporal-Spike Time Stimulation (FTSTS)” strategy that harnesses the spike-timing dependent plasticity to control the synchronization of neural activity in the network by forcing the neurons in the network to artificially fire in a specific temporal pattern. Our strategy modulates the synaptic strengths of selective synapses to achieve a desired synchrony of neural activity in the network. Our simulation results show that the FTSTS strategy can effectively synchronize or desynchronize neural activity in large spiking neuron networks and keep them in the desired state for a long period of time after the removal of the external stimulation. Using simulations, we demonstrate the robustness of our strategy in desynchronizing neural activity of networks against uncertainties in the designed stimulation pulses and network parameters. Additionally, we show in simulation, how our strategy could be incorporated within the existing desynchronization strategies to improve their overall efficacy in desynchronizing large networks. Our proposed strategy provides complete control over the synchronization of neurons in large networks and can be used to either synchronize or desynchronize neural activity based on specific applications. Moreover, it can be incorporated within other desynchronization strategies to improve the efficacy of existing therapies for numerous neurological and psychiatric disorders associated with pathological synchronization.


INTRODUCTION
Most of the existing therapies for neurological disorders focus on suppressing symptoms and ignore the dynamical aspects of the underlying network that create the pathological symptoms. Therapies that purely focus on suppressing the pathological symptoms instead of addressing aberrant dynamic characteristics of the network could fail to achieve a long term behavioral recovery (Popovych and Tass, 2014). One such example is the treatment of Parkinson's Disease motor symptoms by suppressing the strong oscillatory firing pattern of neurons in the brain using high frequency stimulation (HFS) (Temperli et al., 2003;Hammond et al., 2007;Benabid et al., 2009;Postuma et al., 2015;Singh, 2018). The applied HFS suppresses the oscillatory activity and the Parkinsonian motor symptoms while the stimulation is on but when the stimulation is off, the oscillations and motor symptoms reemerge (Temperli et al., 2003;Deuschl et al., 2006). Additionally, epilepsy is also characterized by an excessive synchronous firing pattern of neurons, which causes epileptic seizures (Duncan et al., 2006). Similar to Parkinson's diseases, HFS helps to control the excessive synchronization of neurons and epileptic seizures (Klinger and Mittal, 2016). To improve upon open-loop stimulation protocol, new neural stimulation strategies have been developed to more effectively and efficiently desynchronize pathologically synchronous neuronal networks throughout the brain.
Recently, control-theoretic approaches have been used to design more effective and energy efficient desynchronization strategies (Popovych and Tass, 2014). These optimal stimulation protocols include a single pulse minimum energy desynchronizing control input (Deuschl et al., 2006;Nabi et al., 2013a,b;Mauroy et al., 2014;Wilson and Moehlis, 2014;Monga et al., 2018) and closed-loop delayed feedback Popovych et al., 2005Popovych et al., , 2017Kiss et al., 2007;Vlachos et al., 2016) approaches. The minimum energy control input approach designs a pulse that pushes the state of the network to a phaseless set-point, which is the point where all the isochrons of the system converge (Nabi et al., 2013a,b). Here, the network's inherent noise randomly pushes each individual neuron onto a different isochron with its own phase. The net effect of this random reset is an asynchronous population activity. Another approach is closed-loop delayed feedback control. In the closed-loop delayed feedback control approach (Vlachos et al., 2016), the time-delayed average population activity is used as a feedback to desynchronize the network. Since this approach only feeds the past population activity, the applied desynchronizing input is only active whenever the network becomes synchronous. While these approaches provide an optimal desynchronization strategy, most of them assume that the network connections are static and ignore the inherent plastic nature of neuronal synapses (Abbott and Nelson, 2000).
Hebbian plasticity is a well-known form of activity-dependent synaptic plasticity (Abbott and Nelson, 2000). This form of plasticity enforces productive connections between neurons that produce action potentials and depresses unproductive connections that do not elicit action potentials (Hebb, 1949). One form of activity-dependent synaptic plasticity is spiketime dependent plasticity (STDP) (Song et al., 2000). This rule increases the weight of a synaptic connection when the presynaptic neuron fires before the post-synaptic neuron within a given time window and decreases the weight when the order is reversed (Song et al., 2000). An increase or decrease in the synaptic weight is coined long-term-potentiation (LTP) or long-term-depression (LTD), respectively. The introduction of plasticity into a neuronal network creates multiple stability points with different levels of synchronous activity (Tass and Hauptmann, 2007;Pfister and Tass, 2010;Popovych and Tass, 2014). Since the connections are plastic, an external stimulus can move the network from one stability point to another in order to drive the network from a synchronous to an asynchronous state.
To incorporate the synaptic plasticity within the desynchronization stimulation protocol, Coordinate Reset (CR) based stimulation strategies have been developed, which temporarily makes the pathologically synchronous stability point unstable to move the network into the asynchronous state by harnessing synaptic plasticity (Tass, 2003a,b;Tass and Majtanik, 2006;Tass and Hauptmann, 2007;Pfister and Tass, 2010;Ebert et al., 2014;Zeitler and Tass, 2015). In this approach, the network is driven to the asynchronous regime in the presence of external stimulation because of the inherent domination of LTD in asynchronous bistable networks. Once CR stimulation moves the network from the synchronous to the asynchronous stability point, the stimulation input is no longer required (Tass and Hauptmann, 2007;Pfister and Tass, 2010). This results in a long-lasting desynchronization of the network after the stimulation protocol is turned off (Tass and Majtanik, 2006). While this protocol is effective, it has its own limitations. CR based stimulation strategies rely on acute desynchronization of the networks regardless of the plasticity rules during the stimulation. If the network is not bistable, the acute-desynchronization of the network may not depress the synaptic weight of the network and the desynchronization of the network will be brief (Pfister and Tass, 2010). This lack of bistability may arise in a network where LTP dominates LTD (i.e., pathological networks). Moreover, this strategy only provides a way to desynchronize neural activity and not to control the synchronization level of the network.
In this work, we have developed a novel stimulation strategy "Forced Temporal Spike-Time Stimulation (FTSTS)" which addresses above shortcomings of the CR-based stimulation approach. While all other stimulation strategies focus on desynchronizing neural activity within a network, our strategy focuses on harnessing the underlying synaptic plasticity of the network to control the average network synaptic strength by forcing the spiking neurons to fire in specific temporal patterns. Thus, our strategy provides complete control over the synchrony level of networks for a long period of time, not just desynchronization. We demonstrate the efficacy of FTSTS strategy in controlling the desired synchrony level in large excitatory-inhibitory (E-I) networks. We show in simulation that the FTSTS strategy can effectively desynchronize the neural activity in networks where LTP dominates LTD on average. Further, we combine the FTSTS strategy with the CR stimulation strategy to demonstrate how this can enhance the overall performance of the CR stimulation strategy in desynchronizing large networks.
The paper begins with a description of models used to describe the spiking E-I networks dynamics, STDP rules, and a measure of synchrony as well as the stability analysis of E-I networks in section 2. In section 3, we first provide a mechanistic understanding of the FTSTS strategy by considering control of synchrony in a two neuron E-I network. We then demonstrate the efficacy of the FTSTS strategy in desynchronizing large E-I networks subject to different plasticity rules, and uncertainties in the network parameters and the designed stimulation parameters. Finally, we show how the FTSTS strategy can be incorporated within the CR stimulation strategy to improve the overall performance of the CR stimulation strategy. The paper ends with a detailed discussion on the comparison of our approach with the existing stimulation strategies for desynchronization of spiking neural networks as well as the limitations of our strategy in section 4.
Here, v E (t) and v I (t), in millivolts (mV), represent the membrane potential of the excitatory and inhibitory neurons, respectively. τ m (in ms) is the membrane time constant. Z i (t) denotes the synaptic input to the i th population of neurons where i ∈ {E, I}.
The synapses between the excitatory and inhibitory populations are randomly connected with a probability of ǫ. The synaptic input function: defines the input to the i th neuron population. In Equation (3), J ij represents the synaptic strength between a presynaptic neuron in population j and postsynaptic neuron in population i, in mV, where i ∈ {E, I} and j ∈ {E, I}. For example, the synaptic strength of a I-to-E synapse is J EI . C ij = 0.3N tot denotes a scaling factor where N tot is the total number of neurons in the network. S ij (t) is the synaptic function. The Gaussian distributed baseline current to the i th type neuron is denoted as: with a mean baseline current of µ i and a variance of σ 2 i τ m . χ(t) is white noise with a mean of 0 and a variance of 1. Finally, V i stim (t) denotes the external stimulation input to the i th neuron population.
The synaptic function S ij (t) is modeled as (Brunel and Hansel, 2006): Here, X ij describes the input to the i th population of neurons from the j th population of neurons. The time constants governing the decay and rise time are τ d (in ms) and τ r (in ms), respectively. Synaptic connections between the i th and j th neuron populations are randomly connected with a probability of ǫ. The weight of each synaptic connection is defined as W ij . Throughout the work, we assume that E-to-I connections (W IE (t)) are plastic and the I-to-E connections (W EI ) are static except in section 3.9 and Figure 11 where we consider both connections to be plastic. Unless otherwise specified, we further assume no synaptic connectivity among neurons in excitatory or inhibitory populations. The Dirac-Delta function δ(t − t pre + t delay ) models the synaptic input to a postsynaptic neuron from a presynaptic neuron when the presynaptic neuron fires at time t pre (in ms) with a synaptic delay of t delay (in ms).

Spike-Timing Dependent Plasticity (STDP) Model
The coupling value of the plastic E-to-I synapse (W IE (t)) is governed by STDP (Song et al., 2000), which is defined as follows: where W IE (t) is given as: Here, W IE (t) defines the change in the synaptic weight determined by the spike-time of a presynaptic (t pre ) and postsynaptic (t post ) neuron. The rate at which the E-to-I synaptic coupling changes is governed by the learning rate η e . Additionally, the relative contribution of LTD and LTP to W IE (t) is denoted by a LTP and a LTD . (Song et al., 2000;Ebert et al., 2014). Since LTD generally dominates LTP, a LTD is 10% larger than a LTP . A post (t) and A pre (t) are described by the following two exponential functions: Here, the size of the LTP and LTD time window is defined by the STDP time constants τ LTP and τ LTD . In a similar fashion to a LTP and a LTD , τ LTD is set to be 10% greater than τ LTP . Upon the firing of a presynaptic or postsynaptic neuron, a small value A 0 is added to the appropriate exponential STDP decay function. The E-to-I synaptic weight is defined as the coupling value W IE (t) multiplied by the synaptic strength J IE (i.e., J IE W IE (t)).

Synchrony Measurement
We measure the synchrony level of the network by computing the Kuramoto order parameter (R(t)) based on the spike times of neurons in the excitatory population (Kuramoto, 1984;Daido, 1992;Tass, 2007;Ebert et al., 2014). The phase (φ k (t)) of the k th neuron in the excitatory population is calculated using the following equation: where t k,i is the i th spike time for the k th neuron. The Kuramoto order parameter (R(t)) and the average phase of the neurons (ψ(t)) are calculated using: Here, N E represents the number of excitatory neurons. A highly synchronous network has a Kuramoto order parameter of R(t) = 1 and a completely asynchronous network has a value of R(t) = 0.

Determination of Synchronous and Asynchronous Regimes
It is well-known that plastic neural networks exhibit multiple stability points (Song et al., 2000;Tass and Hauptmann, 2007;Pfister and Tass, 2010;Popovych and Tass, 2014). Similar to other networks, our E-I network exhibits two stability points at a high and low average E-to-I synaptic weight value. Figure 1A shows the average E-to-I synaptic weight converging to either J IE W IE = 10 mV or J IE W IE = 290 mV. The average synaptic weight converges to J IE W IE = 290 mV if the initial average synaptic weight is >100 mV and converges to J IE W IE = 10 mV when it <100 mV. Additionally, we find that the network becomes more synchronous as the average E-to-I synaptic weight increases, which is shown in Figure 1B in terms of the Kuramoto order parameter R(t). Therefore, the network exhibits a high level of synchrony at high synaptic weights and a low level of synchrony at low synaptic weights.

Model Parameters
We use the model parameters defined in Table 1 unless stated otherwise. All the simulation are performed in Matlab R2016b. The differential equations are solved using Euler's method with a step size of 0.1 ms. All our codes and relevant simulation parameters to replicate the presented results in this manuscript are available on our research webpage (https://webpages.uidaho. edu/gkumar) and in the Supplementary Material.

RESULTS
We begin by providing an insight into the underlying mechanism of our stimulation strategy "Forced Temporal Spike-Time Stimulation" (FTSTS) using an illustrative example of a two FIGURE 2 | The FTSTS protocol for a two neuron E-I network. (A) Shows an excitatory-inhibitory network. (B,C) Show the FTSTS input pattern for the inhibitory and excitatory neuron, respectively. The FTSTS pulse parameters are U stim = 300 mV, T stim = 1 ms, and T neutral = 10 ms. This FTSTS protocol depresses the E-to-I synaptic weight as shown in (D). The drop in synaptic weight is due to FTSTS inducing a post-before-pre spiking pattern in the E-I network, which is shown in (E).
On the other hand, (F) Shows how swapping the FTSTS inputs to the excitatory and inhibitory neuron increases the E-to-I synaptic weight. This induces the pre-before-post spiking pattern shown in (G).
neuron E-I network in section 3.1. Next, we demonstrate the efficacy of the FTSTS strategy in controlling synchronized activity of a 2,000 and 10,000 neuron E-I networks (see sections 3.2, 3.3). We then show the robustness of the FTSTS strategy in the presence of uncertainties in the designed stimulation pulses, model parameters and network connectivity (see sections 3.4, 3.5, 3.6, 3.7). Next, in section 3.8, we combine the FTSTS strategy with the existing coordinate reset (CR) stimulation strategy to show the efficacy of the FTSTS-CR strategy over the CR stimulation strategy. Finally, in sections 3.9 and 3.10, we demonstrate that FTSTS strategy can desynchronize E-I networks with additional plastic synapses and symmetric spiketime plasticity rules.

Control of E-to-I Synaptic Weight in a Two Neuron Network
We considered an excitatory-inhibitory (E-I) network of two neurons to develop our "Forced Temporal Spike-Time Stimulation" (FTSTS) strategy (Figure 2A). We set a scaling factor of C ij = 10N i and probability of connectivity of ǫ = 1 (see section 2.1 for the meaning of these variables). Based on the STDP rule of activity-dependent plasticity, we designed stimulation inputs for both the inhibitory (V I stim (t)) and excitatory (V E stim (t)) neuron that forced the postsynaptic inhibitory neuron to spike before the presynaptic excitatory neuron, as shown in Figures 2B,C, respectively. The protocol stimulated the postsynaptic inhibitory neuron and the presynaptic excitatory neuron using charge-balanced rectangular pulses with an equal and opposite amplitude (U stim ). Figure 2E shows the induced firing patterns in neurons, which decreased the average E-to-I synaptic weight of the network as shown in Figure 2D. On the other hand, the average E-to-I synaptic weight increased when V E stim (t) and V I stim (t) were switched such that the presynaptic excitatory neuron fired before the postsynaptic inhibitory neuron. The increased synaptic weight observed from the induced spiking pattern in Figure 2G is shown in Figure 2F. (E-G) Show the spiking patterns before, during, and after the FTSTS protocol, respectively. The FTSTS pulse parameters are U stim = 100 mV, T stim = 1 ms, and T neutral = 10 ms.

FTSTS Effectively Controls the Neuronal Synchronization in a 2000 Neuron Network
We applied the FTSTS strategy to a E-I network of 1,600 excitatory and 400 inhibitory neurons to demonstrate how our strategy can be used to control the synchrony of neuronal activity in large networks. In a larger network of neurons, our strategy forces the postsynaptic inhibitory population of neurons to spike before the presynaptic excitatory neuron population. We assumed that all the neurons in each specific population receives the same input. The applied FTSTS inputs to each neuron population are shown in Figures 3A,B. These inputs induced a specific spiking pattern, as shown in Figure 3F, which depressed the average E-to-I synaptic weight (shown in Figure 3C). The period of stimulation is highlighted with a solid black line in Figure 3C. Since the network has an asynchronous regime that converges to a low average E-to-I synaptic weight, we only required enough input to drive the network into the asynchronous regime. Therefore, we provided enough input to depress the synaptic weight of the network to 75 mV, which is slightly over the synchronous-asynchronous regime boundary (see Figure 1). As a result, the system naturally converged to the low synaptic weight stability point when the FTSTS protocol was turned off. The synchronous and asynchronous firing patterns before and after the stimulation protocol are displayed in Figures 3E,G, respectively. The synchrony level of the network as it transitioned from the synchronous to asynchronous regime is shown in Figure 3D. The activity of the E-I network prior to the applied stimulation (the first 2 s) was measured around R(t) = 0.75. When the stimulation protocol was turned on, the measured synchrony level became low (see Figure 3D). This is due to the asynchronous firing between each FTSTS pulse. When the FTSTS protocol was turned off, the network remained in the asynchronous regime at the measured network synchrony level of R(t) = 0.05 (see Figure 3D).
Next, we demonstrate how our FTSTS strategy can also be used to synchronize the asynchronous network activity. To do so, we swapped the inputs to the inhibitory and excitatory neurons used in the desynchronization case, which are shown in Figures 4A,B. This stimulation protocol forced the presynaptic excitatory neuron population to fire prior to the postsynaptic inhibitory neuron populations, which is shown in Figure 4F. Similar to the two neuron case, swapping the inputs to the excitatory and inhibitory neurons induced LTP in the network and increased the average E-to-I synaptic weight of the network, as shown in Figure 4C. Again we are only required to drive the network into the synchronous regime (i.e., the average E-I synaptic weight above 100 mV) to synchronize the network. Therefore, the stimulation protocol was turned off when the average E-to-I synaptic weight reached 125 mV, which was slightly over the threshold. Here, the network will remain in the synchronous regime and the average E-to-I synaptic weight will converge to the high synaptic weight stability point. The spiking patterns of the E-I network before, during, and after the FTSTS protocol are shown in Figures 4E-G, respectively. Figure 4D shows the changes in the network synchrony level before, during, and after the FTSTS protocol. As shown in Figure 4D, the network synchrony increased from R(t) = 0.05 to R(t) = 0.6 after the removal of the stimulation. It should be noted that the synchrony level increases significantly during the stimulation in this case compared to the case of desynchronization (see Figure 3D for comparison). This is due to the large input of U stim = 200 mV compared to the input used in Figure 3, which forces the neurons to fire in a highly synchronous firing pattern. If a smaller input is used to resynchronize the network, it would be more noisy during the FTSTS protocol and the FTSTS would be required for a longer period of time to resynchronize the network.

Desynchronization of Neural Activity in a Large E-I Network
We applied our FTSTS protocol to demonstrate its applicability in larger networks. For demonstration purpose, we considered a E-I network with 8,000 excitatory and 2,000 inhibitory neurons. We set the probability of connectivity of the E-to-I and I-to-E synapses ǫ to 0.01. The FTSTS protocol induced the same postbefore-pre firing patterns in the larger network which decreased the average E-to-I synaptic weight, as shown in Figure 5A. The stimulation protocol desynchronized the network in ∼22 s, which is comparable to the desynchronization time of the 2,000 neuron network. The changes in the network synchrony level before, during, and after the FTSTS protocol are shown in Figure 5B. As noted in Figure 5B, the initial synchrony level of R(t) = 0.8 is reduced to approximately R(t) = 0.05 after the FTSTS protocol. Once the stimulation protocol reduced the average E-to-I synaptic weight below 75 mV (i.e., the asynchronous regime), we no longer required the external inputs to keep the network in the asynchronous regime.

Robustness to Uncertainties in the FTSTS Pulse Parameters
Here, we demonstrate the robustness of our protocol in desynchronizing a 2,000 neuron E-I network against uncertainties in the FTSTS pulse parameters. In particular, we considered uncertainty in the FTSTS pulse amplitude U stim , which we modeled in the form of a Gaussian distribution with mean U stim and variance U stim 10 . Each of the applied pulse amplitude during stimulation was randomly chosen from this distribution. As shown in Figure 6B, the FTSTS strategy efficiently desynchronized the network by driving the network into the asynchronous regime. Figure 6A shows the changes in the average synaptic weight of the network before, during, and after the FTSTS protocol.

Robustness to Uncertainties in the Network Model Parameters
In this section, we show the robustness of our FTSTS strategy against uncertainties in the network model parameters. For demonstration, we considered variations in the membrane time constant τ m of neurons in the network. We randomly assigned the membrane time constant τ m of individual neurons in the 2,000 neuron E-I network from a uniform distribution U(8, 12) to show the efficacy of our FTSTS strategy in desynchronizing the network activity. Figure 7A shows our simulation results for τ m ∈ U(8, 12). Here, the FTSTS protocol forced the average E-to-I synaptic weight of the network into the asynchronous regime within ∼15 s of stimulation which led to desynchronization of the network activity after the removal of the stimulation, as shown in Figure 7B. As noticed here, the stimulation protocol desynchronized the network faster in this case compared to the case in Figure 3C where there is no variation in the membrane time constant. This is not surprising as an increase variability in the membrane time constant would induce more noise and desynchronize the firing pattern of the neurons initially, which is seen in Figure 7B with a lower initial Kuramoto order parameter value.   U (8, 12). The applied FTSTS pulse parameters are U stim = 100 mV, T stim = 1 ms, and T neutral = 10 ms.

Addition of E-to-E and I-to-I Synaptic Connections
In this section, we show the efficacy of the FTSTS strategy in desynchronizing 2000 neuron E-I network in the presence of E-to-E and I-to-I synaptic connectivity. We assumed that the synaptic strength of all synapses within the network are static except E-to-I synapses. We set the synaptic strength of the static I-to-E, E-to-E and I-to-I synapses as J EI = 90 mV, J EE = 50 mV, and J II = 50 mV, respectively, with scaling factors of C EE = N tot and C II = N tot , where N tot = N E + N I . The probability of connectivity of the E-to-E and I-to-I was 0.1. The addition of E-to-E and I-to-I synapses within the E-I network didn't change the bifurcation of the regime into synchronous and asynchronous with respect to the network average E-I synaptic weight qualitatively. Our simulation results show that the FTSTS strategy effectively desynchronized the network activity, shown in Figure 8B, in the presence of E-to-E and I-to-I synapses by driving the average E-to-I synaptic weight of the network into the asynchronous stability regime, as shown in Figure 8A.

Robustness to Partial Spatially Inseparable Excitatory and Inhibitory Neuron Population
In this section, we show the robustness of the FTSTS strategy for a case where E-I populations are not well separated. We assume that 25% of the excitatory and 25% of the inhibitory population are not spatially separable. Therefore, this inseparable population of E-I neurons receives inputs designed for the excitatory and the inhibitory population. Figure 9A shows that the FTSTS strategy is still able to push the average synaptic weight of the network into the asynchronous regime. The change in slope during the FTSTS protocol in Figure 9A is most likely due to the synaptic weight of the separable populations reaching the minimum weight value. At this point, the disruption from the stimulation protocol and the low average synaptic weight of the network helped to further depress the synaptic weight of the neurons. Figure 9B shows the reduction of the network synchrony with the Kuramoto order parameter dropping from R(t) = 0.72 to R(t) = 0.05.

Integration of FTSTS With the Coordinate Reset Strategy
In this section, we demonstrate how our FTSTS strategy could be incorporated within the standard coordinate reset (CR) stimulation protocol to effectively stimulate a large population of neurons. One way to implement the CR stimulation protocol is to divide the synchronous population of neurons into four subpopulations, which receive separate but identical inputs at different times over the course of period T (Tass, 2003a,b). T is the overall period of the synchronous neuron population without input. If the neuron population is divided into four subpopulations, then each subpopulations approximately will receive input every T/4. The order that each subpopulation receives input is randomly assigned at every period. In order to compare the CR to the FTSTS strategy, we individually divide the excitatory and inhibitory population into four subpopulations (8 subpopulations for the E-I network). Figures 10E,G show the CR stimulation pattern applied for one period to the excitatory and inhibitory population, respectively. These figures show the stimuli provided to the first quarter of neurons in the E-I network with one pulse to excitatory subpopulation 1 (Figure 10E, blue) and one pulse to the inhibitory subpopulation 2 (Figure 10G, orange). The stimuli to the second quarter of the E-I network is a pulse delivered to the excitatory subpopulation 3 (Figure 10E, yellow) and a pulse to the inhibitory subpopulation 1 (Figure 10G, blue). This is repeated for the remaining subpopulations in the excitatory and inhibitory populations. After every subpopulation has been stimulated over period T, a new random stimulation order is assigned for each subpopulation. The efficiency of the CR approach is shown in Figure 10A. Here, the CR stimulation depresses the average synaptic weight of the network over the course of ∼100 s. This causes a drop in synchrony from R(t) = 0.7 to R(t) = 0.05, which is shown in Figure 10C.
We integrated this strategy with our FTSTS strategy and applied it to an E-I network, consisting of 2,000 neurons, in the presence of E-to-E and I-to-I synaptic connectivity. We randomly divided each excitatory and inhibitory population into four subpopulations. Then, we adjusted the CR stimulation pattern to incorporate our FTSTS protocol (FTSTS-CR), such that each randomly selected pair of excitatory and inhibitory subpopulations are forced to spike post-before-pre. The FTSTS-CR stimulation pattern for one period is shown in Figures 10F,H. For a single FTSTS-CR pulse, we set U stim = 100 mV, T stim = 1 ms and T neutral = 7 ms. Similar to CR, we repeated this stimulation protocol for the other subpopulations in a random sequence for one period. Then, a new stimulation order is assigned for the next period T. Our protocol forces a randomly selected inhibitory subpopulation to fire prior to a randomly selected excitatory subpopulation, which over the course of 60 s depresses the average E-to-I synaptic weight of the entire population into the asynchronous regime as shown in Figure 10B. When the average E-to-I synaptic weight reached a preset value of 75 mV (asynchronous regime), we turned off the stimulation protocol. Figure 10D shows that the network remained in the desynchronized state at a synchrony level of R(t) = 0.05 for the rest of the simulation after removal of the FTSTS-CR protocol.
FIGURE 9 | Efficacy of the FTSTS strategy in desynchronizing 2,000 neuron E-I network where 25% of the population is inseparable and receives both the excitatory and inhibitory population input. (A) Shows the change in the average E-to-I and I-to-E synaptic weight of the network (black line), which is compared to the completely separable case considered in Figure 3 (red line). (B) Shows the network synchrony level during the simulation. The FTSTS pulse parameters are U stim = 100 mV, T stim = 1 ms, and T neutral = 10 ms.

Robustness to Additional Plastic Synapses
In previous sections, we presented our results for networks where we considered only E-to-I plastic synapses. In this section, we demonstrate the efficacy of our FTSTS strategy in desynchronizing networks where both E-to-I and I-to-E synaptic connections are plastic. We modeled the plasticity dynamics of I-to-E synapses using an anti-Hebbian STDP plasticity rule (Bell et al., 1997;Luz and Shamir, 2012). For anti-Hebbian STDP, we used (Equations 8,9) with the changed parameters a LTD = 1 and a LTP = −1.1 so that pre-before-post spike times decrease and post-before-pre spike times increases the synaptic weight. Our simulation results (see Figures 11A,B) show that the FTSTS strategy can potentially desynchronize the network by depressing both the E-to-I and I-to-E synaptic weights. Additionally, Figure 11C shows a decrease in the network synchrony from R(t) = 0.82 to R(t) = 0.05.

Robustness to a Symmetric Plasticity Rule
We demonstrate how a modified FTSTS protocol is able to control the synaptic weight of an E-I network with a symmetric plasticity rule. We use the same network described in section 2 with the following changes to the plasticity rule. The E-to-I synaptic weight is govern by the following equation Tass and Hauptmann, 2009): where ISI IE is the inter-spike-interval between spike-times of inhibitory and excitatory neurons in an E-to-I synapse. c P and c D are the potentiation and depotentiation learning rates, respectively. τ P is the potentiation time constant and τ D is the depotentiation time constant. The symmetric plasticity parameters can be found in Table 1.
In order to apply our FTSTS protocol for desynchronizing E-to-I network with a symmetric plasticity rule, we modified the FTSTS protocol by offsetting the pulse to the inhibitory population by T offset . This forced ISIs that promote either an increase or decrease in the synaptic weight. Using our modified FTSTS protocol, we show in Figures 12A,C that our FTSTS strategy can efficiently depress the synaptic weight of the synchronous network and desynchronize the network by forcing larger ISI IE values. Additionally, our protocol can also increase the average synaptic weight of the network to resynchronize the network by forcing short ISI IE values. Figure 12B shows the increase in the average E-to-I synaptic weight. The subsequent increase in the synchrony level is shown in Figure 12D, where the order parameter increases from R(t) = 0.05 to R(t) = 0.6.

DISCUSSIONS
In this paper, we developed and presented a novel stimulation strategy "Forced Temporal Spike-Time Stimulation" (FTSTS) for controlling synchronous activity of neurons in large spiking neural networks. Compared to other desynchronization strategies for large-scale spiking neural networks reported in the literature, our strategy focuses on controlling the average network synaptic weight by harnessing synaptic plasticity using a Hebbian-based spike-timing dependent plasticity (STDP) protocol that as a result controls the synchronization of neurons within the network. We presented a two neuron excitatoryinhibitory (E-I) network as an example to provide a mechanistic understanding of our approach. We later demonstrated the efficacy and robustness of the FTSTS strategy on large networks by varying the model parameters, synaptic connectivity and noisy inputs to the network. These results are also summarized in Figure 13 for clarity. While we only considered a LIF network, our method will lead to similar outcomes qualitatively for other neuronal models, since our approach is based on Hebbian activity-dependent plasticity. One of the prominent features of our FTSTS strategy is that it allows both synchronization and desynchronization of network activity by reversing the stimulation protocol (see Figures 3, 4), thus provides a complete control over the synchronization level of neural activity within a given network.
Our FTSTS strategy differs from existing stimulation strategies for desynchronizing spiking neural networks in many ways. Our strategy is based on harnessing the underlying synaptic plasticity compared to most of the desynchronization strategies reported in literature Popovych et al., 2005Popovych et al., , 2017Deuschl et al., 2006;Kiss et al., 2007;Nabi et al., 2013b;Mauroy et al., 2014;Wilson and Moehlis, 2014;Vlachos et al., 2016;Monga et al., 2018). Most of these strategies ignore the inherent synaptic plasticity among neurons in the network in designing the stimulation protocol for desynchronizing the network activity  (D) Shows the change in the synchrony of the network before, during, and after stimulation. The FTSTS parameters used to decrease the average E-to-I synaptic weight are U stim = 200 mV, T stim = 1 ms, T neutral = 22 ms, and T offset = 11 ms. The FTSTS parameters used to increase the average E-to-I synaptic weight are U stim = 100 mV, T stim = 1 ms, T neutral = 10 ms, and T offset = 5 ms.
(One exception is "Coordinate Reset" (CR) (Tass, 2003a,b;Tass and Majtanik, 2006;Tass and Hauptmann, 2007;Pfister and Tass, 2010;Ebert et al., 2014;Zeitler and Tass, 2015). As a result, these strategies effectively desynchronize the network activity if the stimulation protocol is active. Once the stimulation protocol is turned off, the network resynchronizes rapidly because of the disappearance of the asynchronous regime in the absence of stimulation. Our strategy alleviates this problem, like CR, by explicitly incorporating and harnessing Hebbian-based STDP, which allows the network to stay in the asynchronous regime for a longer-time period after the stimulation is turned off (see Figure 3).
Almost all the stimulation strategies focus on desynchronizing the network activity by randomizing the firing patterns (D) Efficacy of the FTSTS strategy in desynchronizing 2,000 neuron E-I network in the presence of E-to-E and I-to-I synaptic connectivity. (E) Efficacy of the FTSTS strategy in desynchronizing 2,000 neuron E-I network where 25% of the population is inseparable and receives both the excitatory and inhibitory population input. The FTSTS pulse parameters for all the studies are U stim = 100 mV, T stim = 1 ms, and T neutral = 10 ms. of neurons through direct stimulation. In comparison, our FTSTS strategy focuses on decreasing the average synaptic weight of the network by taking advantage of the Hebbianbased STDP protocol, which leads to the desynchronization of network activity. For example, the CR-based stimulation strategy desynchronizes the network activity by forcing different subpopulations of neurons to fire out of phase with each other, which resets the phase and desynchronizes the network (Tass, 2003a,b;Tass and Majtanik, 2006;Tass and Hauptmann, 2007;Pfister and Tass, 2010;Ebert et al., 2014;Zeitler and Tass, 2015). This generates an artificial asynchronous firing pattern that increases the basin of attraction of the asynchronous regime (i.e., lower synaptic weight stability point) (Pfister and Tass, 2010;Popovych and Tass, 2014). The underlying synaptic plasticity within the network then drives the average synaptic weight of the network toward the lower synaptic weight stability point (see Figures 10, 14 for comparison of our approach to the CR-based stimulation strategy).
Our developed framework can be incorporated into other desynchronization strategies, such as CR, to improve their efficacy. Figure 10 shows a comparison between the FTSTS-CR and CR performances in desynchronizing a E-I network consisting of 2,000 neurons with E-I synaptic plasticity (see section 3.8 for details of model parameters and specifics about the design of FTSTS-CR stimulation strategy). Since the FTSTS-CR stimulation strategy focused on decreasing the average synaptic weight of network, which as a result desynchronized the neural activity (see Figures 10B,D), this strategy outperformed the CR stimulation strategy shown in Figures 10A,C. One of the limitations of the CR stimulation strategy is that the long-lasting effects occur only in networks where the long-term depression (LTD) dominates the long-term potentiation (LTP) of the synapses on average so that the network exhibits bistability (Pfister and Tass, 2010). It has been found that LTP dominates in specific aberrant neuronal pathways and brain regions such as the striatum indirect pathway underlying Parkinson's disease and hippocampus underlying epilepsy (Mathern et al., 1998;Johnston, 2004;Shen et al., 2008). In such brain networks, the CR and FTSTS-CR stimulation strategy would both fail to produce long-lasting desynchronization of the network activity but the FTSTS-CR will have a longer acute desynchronization effect compared to the CR protocol. To demonstrate this, we applied the CR stimulation strategy to an E-I network consisting of 2,000 neurons where LTP dominates LTD. As shown in Figure 14D, the network synchrony level is disrupted during the stimulation period of 140 s but increased to a synchrony level of R(t) = 0.7 after removal of the stimulus. Since CR stimulation, in this scenario, only induces acute desynchronization of the network and does not reduce the average synaptic weight, the average synaptic weight of the network remains in the synchronous regime the entire time, as shown in Figure 14B. As a result, the network resynchronized rapidly after the removal of the CR stimulus. We compare the desynchronization efficacy of the FTSTS-CR and CR approach in Figures 14A,B, respectively. During the period of FTSTS-CR stimulation, there is a decrease in the average synaptic weight (see Figure 14A). This results in a reduction in the synchrony level to R(t) = 0.05 when the FTSTS-CR stimulus is removed. While this desynchronization is transient due to the domination of LTP, the network remains desynchronized for a longer period of time compared to CR as shown in Figure 14C.
In this work, we have considered excitatory-inhibitory (E-I) networks with plastic E-to-I synapses. In general, our approach is applicable to other types of spiking neural networks such as purely excitatory or inhibitory networks as well as to networks with other plastic synapses such as E-to-E or I-to-I synapses. One of the limitations of our approach is that it assumes the same stimulus waveform is delivered to individual neurons within a subpopulation. Although we have demonstrated in simulation that our FTSTS strategy effectively desynchronizes the neuronal firings in a network even when the stimulation waveform parameters for individual neurons are drawn randomly from a given distribution (see Figure 6) and are contaminated with input designed for the opposite population (Figure 9), it is still able to utilize the relationship between pre and post firings to effectively harness the synaptic plasticity. Multi-laser optogenetics and recent development in optogenetics to excite or inhibit the same neuron using two different light wavelengths could potentially alleviate this limitation for experimental implementation of our strategy (Forli et al., 2018). Additionally, we assumed the majority of the excitatory and inhibitory neuron populations were spatially separate, which allows for the neuron populations to be separately stimulated. Two examples of spatially separate excitatory and inhibitory neuron populations are the striatum and cortex or the globus pallidus external segment (GPe) and subthalamic nucleus (STN) (Lanciego et al., 2012;Hegeman et al., 2016). The GPe-STN network has traditionally been targeted for DBS-HFS to treat PD and could be potential area to test this hypothesis. Although we have not optimized the FTSTR pulses to achieve a better performance or to make it more energy efficient, it is not difficult to formulate optimization problems that minimizes the average synaptic weight, network synchrony level, and applied stimulation energy simultaneously to achieve a better overall performance.

DATA AVAILABILITY
All datasets generated for this study are included in the manuscript and/or the Supplementary Files.

AUTHOR CONTRIBUTIONS
GK and JS formulated the problem, discussed the results, and wrote the manuscript. JS developed the model, designed the FTSTS protocol, and ran the simulations.