(C) PLOS One This story was originally published by PLOS One and is unaltered. . . . . . . . . . . Synaptic reshaping of plastic neuronal networks by periodic multichannel stimulation with single-pulse and burst stimuli [1] ['Justus A. Kromer', 'Department Of Neurosurgery', 'Stanford University', 'Stanford', 'California', 'United States Of America', 'Peter A. Tass'] Date: 2022-11 In the present paper, we study synaptic reshaping by periodic multichannel stimulation (PMCS) in networks of leaky integrate-and-fire (LIF) neurons with spike-timing-dependent plasticity (STDP). During PMCS, phase-shifted periodic stimulus trains are delivered to segregated neuronal subpopulations. Harnessing STDP, PMCS leads to changes of the synaptic network structure. We found that the PMCS-induced changes of the network structure depend on both the phase lags between stimuli and the shape of individual stimuli. Single-pulse stimuli and burst stimuli with low intraburst frequency down-regulate synapses between neurons receiving stimuli simultaneously. In contrast, burst stimuli with high intraburst frequency up-regulate these synapses. We derive theoretical approximations of the stimulation-induced network structure. This enables us to formulate stimulation strategies for inducing a variety of network structures. Our results provide testable hypotheses for future pre-clinical and clinical studies and suggest that periodic multichannel stimulation may be suitable for reshaping plastic neuronal networks to counteract pathological synaptic connectivity. Furthermore, we provide novel insight on how the stimulus type may affect the long-lasting outcome of conventional DBS. This may strongly impact parameter adjustment procedures for clinical DBS, which, so far, primarily focused on acute effects of stimulation. Synaptic dysfunction is associated with several brain disorders, including Alzheimer’s disease, Parkinson’s disease (PD) and obsessive compulsive disorder (OCD). Utilizing synaptic plasticity, brain stimulation is capable of reshaping synaptic connectivity. This may pave the way for novel therapies that specifically counteract pathological synaptic connectivity. For instance, in PD, novel multichannel coordinated reset stimulation (CRS) was designed to counteract neuronal synchrony and down-regulate pathological synaptic connectivity. CRS was shown to entail long-lasting therapeutic aftereffects in PD patients and related animal models. This is in marked contrast to conventional deep brain stimulation (DBS) therapy, where PD symptoms return shortly after stimulation ceases. In this computational and theoretical study, we study periodic multichannel stimulation, a stimulation technique that allows for manipulating selected synaptic populations. Stimulus trains are delivered to multiple neuronal subpopulations in order to trigger neuronal responses. Using a model network of leaky integrate-and-fire neurons with spike-timing-dependent plasticity and theoretical analysis, we show how the relative timings between and the shape of delivered stimuli can be tuned to down-regulate certain synaptic connections while up-regulating others. Single-pulse stimuli triggered precise neuronal responses and were suitable for inducing a variety of synaptic network structures. When burst stimuli were employed, tuning the intraburst frequency allowed for distinguishing between down- and up-regulation of synaptic connections within individual neuronal subpopulations. Our work provides a theoretical basis for selecting suitable stimulation parameters for inducing long-lasting therapeutic effects in patients suffering from neurological disorders. Competing interests: I have read the journal’s policy and the authors of this manuscript have the following competing interests: PAT works as a consultant for Boston Scientific Neuromodulation. JAK declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. Funding: We gratefully acknowledge funding support by Boston Scientific Neuromodulation (Stanford Project 127674, url: https://www.bostonscientific.com/en-US/about-us/core-businesses/neuromodulation.html ), by the John A. Blume Foundation, and by the Foundation for OCD Research (New Venture Fund 011665-2020-08-01, url: https://www.ffor.org/ ). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The present paper is organized as follows. First, we present our analysis of the statistics of stimulus-triggered spikes for single-pulse and burst stimuli in networks of LIF neurons. Then, we provide a systematic theoretical and computational analysis of synaptic reshaping during PMCS. Based on this analysis, we present stimulation strategies for inducing a variety of network structures and test them in the LIF network model. Then, we discuss our results and their potential implications for DBS. Finally, we present the details of the LIF network model and methods used throughout the paper as well as the derivation of our theoretical approximations. We provide a detailed computational and theoretical analysis of PMCS, which led to promising stimulation strategies for possible clinical applications. We derive a theoretical basis for studying synaptic reshaping for a broad class of stimuli that may trigger complex neuronal spiking responses, such as bursts. Our theoretical predictions are compared to numerical simulations of PMCS of networks of leaky integrate-and-fire (LIF) neurons with STDP using electrical single-pulse and burst stimuli. We found that the same PMCS pattern administered using different types of stimuli leads to qualitatively different network structures. This suggests that the selection of the stimulus shape may not only impact the recruitment of neuronal subpopulations but also influence how stimulation affects plastic synaptic connections. Furthermore, we show how phase shifts between stimulus trains can be tuned to down-regulate synaptic connections between certain subpopulations while up-regulating others. Our results provide testable hypotheses for possible pre-clinical and clinical studies. For instance, PMCS could be delivered to down-regulate pathological synaptic connections in order to induce long-lasting therapeutic effects in patients suffering from neurological disorders, e.g., PD [ 23 , 26 ]. A: During PMCS, periodic stimulus trains (red) are delivered to multiple neuronal subpopulations (1,2,3) with phase shifts (Δα 1 ,Δα 2 ). D: Schematic of PMCS delivery using a multisite DBS electrode where phase-shifted stimulus trains are delivered to three contacts (1,2,3). Neurons are mostly affected by stimulus trains that are delivered to nearby stimulation contacts, neuron i by stimuli delivered to contact 1 and neuron j by stimuli delivered to contact 2. B: For illustration, in this figure we assume that neurons respond with a single spike (black vertical bar) of high fidelity to a stimulus (red). During stimulation, the spike train of the postsynaptic neuron, j, and that of the presynaptic neuron, i, become time shifted by Δα 1 /f, where f is the stimulation frequency. Due to STDP, time lags, Δt, between postsynaptic spikes and delayed presynaptic spike arrivals (gray) lead to weight updates, W, according to the STDP function (C, see Eq 10 ). Negative weight updates result from pairings of presynaptic spike arrivals with earlier postsynaptic spikes (Δt − < 0) and positive updates from pairings with later postsynaptic spikes (Δt + > 0). E: The sign of the overall weight update per inter-stimulus interval (ISI), , determines whether the synapse between neurons i and j is up-regulated (thick arrow) or down-regulated (thin arrow). F: If the neurons are bidirectionally coupled, PMCS may either induce effective unidirectional coupling (a,b), bidrectional coupling (d), or lead to decoupling (c). Blue dashed lines in panels E and F separate parameter regions with from regions with . Dark blue dashed lines separate corresponding regions for . In the present paper, we deliver periodic multichannel stimulation (PMCS) with different types of stimuli to induce a desired network structure in neuronal networks with STDP. During PMCS, phase-shifted periodic stimulus trains are delivered to multiple neuronal subpopulations. Depending on the phase shifts between stimulus trains delivered to synaptically interconnected subpopulations, postsynaptic neurons tend to spike either before or after presynaptic ones. Due to STDP, stimulation may then either induce bidirectional coupling, unidirectional coupling, or lead to a decoupling of the neuronal subpopulations. This basic mechanism of PMCS is illustrated in Fig 1 . The PMCS-induced synaptic network structure is determined by the phase shifts between the stimulus trains. PMCS can be understood as a generalization of classic CRS, with fixed stimulus sequence [ 17 , 19 ], as it allows for arbitrary phase shifts between stimulus trains. In order to identify stimulation patterns that may induce favorable changes of synaptic connectivity, a recent theoretical study presented an approximation of the mean rate of synaptic weight change during ongoing delivery of spatio-temporal stimulus patterns [ 60 ]. The results were derived in the limit of stimulation-controlled spiking activity, where neuronal spiking follows the stimulus pattern [ 60 ]. In that limit, synaptic reshaping is controlled by two main aspects of the stimulation pattern: the spatio-temporal correlations between stimulus deliveries and the statistics of neuronal spiking responses to individual stimuli [ 60 ]. The former aspect is controlled by the delivered stimulus pattern as characterized by the timings and locations of individual stimulus deliveries. Possible patterns include periodic stimulation, CRS, or randomized stimulus patterns. Previous studies have investigated synaptic reshaping for several stimulus patterns by assuming sharp distributions of stimulus-triggered neuronal spikes [ 60 – 64 ]. Less is known about how the statistics of stimulus-triggered spikes impact synaptic reshaping. These statistics are strongly affected by the stimulus type. In previous computational and theoretical studies, stimulation was delivered using charge-balanced electrical single-pulse stimuli that resembled direct electrical stimulation of the neurons’ somata [ 60 – 63 ]. Other computational studies analyzed synaptic reshaping during stimulation employing electrical burst stimuli [ 18 , 54 – 56 , 58 ], sensory stimuli, which resembled stimulation of neuronal fibers [ 56 , 65 ], and vibrotactile stimuli, modeling neuronal responses to vibratory stimulation of the fingertips [ 26 ]. It is currently unknown how such complex stimuli affect the stimulation-induced network structure and whether certain stimulus types are more favorable to induce synaptic reshaping than others. Of particular interest is the response of plastic neuronal networks to external stimulation. Theoretical and computational studies investigated synaptic reshaping in response to periodic input [ 52 , 53 ]. Other studies found that collective spiking events, for instance, in response to a stimulus, strongly impact synaptic reshaping. These events led to qualitatively different synaptic reshaping in neuronal networks with dendritic delays [ 28 ] and neuronal networks with axonal delays [ 29 ]. Such reshaping was also analyzed in networks with both axonal and dendritic delays [ 30 ]. In networks with both types of delays, stimulation resulted in the formation of distinct network motifs [ 44 ]. Another series of computational studies focused on the capability of stimulation to induce long-lasting desynchronization by down-regulating pathological connectivity (decoupling). This decoupling of neuronal subpopulations may drive the network into the attractor of a stable desynchronized state and lead to long-lasting desynchronization effects that outlast stimulation [ 18 , 54 – 59 ]. In order to understand synaptic reshaping of plastic neuronal networks, theoretical and computational papers studied neuronal networks with spike-timing-plasticity (STDP) [ 27 – 33 ]. Synaptic reshaping due to STDP depends on the time lags between postsynaptic and presynaptic spikes [ 34 , 35 ], or higher-order statistics of the spike times [ 31 , 36 ]. The type of plasticity determines the interplay of neuronal activity and network connectivity [ 32 , 33 , 37 , 38 ]. In many brain areas, STDP strengthens synapses if the postsynaptic neuron repetitively spikes shortly after the presynaptic one and weakens synapses if the postsynaptic neuron spikes shortly before the presynaptic neuron [ 35 ]. STDP may stabilize certain neuronal activity patterns [ 39 , 40 ], and it may lead to the formation of several network motifs [ 32 , 41 – 45 ], including neuronal assemblies. The latter are believed to play a key rule in memory formation [ 46 ]. STDP may also lead to the coexistence of different stable states, such as synchronized, desynchronized, and cluster states [ 18 , 38 , 47 – 51 ]. One stimulation technique that was found to induce long-lasting symptom relief in PD is coordinated reset stimulation (CRS). CRS was computationally developed to specifically induce long-lasting desynchronization effects [ 17 , 18 ]. It is a multisite stimulation technique during which spatio-temporal stimulus patterns are delivered to multiple neuronal subpopulations [ 17 , 19 ]. CRS was tested in the 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP)-monkey model of PD using downscaled DBS electrodes [ 20 ]. Two hours of CRS delivered to the STN on five consecutive days led to a significant reduction of akinesia, outlasting stimulation by several weeks [ 20 ], see also [ 21 , 22 ]. In PD patients, CRS unilaterally administered to the STN led to significant bilateral motor improvement, accumulating over consecutive days [ 23 ]. In contrast, symptoms return shortly after cessation of high-frequency DBS, the current standard of care for medically refractory PD [ 20 , 24 ]. Recently, CRS was also delivered using non-invasive vibrotactile fingertip stimulation and clinically significant long-lasting therapeutic effects were observed [ 25 , 26 ]. Currently, it is hypothesized that the long-lasting therapeutic effects of CRS result from stimulation-induced synaptic reshaping in related brain areas, driving the neuronal network into more favorable attractors, in this way restoring physiological function after cessation of stimulation [ 18 ]. Brain stimulation may induce long-lasting changes of synaptic connectivity by utilizing synaptic plasticity. Corresponding direct and indirect evidence has been provided by several experimental studies [ 10 – 16 ]. For instance, in rat brain slices, high-frequency stimulation of the subthalamic nucleus (STN), a common target area for deep brain stimulation (DBS) in PD, induced changes of synaptic weights that outlasted stimulation for at least about half an hour (the maximum observation window considered in the experiments) [ 10 , 14 ]. Another study showed that synaptic weight changes are strongly affected by the applied pattern of electrical stimuli and the concentration of neuromodulators, such as dopamine [ 11 ]. Indirect evidence of stimulation-induced synaptic plasticity has been provided by experiments on transcranial magnetic stimulation of the visual cortex [ 12 , 13 ]. There, aftereffects remained stable for hours after cessation of stimulation. Counteracting pathological synaptic connectivity using brain stimulation may improve current treatments and enable to develop novel treatments that aim at long-lasting symptom relief for patients suffering from neurological disorders. Several brain disorders are accompanied by synaptic dysfunction and/or pathological synaptic connectivity [ 1 ]. Such disorders include Alzheimer’s disease [ 2 ], attention-deficit/hyperactivity disorder [ 3 , 4 ], obsessive compulsive disorder (OCD) [ 5 , 6 ] and Parkinson’s disease (PD) [ 7 , 8 ]. Pathological connectivity may arise from different biological mechanisms that may, for instance, lead to the loss of synapses in disease-specific brain areas and/or to an impairment of synaptic plasticity [ 1 , 7 , 9 ]. A: Raster plot of simulated spiking activity of LIF neurons 1000 sec after onset of PMCS with burst stimuli with high intraburst frequency and short phase lags between subsequent subpopulations. B: Snapshots of simulated synaptic weight matrix after 1000 sec of stimulation (sim) and the corresponding theoretical prediction obtained from Eq 1 (theory). C: Schematic of stimulation-induced network structure (see caption of Fig 9 for symbols). A’-C’: Same as A-C but for single-pulse PMCS. One phase lag was chosen to be significantly smaller than the others in order to up-regulate only one population of synapses. Parameters: d e = 1 (all panels). In panels A-C, we used A stim = 0.8, an intraburst frequency of 120 Hz, five pulses per burst, f = 5 Hz, and Δα 1/2 = 0.1. In panels A’-C’, we used single-pulse stimuli with A stim = 0.4, a stimulation frequency of f = 10 Hz, and the phase lags Δα 1 = 0.1 and Δα 2 = 0.5. A: Raster plot of simulated spiking activity in the LIF network model 1000 sec after onset of PMCS with single-pulse stimuli and small positive phase lags between subsequent subpopulations. B: Snapshots of simulated synaptic weight matrix after 1000 sec of stimulation (sim) and the corresponding theoretical prediction obtained from Eq 1 (theory). C: Schematic of stimulation-induced feed-forward structure (see caption of Fig 9 for symbols). A’-C’: Same as A-C but for slow PMCS with long single-pulse stimuli. The phase lags were adjusted to cause a circular network structure. Parameters: A stim = 0.4 and Δα 1/2 = 0.3 (A-C), and Δα 1/2 = 0.33 (A’-C’) and d e = 1 and f = 10 Hz (A-C), and d e = 20 and f = 2.5 Hz (A’-C’). A: Raster plot of the spiking activity in the LIF network model 1000 sec after onset of PMCS with single-pulse stimuli and small positive phase lags between stimuli delivered to subpopulation two and subpopulations one and three. B: A snapshot of the synaptic weight matrix after 1000 sec of stimulation (sim) and the corresponding theoretical prediction obtained from Eq 1 (theory). C: Schematic of stimulation-induced tree network structure (see caption of Fig 9 for symbols). A’-C’: Same as A-C but for single-pulse PMCS with small negative phase lags between stimuli delivered to subpopulation two and subpopulations one and three. Parameters: A stim = 0.4, f = 10 Hz, and d e = 1 (all panels) and Δα 1 = 0.7, and Δα 2 = 0.3 (A-C), and Δα 1 = 0.3, and Δα 2 = 0.7 (A’-C’). A: Raster plot of simulated spiking activity for PMCS with burst stimuli with high intraburst frequency. B: Snapshots of simulated synaptic weight matrix after 1000 sec of stimulation (sim) and the corresponding theoretical prediction obtained from Eq 1 (theory). C: Schematic of stimulation-induced network structure (see caption of Fig 9 for symbols). A’-C’: Same as A-C but for burst stimuli with low intraburst frequency. Parameters: A stim = 0.8, f = 5 Hz, and d e = 1 (all panels). In panels A-C, we used an intraburst frequency of 120 Hz, three pulses per burst, and Δα 1/2 = 0.33 (A-C), In panels A’-C’, we used 60 Hz, eight pulses per burst, and Δα 1/2 = 0.31 (A’-C’). Next, we analyzed the dynamics of weights of inter-population synapses. To this end, we chose non-zero phase lags (Δα k ≠ 0) between stimulus trains delivered to adjacent neuronal subpopulations. We tuned the phase lags to induce several network structures in the network of LIF neurons. We present the results for the following network structures: A: Raster plot of simulated spiking activity for PMCS with Δα 1/2 = 0 and burst stimuli with low intraburst frequency. The scale bar indicates a time window of 100 ms. Stimuli are illustrated by red blocks. B: Snapshots of the simulated synaptic weight matrix after 1000 sec of stimulation (sim) and the corresponding theoretical prediction obtained from Eq 1 (theory). C: Schematic of stimulation-induced network structure. Circles represent subpopulations one, two, and three. Dashed arrows illustrate down-regulated synaptic connections. A’-C’: Same as A-C but for burst stimuli with high intraburst frequency. Thick arrows in panel C’ illustrate up-regulated synaptic connections. Parameters: A stim = 0.8, f = 5 Hz, and d e = 1 (all panels). An intraburst frequency of 60 Hz and five pulses per burst were used in panels A-C, and an intraburst frequency of 120 Hz and five pulses per burst were used in panels A’-C’. Thus, periodic stimulation may either lead to a down-regulation of all (intra-population) synapses or to an up-regulation of all (intra-population) synapses, depending on the shape of employed stimuli. If burst stimuli are used, the intraburst frequency plays a key role for the dynamics of (intra-population) synapses: low intraburst frequencies result in a down-regulation of synapses, whereas high intraburst frequencies result in an up-regulation. Raster plots of exemplary spiking activity and snapshots of the synaptic weight matrix for the network of LIF neurons during periodic stimulation with two types of burst stimuli are shown in Fig 9 . Next, we compared results for single-pulse and burst stimuli with different numbers of pulses per burst. Corresponding simulation results are compared to theoretical approximations in Fig 8B and 8C . The dynamics of the mean synaptic weight strongly depended on the intraburst frequency. Specifically, for low intraburst frequencies bursts with more pulses per burst led to faster decoupling ( Fig 8B , panel for t = 20 sec). The opposite trend occurred for high intraburst frequencies. Here, only single-pulse stimuli and bursts with two pulses led to a reduction of the mean weight whereas more pulses per bursts led to a strengthening of synapses. First, we considered single-pulse stimuli. Differences between periodic stimulation employing either short or long single-pulse stimuli can be inferred from the data points for Δα 1/2 = 0 in Fig 6 . Short single-pulse stimuli often led to faster decoupling, i.e. a faster down-regulation of intra-population synapses, than long single-pulse stimuli (data points for t = 20 sec in Fig 6 ). In an earlier study, it was suggested that the dynamics of intra-population synapses is strongly affected by the shape of administered stimuli [ 60 ]. This effect was referred to as stimulus-induced reshaping of synaptic connectivity. In order to study the impact of the stimulus shape on the synaptic weight dynamics of intra-population synapses, we considered the case of periodic stimulation, i.e. Δα 1/2 = 0 (modulo one). During periodic stimulation all synapses can be considered as intra-population synapses as the entire neuronal population receives stimuli simultaneously. Of particular interest is how PMCS can be employed to induce a desired network structure, i.e., to down-regulate certain synaptic populations while up-regulating others. Following, we study how the PMCS parameters can be adjusted to modulated intra- and inter-population synapses such that a particular synaptic network structure is induced. For comparison, we also performed simulations for CRS with rapidly varying sequence (CRS RVS). During CRS RVS each subpopulation receives exactly one stimulus during a cycle period of duration 1/f. In contrast to CRS with fixed sequence, the sequence at which individual subpopulation receive stimuli is shuffled after each cycle period (see Fig 8A’’’ ). CRS RVS has been used in preclinical and clinical studies on CRS as a treatment for PD [ 20 , 23 ]. Its decoupling properties for single-pulse stimuli have been analyzed in great detail in Ref. [ 61 ]. For single-pulse stimuli, CRS RVS led to decoupling for a wide frequency range ( Fig 8B’’’ and 8C’’’ ). Its decoupling effects were more robust with respect to changes of the stimulation frequency than those of PMCS with Δα 1/2 = 1/M, which corresponds to CRS with fixed sequence (compare Fig 8 panels B”, C” and B’’’, C’’’). However, besides single-pulse stimuli, only burst stimuli with two pulses per burst led to a complete decoupling in the LIF network. Note that only a limited range of stimulation frequencies was studied for CRS RVS in order to avoid an overlap of stimuli delivered to the same subpopulation. For the PMCS pattern with Δα 1/2 = 1/M, decoupling was even less robust with respect to changes of the stimulation frequency. For high intraburst frequencies, we only found decoupling for single-pulse stimuli and burst stimuli with two pulses per burst and only in limited frequency intervals ( Fig 8C” ). For low intraburst frequencies, it was possible to fine-tune the stimulation frequency such that complete decoupling was achieved. Overall, complete decoupling occurred for less frequencies the more pulses per burst stimulus were used ( Fig 8B” ). PMCS employing other stimuli typically led to non-zero mean synaptic weight ( Fig 8B” and 8C” ). For the PMCS pattern with Δα 1/2 = 0.5, complete decoupling was obtained over a wide frequency range if either single-pulse stimuli or burst stimuli with two pulses per burst and high intraburst frequency were used ( Fig 8C’ ). For other stimuli, we found decoupling only in finite frequency intervals. For other stimulation frequencies, the mean synaptic weight approached non-zero values (see Fig 8B’ and 8C’ ). Thus, the decoupling effects of this PMCS pattern were less robust with respect to variations of the stimulation frequency than those of periodic stimulation. For periodic stimulation, a complete decoupling was achieved for single-pulse stimuli and for burst stimuli with two pulses per burst for a wide range of stimulation frequencies ( Fig 8A, 8B and 8C ). For bursts with more pulses, decoupling was only obtained for a low intraburst frequency (60 Hz, Fig 8B ). In contrast, for a high intraburst frequency (120 Hz), burst stimuli with more pulses per burst led to an increase of the mean synaptic weight during stimulation. The fastest decoupling was achieved by PMCS with burst stimuli and a low intraburst frequency ( Fig 8B , panel for t = 20 sec). A: Raster plot of spiking activity of LIF neurons after 1000 sec of periodic stimulation (Δα 1/2 = 0) using burst stimuli of three pulses per burst and an intraburst frequency of 60 Hz. Furthermore, snapshots of the synaptic weight matrix of the LIF network (sim) and the corresponding theoretical approximation obtained from Eq 1 (theory) are shown. Successful stimulation-induced decoupling led to synaptic weights that were close to zero (white). B, C: The mean synaptic weight after 20 sec and 1000 sec of periodic stimulation as function of the stimulation frequency, f, for single-pulse stimuli (1P) and burst stimuli with two (2P), three (3P), and four (4P) pulses per burst (see legend). Panels B and C show results for an intraburst frequency of 60 Hz (B) and 120 Hz (C), respectively. Markers show simulation results (averaged over five network realization) and curves theoretical approximations ( Eq 4 ). Vertical lines mark the frequency of the synchronous rhythm prior to stimulation (red dashed line) and the set of stimulation parameters corresponding to the raster plot shown in panel A (colored dotted line), respectively. A’-C’ and A”-C”: Same as A-C but for PMCS with Δα 1/2 = 0.5 (A’-C’) and Δα 1/2 = 0.33 (A”-C”), respectively. (a’’’-c’’’): Similar plots for CRS RVS. Parameters: d e = 1 (all panels). A stim = 0.8 for PMCS with burst and CR RVS with single-pulse and burst stimuli. A stim = 0.4 for PMCS with single-pulse stimuli (B, C, B’, C’, and B”, C”). In Fig 8 , we compare simulation results for the network of LIF neurons to theoretical predictions for the mean synaptic weight obtained from Eq 4 . Complete decoupling, 〈w(t)〉 ≈ 0, was possible for all three PMCS patterns if the stimulus type and the stimulation frequency were adjusted accordingly. In Figs 6 and 7 , fast PMCS employing single-pulse stimuli led to an overall reduction of the mean synaptic weight across the (Δα 1 , Δα 2 )-plane ( Fig 6A, 6B and 6D ). We also observed a reduction of the mean synaptic weight for burst stimuli with low intraburst frequency; however, this reduction was sensitive to the phase lags Δα 1 and Δα 2 ( Fig 7B and 7D ). For both stimuli, the fastest reduction of the mean synaptic weight was obtained for either periodic stimulation (Δα 1/2 = 0) or for PMCS patterns where two subpopulations received stimuli simultaneously, i.e., Δα 2 ≈ 1 − Δα 1 . We also want to highlight the PMCS pattern with Δα 1/2 = 1/M, which corresponds to the original CRS pattern with fixed sequence introduced in Refs. [ 17 , 19 ]. For these three types of PMCS patterns, we analyzed how the mean synaptic weight during stimulation depended on the stimulation frequency and the stimulus type. Motivated by a recent study on decoupling stimulation [ 60 ], we analyzed to which extend PMCS may be employed as decoupling stimulation, i.e., was able to reduce the mean synaptic weight substantially. We found that the decoupling properties of PMCS depended on the phase lags, the stimulation frequency, and the stimulus type. In a portion of the (Δα 1 ,Δα 2 )-plane, PMCS with burst stimuli with a large number of pulses led to larger 〈w(t)〉 than PMCS with single-pulse or burst stimuli with low numbers of pulses (compare Figs 6 and 7 ). This effect was most pronounced in the panel corners, where phase lags were small but non-zero and even the reduction of 〈w(t)〉 by PMCS with single-pulse stimuli was weak (see Fig 6 ). Panels A,C and B,D show results for PMCS with burst stimuli with high (120 Hz) and low (60 Hz) intraburst frequencies, respectively. Rows show results for burst stimuli with three (top, 3P) and five (bottom, 5P) pulses, respectively. Columns show the mean synaptic weight 20 sec (left), 100 sec (middle), and 1000 sec (right) after onset of PMCS. Color maps show theoretical predictions, obtained from Eq 4 , and filled circles resulted from simulations of the LIF network model. Simulated spike trains, snapshots of the synaptic weight matrix, and corresponding theoretical predictions of the block structure of the synaptic weight matrix, obtained from Eq 1 , are shown below for the parameter sets marked by colored crosses (see corresponding figure labels). Parameters: f = 5 Hz, A stim = 0.8, d e = 1. Results for 〈w(t)〉 during PMCS with burst stimuli are shown in Fig 7 . 〈w(t)〉 depended on the phase lags, Δα 1 and Δα 2 . Furthermore, the shape of burst stimuli strongly affected 〈w(t)〉. For bursts with large numbers of pulses, small variations of the phase lags strongly affected 〈w(t)〉 ( Fig 7 ). This parameter sensitivity was most pronounced for phase lags that where close but not equal to zero (modulo one) (panel corners in Fig 7 ). Faster stimulation typically led to a faster reduction of 〈w(t)〉 as more stimuli per time were delivered (compare columns for t = 20 sec in Fig 6 ). The theoretical approximations for 〈w(t)〉 ( Eq 4 ) are shown as color maps in Fig 6 and show good quantitative agreement with simulation results for fast stimulation and short pulses. For long pulses and fast stimulation, theoretical approximations accurately predicted 〈w(t)〉 after long stimulation, which depends on the signs of the mean rates of weight change for the individual synaptic populations rather than the exact values ( Eq 1 ). Deviations occurred for long pulse durations and slow stimulation for which stimulus-triggered spiking responses were more susceptible to synaptic input from other subpopulations. This led to deviations between simulation results and theoretical predictions. The latter were derived under the assumption of stimulation-controlled spiking. Theoretical predictions (color maps) obtained from Eq 4 after 20 sec and 1000 sec of PMCS are compared to simulation results for the LIF network model (filled circles). w 0 = 0.38 is the mean synaptic weight prior to stimulation. Labeled columns show results for different stimulation frequencies: f = 10.0 Hz (A, D), f = 5.0 Hz (B, E), and f = 2.5 Hz (C, F). Rows show results for short pulses (d e = 1, top) and long pulses (d e = 20, bottom) ( Methods ). Simulated spike trains, snapshots of the connectivity matrix, and theoretical predictions of the block structure of the synaptic weight matrix are shown below for the parameter sets marked by red crosses (see red figure labels). Parameters: A stim = 0.4. Results for single-pulse stimuli with different pulse durations and PMCS frequencies are shown in Fig 6 . Stimulation led to either a reduction or an increase of 〈w(t)〉 depending on the phase lags Δα 1 and Δα 2 and the pulse duration ( Fig 6 ). PMCS with short single-pulse stimuli typically reduced 〈w(t)〉 more than PMCS with long single-pulse stimuli (compare simulation results in different rows of Fig 6 ). However, this effect was rather weak for the considered range of pulse durations. Next, we compared the theoretical predictions for 〈w(t)〉 to simulation results for networks of LIF neurons with STDP. Prior to stimulation, the LIF network was prepared in a synchronized state with a mean synaptic weight of 〈w(t)〉 = w 0 . Thus, a reduction of 〈w(t)〉 below w 0 indicates an overall decoupling whereas 〈w(t)〉 > w 0 indicates overall strengthening of synapses. Note that stimulation may up-regulate certain groups of synapses while down-regulating others ( Fig 5 ). Theoretical predictions for 〈w(t)〉 after long stimulation durations for a single-pulse stimulus and a burst stimulus with three pulses per burst are presented in Fig 5B and 5D , respectively. lim t→∞ 〈w(t)〉 was invariant under the transformations (Δα 1 , Δα 2 ) → (Δα 2 , Δα 1 ) and (Δα 1 , Δα 2 ) → (1 − Δα 1 , 1 − Δα 2 ). These symmetries result from the homogeneous network structure, equally sized subpopulations, and the invariance of the mean synaptic weight under shuffling of the weight matrix’s components. Furthermore, lim t→∞ 〈w(t)〉 only depended on the phase lags modulo one as the addition of integers to the phase lags results in a shift of the PMCS pattern by multiples of the stimulation period, which did not affect 〈w(t)〉 after long stimulation durations if STDP was sufficiently slow. Stimulation-induced changes of 〈w(t)〉 are of particular interest in the context of decoupling stimulation [ 59 , 60 ]. Decoupling stimulation was suggested to down-regulate pathological synaptic connectivity to drive multistable plastic neuronal networks into the attractor of a stable weakly connected desynchronized state, such that desynchronization effects persist after cessation of stimulation. In our LIF network model, this would correspond to a weakening of excitatory connections [ 60 , 61 ]. Of particular interest in this context are PMCS patterns that lead to a substantial reduction of 〈w(t)〉. We systematically varied the phase lags, Δα k , characterizing the PMCS pattern, and analyzed the induced network structure at different times after stimulation onset. First, we focused our analysis on the mean synaptic weight, 〈w(t)〉 = ∑ i,j∈syn. w i→j (t)/N syn . The sum runs over all synapses. The indices i and j refer to the presynaptic neuron and the postsynaptic neuron, respectively. N syn is the total number of synapses in the network. Note that 〈w(t)〉 carries less information about the network structure than the block structure of the synaptic weight matrix. This is because different combinations of blocks of up- and down-regulated synaptic connections may result in similar values of 〈w(t)〉. Nevertheless an analysis of 〈w(t)〉 can uncover which PMCS patterns lead to the up-regulation (down-regulation) of a certain portion of these blocks. We showed that the PMCS pattern and the shape of employed stimuli affected the stimulation-induced network structure. This suggests that PMCS can be tuned to induce a desired network structure, i.e., it may up-regulate certain synapses while down-regulating others. A more detailed analysis is presented in Fig 5 . Here, simulation results are compared to theoretical estimates of the synaptic weight matrix after long stimulation durations. Single-pulse and burst stimuli resulted in qualitatively different network structures. Intra-population synapses were down-regulated for single-pulse stimuli ( Fig 5A ). In contrast, they were up-regulated for burst stimuli with three pulses per burst and an intraburst frequency of 120 Hz ( Fig 5C ). Furthermore, for the employed single-pulse stimuli and stimulation frequency, it was not possible to choose the phase lags Δα 1 and Δα 2 such that bidirectional coupling between different subpopulations was induced (up-regulation of both and ) ( Fig 5A ). However, this was possible for certain combinations of phase lags when burst stimuli were used instead. For instance, for Δα 1/2 ≈ 1 bidirectional coupling between subpopulations one and two was induced ( Fig 5C ). For the same spatio-temporal PMCS pattern, single-pulse and burst stimuli resulted in qualitatively different network structures. If the pattern illustrated in Fig 4A was delivered using single-pulse stimuli, all synaptic connections were down-regulated ( Fig 4F ). Thus, such PMCS led to a complete decoupling of the LIF neurons (〈w x→y 〉 ≈ 0). In contrast, if the same PMCS pattern was delivered using burst stimuli with three pulses per burst, a qualitatively different synaptic network structure was obtained ( Fig 4I ). Intra-population synapses (diagonal blocks of the synaptic weight matrix) and synapses between subpopulations one and three were up-regulated, while synapses between other subpopulations were down-regulated. For the PMCS pattern illustrated in Fig 4A’ , all but the synapses in the block were down-regulated during PMCS with single-pulse stimuli ( Fig 4F’ ). In contrast, most synapses were up-regulated when burst stimuli were employed ( Fig 4I’ ). Predicted block structure of the synaptic weight matrix after long PMCS employing single-pulse stimuli (A, B) and PMCS employing burst stimuli consisting of three pulses (C, D). A: lim x→y 〈w x→y (t)〉, obtained from Eq 1 , (color map) is compared to simulation results for the LIF network model after 1000 sec of stimulation (filled circles). B: Predicted mean synaptic weight lim t→∞ 〈w(t)〉 obtained from Eq 4 . C,D: Same as A,B but for burst stimuli with three pulses per burst and an intraburst frequency of 120 Hz. Parameters: A stim = 0.4, f = 10 Hz, d e = 1 (A,B), and A stim = 0.8, f = 5 Hz, d e = 1, and an intraburst frequency of 120 Hz (C,D). The PMCS-induced block structure of the synaptic weight matrix depended on the phase lags Δα 1 and Δα 2 . Simulation results for the LIF network model are compared to theoretical predictions for the limit of long stimulation durations (lim t→∞ 〈w x→y (t)〉, Eq 1 ) in Fig 5A (single-pulse stimuli) and Fig 5C (burst stimuli with three pulses per burst). The long-term limit for weights of synapses in the blocks along the diagonal of the synaptic weight matrix—intra-population synapses—were insensitive to the phase lags Δα 1 and Δα 2 . In contrast, the long-time limit of weights of synapses in the off-diagonal blocks—inter-population synapses—showed a pronounced dependence on the phase lags. For a more detailed analysis of the stimulation-induced network structure, we studied the block structure of . Following, refers to the block containing weights of synapses between the presynaptic neuronal subpopulation x and the postsynaptic subpopulation y (see Fig 4F ). Depending on the phase lags Δα 1 and Δα 2 , PMCS either up-regulated or down-regulated synapses between individual neuronal subpopulations. This resulted either in high or low weights in the corresponding blocks . Theoretical estimates of lim t→∞ 〈w x→y (t)〉 ( Eq 1 ) are presented in Fig 4C and 4C’ for the two PMCS patterns, respectively. Panels A-I show results for the PMCS pattern illustrated in panel A. Red rectangles mark stimuli. Individual stimuli were either charge-balanced single-pulse (D-F) or burst stimuli (G-I). B: Phase shifts, Δϕ xy , between stimulus trains delivered to the postsynaptic subpopulation y and presynaptic subpopulation x. C: lim t→∞ 〈w x → y (t)〉 obtained from Eq 1 for single-pulse stimuli (1P) and burst stimuli with three pulses (3P) (see panel I for color code). Panels D-F show a raster plot of simulated spiking activity of the LIF network model 1000 sec after stimulation onset (D), and snapshots of the network’s synaptic weight matrix taken 20 sec (E) and 1000 sec (F) after stimulation onset. After 1000 sec, acute effects of stimulation have fully developed. The scale bar refers to a 100 ms time interval. The block of the synaptic weight matrix is marked red in panel F. G-I: Same as D-F but for burst stimuli with three pulses per burst and an intraburst frequency of 120 Hz. A’-I’: Same as A-I but for the PMCS pattern illustrated in panel A’. Parameters: f = 5 Hz; A stim = 0.4 (single-pulse stimuli) and 0.8 (burst stimuli), d e = 1. Δα 1/2 = 0.5 (A-I), and Δα 1 = 0.1 and Δα 2 = 0.3 (A’-I’). Using the phase lags Δϕ xy and in Eq 1 , we obtained an approximation for the time-dependent mean synaptic weight 〈w(t)〉 during stimulation (4) Here, we assumed a homogeneous network with M equally sized subpopulations, as realized in our LIF network model. Note that while the mean synaptic weights of synapses interconnecting the individual subpopulations, 〈w x→y (t)〉, evolve approximately linear in time until they reach one of the hard bounds ( Eq 1 ), the overall mean synaptic weight, 〈w(t)〉, may possess a more complex dynamics ( Eq 4 ). Next, we considered PMCS with an arbitrary number of subpopulations M. Using , we approximated the mean weight 〈w x→y (t)〉 of synapses interconnecting the subpopulations x and y at time t after stimulation onset. Here and in the following, the indices x and y refer to the presynaptic and postsynaptic neuronal subpopulations, and the indices i and j to the presynaptic and postsynaptic neuron, respectively. Given 〈w x→y (t)〉 at stimulation onset t = t 0 , 〈w x→y (t)〉 at time t during stimulation (t > t 0 ) is approximately given by (1) Here, [x] clip,[0,1] = x for x ∈ [0, 1], [x] clip,[0,1] = 0 for x < 0, and [x] clip,[0,1] = 1 for x > 1 accounts for the hard bounds for individual synaptic weights ( Methods ). The function (2) describes a linear decay/increase of 〈w x→y (t)〉 during stimulation with the mean rate of weight change . Δϕ xy is the phase shift between stimulus deliveries to the postsynaptic subpopulation y and the presynaptic subpopulation x and follows from the phase lags Δα k characterizing the PMCS pattern ( Fig 1A ). Eq 1 was derived based on the assumption that the weights of individual synapses are close to the hard bounds prior to stimulation onset, i.e, either w i→j ≈ 0 or w i→j ≈ 1 at time t 0 . This is typically observed in networks with additive STDP and hard bounds, Ref. [ 37 , 67 ] and was observed in a previous study using a similar LIF model [ 60 ]. Eq 1 either accounts for a linear increase of the weights of weak synapses ( ) or a linear decrease of the weights of strong synapses ( ). In Fig 3 , we compare with . For single-pulse stimuli, synaptic weight changes mainly resulted for small positive and negative phase lags, Δα (modulo 1) ( Fig 3A and 3B ). For slow single-pulse stimulation, larger phase lags also contributed to weight changes ( Fig 3C ). For burst stimuli, weight changes occurred for a wide range of phase lags ( Fig 3A’, 3B’, 3C’, 3A”, 3B” and 3C” ). Faster stimulation and increasing the number of pulses per stimulus resulted in larger positive and negative weight updates. In addition, we derived a theoretical estimate for the mean rate of synaptic weight change, , using the distributions of stimulus-triggered spikes, Λ k (t). is given in Eq 19 in the Methods section. Δα is the phase lag between periodic stimulus trains delivered to the postsynaptic and presynaptic subpopulation and is given by the time shift between stimuli divided by the stimulation period ( Fig 1B ). ξ is the difference between axonal and dendritic delays. During the derivation of , we assumed that Λ k (t) approximated the spiking of neurons in individual subpopulations well, even if spiking of other subpopulations was altered by stimulation. Panels show results for single-pulse (A-C) and burst stimuli with three (A’-C’) and five (A”-C”) pulses per burst. Results for high (120 Hz) and low (60 Hz) intraburst frequencies are shown in red and blue, respectively. Curves show theoretical approximations ( Eq 19 ) and markers obtained from simulations and averaged over five network realizations and initial conditions ( Methods ). Parameters: d e = 1 and A stim = 0.4 (A,B,C), A stim = 0.8 (A’-C’, A”-C”). Note that different intervals are plotted on the y-axes. Fig 3. Effect of PMCS with (M = 2) on the mean weight of synapses with postsynaptic neuron in subpopulation 2 and presynaptic neuron in subpopulation 1 as function of the phase lag, Δα, between stimulus trains delivered to the two subpopulations. We considered the case of M = 2 equally sized subpopulations and delivered PMCS with phase lag Δα ≔ Δα 1 and stimulation frequency f to the network of excitatory LIF neurons. We recorded the trace of the mean weight of synapses between neurons in subpopulation one and neurons in subpopulation two, 〈w 1→2 (t)〉, and estimated the mean rate of weight change, ( Methods ). Simulation results for different types of stimuli and stimulation frequencies are shown in Fig 3 . Next, we studied how PMCS delivered to M neuronal subpopulations affects the strength of synaptic connections. For illustration, we restricted our analysis to M ≤ 3. The spatio-temporal characteristics of the PMCS pattern are controlled by the phase lags Δα k , k = 1, 2, ‥, M − 1, between stimulus deliveries to adjacent neuronal subpopulations ( Fig 1A ). Note that simulations were performed for networks of LIF neurons with homogeneous synaptic connectivity. Thus, spatial adjacency was solely induced by the grouping of neurons into subpopulations according to their indices. Neurons in the same subpopulation received stimuli simultaneously. For statistical analysis, we evaluated the distribution of the kth spike time during an ISI, Λ k (t), ( Methods ) and the corresponding cumulative distribution function . The shape of F k (t) for single-pulse stimuli depended on the pulse duration, d e . Short pulses typically resulted in step-like F k (t) whereas long pulses led to broad Λ k (t) resulting in a slower increase of F k (t), mostly during the excitatory part of the stimulus pulse ( Fig 2C ). For burst stimuli, individual stimulus pulses led to step-like increases of F k (t) ( Fig 2G and Methods ), indicating that spiking occurred shortly after the onset of the kth stimulation pulse ( Fig 2G ). A: Raster plot of spiking activity of a network of 10 3 LIF neurons during periodic single-pulse stimulation (red bars) of 333 neurons. B: Mean number of spikes per ISI and neuron as function of the dimensionless stimulation amplitude, A stim , and pulse width, d e ( Methods ). Gray contour lines correspond to 0.75 and 1.25 spikes per ISI and neuron (from bottom to top). C: Cumulative distribution function of the timing of the first spike of LIF neurons after stimulus delivery ( Methods ) for the stimulus parameters marked orange, black, and brown, respectively, in panel B. The inset zooms into the first 5 ms. The stimulation current, I stim , corresponding to d e = 10 is shown in panel D. E: Same as A but for burst stimuli with five pulses per burst and an intraburst frequency of 120 Hz. F: Mean number of spikes per ISI and neuron as function of the number of stimulus pulses per burst and the stimulation amplitude. G: Cumulative distribution function of the timing of the kth spike per ISI for the parameter set marked white in panel F. Color code: F 1 (t) (black), F 2 (t) (blue), F 3 (t) (green), F 4 (t) (brown), and F 5 (t) (red). H: Corresponding stimulation current, I stim . Parameters: f = 5 Hz (all); A stim = 0.4, and d e = 10 (A,C,D) and d e = 1 (C, orange) and d e = 20 (C, brown); and A stim = 0.8, d e = 1, and 120 Hz intraburst frequency (E, F, G, H). Results are shown in Fig 2 . Spiking of stimulated neurons entrained with stimuli ( Fig 2A and 2E ). Stimulus-triggered spiking responses depended on the shape of administered stimuli ( Fig 2C, 2D, 2G and 2H ). We measured the mean number of spikes per ISI ( Fig 2B and 2F ). Periodic single-pulse stimulation led to at most two spikes per ISI ( Fig 2B ). The first spike occurred shortly after stimulus onset. For sufficiently strong and broad stimuli, a second spike occurred towards the end of the ISI, when the LIF neurons’ refractory periods had past and neurons became susceptible to excitatory synaptic input from the rest of the network. In contrast, strong burst stimuli triggered several spikes per ISI, typically one spike per pulse ( Fig 2F ). We studied the response of networks of excitatory LIF neurons with STDP to PMCS. The dynamics of the LIF neurons’ membrane potentials was adjusted to that of tonically spiking neurons in the rat STN [ 66 ]. The networks were prepared in a stationary state with synchronized spiking and strong synaptic connections with mean synaptic weight 〈w〉 ≈ 0.38 ( Methods ). This resembled a pathological state in which strong synaptic connectivity supports excessive neuronal synchrony, e.g., during PD. Then, stimulation was delivered using either charge-balanced single-pulse stimuli, referred to as single-pulse stimuli in the following, or bursts of charge-balanced pulses, referred to as burst stimuli in the following (see Methods for details). Discussion Several neurological disorders are associated with impaired synaptic connectivity. We suggest PMCS to harness synaptic plasticity to counteract pathological synaptic connectivity. During PMCS, phase-shifted periodic stimulus trains are administered to different neuronal subpopulations. We studied synaptic reshaping during PMCS employing single-pulse and burst stimuli in networks of excitatory LIF neurons with STDP. We showed that: (a) The phase lags between stimulus trains can be tuned to up-regulate certain populations of synapses while down-regulating others. This way, a desired network structure can be induced. (b) The type of administered stimuli strongly affects the PMCS-induced network structure. (c) During PMCS with burst stimuli an adjustment of the intraburst frequency allows for either up- or down-regulation of synapses between neurons within the same subpopulation. (d) PMCS can be employed as a decoupling stimulation. In particular, PMCS patterns in which multiple subpopulations receive stimuli simultaneously led to decoupling effects that were robust with respect to changes of the stimulation frequency if appropriate stimuli were used. In our LIF model such stimuli included single-pulse stimuli and burst stimuli with low intraburst frequency. Our computational analysis of PMCS in networks of LIF neurons and our theoretical approximations revealed stimulation strategies for a controlled modulation of synaptic connectivity. PMCS may be suitable to counteract pathological synaptic connectivity in patients suffering from neurological diseases and may present a step towards novel treatments that lead to symptom relief that outlasts stimulation. Combined computational and theoretical approach to predict PMCS-induced synaptic reshaping To predict synaptic reshaping during PMCS, we presented a combined computational and theoretical approach, consisting of the following steps. (i) The statistics of stimulus-triggered neuronal spikes were estimated based on simulations of LIF neurons during periodic stimulation of one neuronal subpopulation. (ii) These statistics were used to estimate the mean rate of synaptic weight change during stimulation; and (iii), the latter estimate was used to predict the dynamics of the weights of synapses between individual neuronal subpopulations. Steps (i) and (ii) of our approach result in estimates of the mean rate of weight change. These estimates were compared to simulation results for the LIF network model in which two subpopulations received stimuli at a fixed phase lag, which resembled a standard STDP experiment [35]. The theoretical predictions approximated the simulation results well for fast stimulation for various stimulus types (Fig 3). Deviations for slow stimulation resulted from the assumption that the distribution of stimulus-triggered spikes does not change during PMCS, made during the derivation of Eq 19. This implicitly assumed weak synaptic interaction and strong stimulation. Accordingly, deviations of simulation results from our theoretical predictions occurred if spiking responses were perturbed by synaptic input from other subpopulations, mainly when this input resulted from stimulus-triggered collective spiking events (see, for instance, Fig 12A’, 12B’ and 12C’). During slow stimulation, the LIF neurons were close to their spiking thresholds for a substantial part of the ISIs. This made them more susceptible to synaptic input from other neuronal subpopulations. For burst stimuli with large numbers of pulses per burst, individual stimulus pulses triggered sharp collective spiking responses. These led to strong synaptic input to other neuronal subpopulations, which led to deviations (Figs 10A’, 10B’, 10C’, 13A, 13B and 13C). An alternative to estimating the mean rate of weight change based on the statistics of stimulus-triggered spikes (steps (i) and (ii)) would be to estimate the mean rate of weight change experimentally or from simulation data (just as we did in Fig 3). Then, these estimates can be used in step (iii) (Eq 1) to predict the outcome of PMCS in experiments or other computational models. Our approach can be applied to a wide range of computational neuron and network models. Furthermore, PMCS delivered to two subpopulations corresponds to a standard STDP experiment, where postsynaptic and presynaptic spikes are elicited by applying phase-shifted stimulus trains [34, 35]. The case of two subpopulations also applies to experiments on paired associative stimulation during which repetitive low-frequency median nerve stimulation was paired with transcranial magnetic stimulation (TMS) [68–70]. This idea has been extended to pair different types of stimuli in order to study synaptic plasticity in different brain areas. For instance, DBS stimuli were paired with cortical TMS stimuli in Refs. [15] and [16]. This suggests that our results may apply to various types of brain stimulation, including DBS, vibrotactile fingertip stimulation, and TMS. In the present paper, we applied our approach to networks of excitatory LIF neurons with STDP. While the LIF neurons responded with high fidelity to administered stimuli, our general approach assumes that neuronal spiking responses are entrained with the periodic stimulus trains. Such an entrainment of oscillators with periodic stimulation is a general phenomenon and has been observed in various biological systems [71]. For instance, the entrainment of neuronal spike rhythms to HF DBS was studied in target areas for the treatment of PD, i.e., in the basal ganglia and the thalamus. There, it was observed that neuronal activity in brain areas that received synaptic input from the stimulated target brain regions entrained with the stimulation [72, 73]. Individual spike times occurred at fixed phases of the HF DBS rhythm, with millisecond precision [72, 73]. Furthermore, an entrainment of cortical activity to HF DBS was observed, suggesting antidromic activitation of cortical neurons by STN HF DBS [74]. Besides DBS, several other stimulation techniques were found to cause an entrainement of neuronal activity with administered stimuli. For instance, during vibrotactile stimulation an entrainment of neuronal spiking responses with skin indentation oscillations was observed in the human somatosensory thalamus [75] and in the monkey primary somatosensory cortex [76, 77]. An entrainment of neuronal activity with stimulus deliveries was also reported by experimental studies on TMS [78] and transcranial alternating current stimulation [79], as well as visual stimulation [80]. In Ref. [81], tremor entrainment was achieved by periodic burst stimulation. Their computational results suggested that this was caused by an entrainment of tremor-related neurons with the stimulus pattern. Widely observed entrainment of neuronal activity with different types of periodic brain stimulation suggests that PMCS may be delivered using a variety of stimulation techniques and that our combined computational and theoretical approach may be applicable. Statistics of stimulus-triggered spikes and synaptic plasticity determine synaptic reshaping We delivered PMCS to two neuronal subpopulations. Our analysis revealed a strong dependence of the mean rate of synaptic weight change on stimulation parameters (Fig 3). Applying our combined computational and theoretical approach, we were able to explain this dependency by considering the impact of stimulation parameters on the statistics of stimulus-triggered spikes and calculating the resulting mean rate of synaptic weight change in neuronal networks with STDP. A strong dependence of synaptic reshaping on the statistics of stimulus-triggered neuronal responses has been observed experimentally. For instance, in the rat STN, Wang et al. observed that stimulus-triggered rebound bursts with more than two spikes per burst triggered long-term potentiation (LTP) of inhibitory synapses on STN neurons whereas bursts with two or less spikes caused long-term depression (LTD) or no change [82]. In another study on cortical neurons, a strong dependence of synaptic plasticity on the neuronal firing rate was observed [83]. Several types of synaptic plasticity have been observed in the brain. In the LIF network model and in our theoretical derivation of the mean rate of weight change, we considered a nearest-neighbor STDP scheme [84]. This presents a simplification as experimental studies found that the presynaptic and postsynaptic spike patterns affect the dynamics of synaptic weights [31, 36]. In general, weight updates during fast spike patterns, such as bursting, may not simply result from a superposition of weight updates for single spike pairs. One suggestion of how to account for this observation was to use STDP schemes that only consider pairings with the latest/next spike, so-called nearest-neighbor schemes [83, 84], comparable to the one we use in the present paper. Other studies suggested more complex multi-spike STDP schemes in which weight updates depend on the timing of triplets or quadruplets of spikes [85, 86]. In addition to the considered spike times, the effect of STDP strongly depends on the STDP function (Eq 10). An overview of experimentally observed STDP functions can be found in Ref. [87]. In our LIF network model, we considered canonical STDP [35, 67, 88] (see Fig 1C). Canonical STDP has been found in several brain regions including the neocortex [88], the hippocampus [35], and the striatum, where STDP was also modulated by dopamine [89]. While detailed information on the shape of the STDP function is often not available, several experimental and clinical studies reported brain stimulation-induced synaptic reshaping and long-lasting effects. In the context of PD, direct and indirect evidence for synaptic reshaping due to stimulation of target brain areas has been presented. Common target areas for therapeutic DBS in PD patients include the STN [90] and the internal segment of the globus pallidus (GPi) [91]. Stimulation-induced reshaping of synaptic connections in response to STN stimulation was observed in rat brain slices [10, 11, 14, 92]; however, while synaptic reshaping depended on the delivered stimulus pattern [11], to the best of our knowledge, only one study showed a dependence on the relative timings of postsynaptic and presynaptic spikes in these brain areas [14] and thereby provided evidence for STDP. Indirect evidence for STN stimulation-induced synaptic reshaping was presented by a study in PD patients [15]. There, STN DBS induced motor-cortical plasticity. In experiments STN DBS was paired with TMS of the primary motor cortex. Repetitive stimulation entailed long-lasting (≈ 45 min) changes of the amplitudes of motor evoked potentials (MEPs) [15]. Another study analyzed the effect of GPi DBS in a similar setup. Single DBS pulses were paired with motor-cortical TMS at different time lags [16]. The authors reported that repetitive delivery of paired stimuli entailed long-lasting changes (≈ 1 hour) of the amplitude of MEPs. Another line of indirect evidence indicates stimulation-induced synaptic reshaping during CRS. Long-lasting effects of CRS were observed in experiments in rat brain slices [93]; preclinical studies in MPTP monkeys, where CRS was delivered using down-scaled multisite DBS electrodes [20–22]; and in PD patients, where CRS was delivered through standard multisite DBS electrodes [23]. Recently, long-lasting therapeutic effects in PD patients were also induced by vibrotactile fingertip CRS [25, 26, 94]. There, multiple neuronal subpopulations were targeted by delivering phase-shifted stimuli to different fingertips. Neurons that respond to signals from cutaneous receptor afferents of different fingers are arranged somatotopically in the sensory cortex [95–97], which provides evidence that the resulting neuronal responses occurred primarily in segregated neuronal subpopulations related to input from the respective stimulated fingertips. Plasticity has also been induced using other types of brain stimulation. For instance, when phase-shifted stimuli were delivered to the two cortical hemispheres using transcranial alternating current stimulation [98]. The authors observed long-lasting changes in cortical connectivity. These changes were later described by a corresponding computational model incorporating STDP [99]. These lines of evidence suggest that various types of brain stimulation may indeed be suitable to realize the PMCS approach and induce long-lasting changes of synaptic connectivity. This may provide a tool to induce long-lasting therapeutic effects in patients suffering from neurological disorders associated with pathological synaptic connectivity. Stimulus shape impacts synaptic weight dynamics Of particular interest is the impact of the stimulus shape on the synaptic weight dynamics, which was referred to as stimulus-induced synaptic reshaping in Ref. [60]. Our combined computational and theoretical approach revealed that changes of the stimulus type had a similar effect on the mean rate of weight change as if the same stimulus was delivered to networks with different STDP functions (see Fig 3). Consequently, the type of employed stimuli had a strong impact on the PMCS-induced network structure. For instance, in Fig 4, the same PMCS pattern delivered using either single-pulse or burst stimuli induced qualitatively different network structures. Previous computational studies considered different stimuli during the delivery of CRS and studied long-lasting effects of stimulation. In Ref. [56], CRS was delivered using three different types of model stimuli: (a) bursts of charge-balanced electrical pulses, (b) a single excitatory postsynaptic potential (corresponding to stimulation of an incoming fiber of an excitatory neuron), and (c) a single inhibitory postsynaptic potential (stimulation of an incoming fiber of an inhibitory neuron). Case (a) was modeled by delivering electrical current pulses to stimulated neurons, similar to the electrical pulses considered here. The authors found that all three types of stimuli led to long-lasting desynchronization [56]. In their model, all three types of stimuli led to a phase reset of the stimulated neuronal subpopulation, which presumably resulted in a sharp distribution of stimulus-triggered spikes. Note that this distribution was not studied in detail; however, phase resetting was demonstrated numerically. Therefore, we speculate that the resulting distributions of stimulus-triggered spikes were rather similar even though different stimuli were used. Other computational studies delivered CRS using vibrotactile stimulation (see results for the computational model in [26]). Vibrotactile stimuli led to a broad distribution of stimulus-triggered spikes. Their results differed significantly from results for electrical model stimuli that led to a sharp distribution of stimulus-triggered spikes for a similar neuronal network model [61]. This is in line with our results that suggest that these differences are a consequence of the qualitatively different distributions of stimulus-triggered spikes for vibrotactile burst stimuli [26] and for electrical single-pulse stimuli [61]. More studies are needed to understand how different stimulus shapes can be employed to modulate synaptic connections using brain stimulation. In the context of DBS, most research has been devoted to tuning the waveform of DBS pulses for improving the efficiency in initiating neuronal responses [100–103], in entraining neuronal spiking [104], and for improving acute therapeutic effects [105, 106]. One study that, indirectly, analyzed how different stimulus patterns modified cortico-subthalamic synapses was performed by Yamawaki et al. [11]. There, in dopamine-intact conditions, high frequency stimulation hardly affected synaptic weights, whereas low-frequency bursts caused LTP of cortico-subthalamic synapses. Our results suggest that the waveform of DBS stimuli as well as the numbers of spikes per burst and the intraburst frequency during DBS may have a strong impact on synaptic connections and, consequently, potential long-lasting aftereffects, e.g. during CRS of the STN [20–23]. PMCS-induced decoupling of neuronal populations Recent studies suggested decoupling stimulation in order to induce long-lasting desynchronization [60–62]. This is important in the context of brain disorders characterized by abnormal neuronal synchrony, such as PD [107]. We analyzed the capability of PMCS with different types of stimuli to decouple neuronal networks. To this end, we investigated which PMCS parameters led to a reduction of the mean synaptic weight during PMCS and how robust this decoupling was with respect to changes of the stimulation parameters, e.g., the stimulus shape, the PMCS pattern, and the stimulation frequency (Figs 6–8). Decoupling during PMCS with short electrical single-pulse stimuli was robust with respect to variations of the PMCS pattern, characterized by the phase lags between stimulus trains delivered to different neuronal subpopulations (Fig 6). The fastest decoupling was observed when stimuli were simultaneously delivered to multiple neuronal subpopulations. In this case, PMCS utilized “decoupling by synchrony”, which was previously presented in Ref. [29]. Decoupling by synchrony occurs in networks with canonical STDP and longer axonal delays than dendritic delays during sharp collective neuronal spiking events [29, 30]. Based on the observation of a corresponding “coupling by synchrony” effect in networks with longer dendritic than axonal delays by Morrison et al. [28], we would expect that PMCS with simultaneous stimulus delivery to multiple subpopulations would strengthen synapses in such networks. When burst stimuli were employed, simultaneous stimulus delivery to multiple neuronal subpopulations led to even faster weight changes than for single-pulse stimuli. However, it depended on the intraburst frequency whether PMCS led to a weakening or a strengthening of synapses. PMCS with high intraburst frequencies increased the mean synaptic weight whereas low intraburst frequencies led to a reduction of the mean synaptic weight for a wide range of phase lags between stimulus trains (Fig 7). Overall the weight dynamics induced by PMCS with burst stimuli was more sensitive to the phase lags between stimulus trains than the weight dynamics induced by PMCS with single-pulse stimuli. The fastest decoupling was observed when stimuli were simultaneously delivered to all subpopulations (periodic stimulation) and burst stimuli with low intraburst frequency were used. For sufficiently low stimulation frequencies, such a PMCS pattern is comparable to theta burst stimulation. A recent study in PD patients studied therapeutic effects of STN theta burst DBS [108]. There, 5 Hz burst stimulation (100 ms stimulation ON followed by 100 ms stimulation OFF) was delivered to PD patients. The authors varied the intraburst frequency and observed pronounced acute effects on rigidity, tremor, and akinesia. While rigidity improved for low and high intraburst frequencies (50 Hz and 100 Hz, respectively), akinesia responded better to low, whereas tremor responded better to high intraburst frequencies. This provides evidence that the intraburst frequency has a strong impact on therapeutic effects. The authors also tested for long-lasting aftereffects that outlasted stimulation; however, they did not observe significant aftereffects. They discussed that this was due to a too short stimulation duration of about 20–30 min used in the experiments, which was substantially shorter than the 2 h duration used in the proof-of-concept study for CRS in PD patients [23]. In other studies, theta burst stimulation was delivered using repetitive TMS [70]. The authors observed long-lasting aftereffects on the amplitude of motor-evoked potentials and reaction time that depended on the burst duration (corresponding to the number of TMS pulses per burst). Long-lasting effects of periodic burst stimulation were also observed during DBS of the globus pallidus in 6-OHDA mice, an animal model for PD [109]. The burst stimulation was able to restore movement, whereas continuous HF DBS did not have such aftereffects. Showing that the stimulus pattern has a strong impact on therapeutic effects. A previous study by Mastro et al. [110] provided evidence that stimulation may have to induce neuron population-specific responses in order to induce long-lasting effects. Variations of the intraburst frequency of delivered burst stimuli led to modulations of these neuronal responses, suggesting an impact of the intraburst frequency on long-lasting effects of stimulation. However, it was not studied whether the observed long-lasting effects resulted from stimulation-induced synaptic reshaping. Together, these studies provide evidence that periodic burst stimulation may indeed be a promising candidate for inducing long-lasting effects in neuronal networks. While our results point out the importance of the intraburst spike pattern for harnessing STDP, the presented experimental evidence also suggests that it may be important for recruiting neuronal target populations and inducing neuronal responses that are beneficial for long-lasting effects. PMCS compared to CRS RVS The presented PMCS may be considered as a generalization of CRS with fixed sequence [17]. CRS with fixed sequence was used in earlier studies, originally, in networks without plasticity. There, CRS was administered with fixed sequence to induce cluster states and acute desynchronization [17, 19]. Later, CRS RVS was delivered in plastic neuronal networks to avoid the formation of sequence related clusters, which might occur during CRS with fixed sequence [18]. CRS RVS was also successfully administered in clinical and preclinical studies to induce long-lasting therapeutic effects in PD patients [23, 26] and related animal models [20], respectively. In Zeitler et al. CRS with slowly varying sequence was introduced [65]. During CRS with slowly varying sequence, the sequence of stimulus deliveries to the individual neuronal subpopulations is shuffled every n cycles. Thus, the variation of n allows for studying the transition between CRS RVS (n = 1) and CRS with fixed sequence (n → ∞). During our detailed analysis of the decoupling aspects of PMCS, we found that CRS with a fixed sequence (corresponding to PMCS with phase lags Δα k = 1/M) leads to a complete decoupling of the LIF network for intermediate stimulation frequencies when single-pulse or burst stimuli with a small number of pulses are used (Fig 8B” and 8C”). Remarkably, for the same stimuli, CRS RVS led to decoupling in an even broader frequency range (Fig 8B’’’ and 8C’’’). This provides evidence that random shuffling of the stimulation sequence is indeed advantageous for the frequency robustness of long-lasting effects. Our results presented in Fig 8 also suggest that varying the size and number of separately stimulated subpopulations may substantially improve decoupling effects. The PMCS pattern in Fig 8A’ corresponds to a CRS pattern with fixed sequence where one subpopulation is twice as large as the other one. We found that such stimulation leads to robust decoupling in a large portion of the considered interval of stimulation frequencies for single-pulse stimuli and short burst stimuli with high intraburst frequencies, or a variety of burst stimuli with low intraburst frequencies (Fig 8B’ and 8C’). In contrast, if stimulation is delivered to three equally sized subpopulations, the frequency range with pronounced decoupling effects shrinks substantially (Fig 8B” and 8C”). In terms of a controlled modulation of the synaptic network structure, previous theoretical studies found that CRS RVS causes different dynamics of intra- and inter-population synapses. However, all inter-population synapses showed on average the same dynamics [61]. In contrast, CRS with fixed sequence (Δα k = 1/M) leads to different phase lags between separately stimulated subpopulations and potentially allows for inducing a larger variety of network structures. [END] --- [1] Url: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1010568 Published and (C) by PLOS One Content appears here under this condition or license: Creative Commons - Attribution BY 4.0. via Magical.Fish Gopher News Feeds: gopher://magical.fish/1/feeds/news/plosone/