Encoding information into autonomously bursting neural network with pairs of time-delayed pulses

Biological neural networks with many plastic synaptic connections can store external input information in the map of synaptic weights as a form of unsupervised learning. However, the same neural network often produces dramatic reverberating events in which many neurons fire almost simultaneously – a phenomenon coined as ‘population burst.’ The autonomous bursting activity is a consequence of the delicate balance between recurrent excitation and self-inhibition; as such, any periodic sequences of burst-generating stimuli delivered even at a low frequency (~1 Hz) can easily suppress the entire network connectivity. Here we demonstrate that ‘Δt paired-pulse stimulation’, can be a novel way for encoding spatially-distributed high-frequency (~10 Hz) information into such a system without causing a complete suppression. The encoded memory can be probed simply by delivering multiple probing pulses and then estimating the precision of the arrival times of the subsequent evoked recurrent bursts.

the two groups 29 . It was shown that the Δt stimulation protocol was very effective in changing the 'overall excitability' of neuronal cell cultures, which qualitatively represents the likelihood of a PB generation. Incidently, a similar stimulation protocol termed as "desynchronizing delayed feedback stimulation" is found to be very effective in the treatment of several neurological disorders like Parkinson's disease (see 30,31 and the references therein). In this paper, we explain how the particular stimulation protocol could have worked for the system and discuss its novelty in connection with a synaptic memory formation in a PB supporting neuronal system. A detailed explanation is given based on computer simulations of the well-known Izhikevich neuronal network model 32 , which is well-known for its population bursting dynamics. A few years ago, we demonstrated that the Izhikevich model supports the phenomenon of coherence resonance when it is driven by different levels of noise: in other words, there is an optimal level of noise at which the variation of the interburst interval becomes minimal 13 . In this study, we use the same model for a very different purpose, specifically, to explore the role of Δt paired-pulse stimulation as for rendering a meaningful change in the network synaptic morphology and the related PB dynamics. Here, the interplay among the value of Δt, the duration of PB, and the time relevant for the spike-timing-dependent-plasticity (STDP) mechanism, which is built into the model system, is found to be essential.

Mathematical model of pB generating neuronal network
So far, several different mathematical models have been introduced to generate spontaneous neuronal PB oscillations [32][33][34][35] . In this study we use the Izhikevich neural network model 32,36 , since it is easy to simulate on a personal computer to generate some realistic PB dynamics 13 . More importantly, the model incorporates a STDP rulean essential feature for synaptic memory formation and removal. Although the details are rather different, all previous mathematical models of PB generating networks include two competing mean-field factors-activation and inhibition. For the model proposed by Augustin et al. 35 , the two factors are recurrent synaptic excitation and neuronal adaptation (current), respectively. Similarly, the Izhikevich model 32,36 includes recurrent excitatory (and fewer inhibitory) synaptic connections, while burst inhibition is indirectly implemented via the asymmetric STDP curve that is biased toward the suppression side over the potentiation side. Specifically, our model is comprised of 200 (160 excitatory and 40 inhibitory) integrate-fire neurons coupled via a STDP function as described in Fig. 1b. The state of individual neuron evolves according to this set of equations: As the value of v crosses a prescribed threshold value (30.0 mV), the two variables are made to reset: In the above model, the variables, v and u, are the membrane potential and the level of recovery following a reset (i.e., spike firing), respectively. I is the sum of synaptic current inputs (I syn ) from other neurons, an external current (I ext ) and a residual noisy current (I η ). I η = 16 pA was given to only one randomly chosen neuron at every 1 ms. The variable external current was used only for the Δt training or for inducing recurrent PBs for evaluating the time precision of recurrent PBs; otherwise, it was set to be zero. Throughout this paper, the following parameter values are used: v r = −82.7 mV, v t = −42.3 mV, C = 1 nF, k = 0.04 pA/mV 2 , b = 0.2 nS, c = −65 mV, and a = 0.02 (0.1) ms −1 and d = 8 (2) pA for excitatory (inhibitory) neurons.
The number of synaptic connections (to postsynaptic neurons) for each neuron is chosen to be 60. The inhibitory neurons are projected only to excitatory neurons with a fixed synaptic weight of −5 and a fixed conduction delay of 1 ms. On the other hand, the excitatory neurons project to inhibitory as well as excitatory neurons with time-varying (0~10) synaptic weights W ij , and with randomly (but uniformly) selected conduction time delays τ ij (1~20 ms), which are what we have typically measured for the inter-channel conduction delays within cultured networks of neurons. The synaptic current to the ith neuron is: in which t j m is the time of mth spiking of jth neuron. The weight W ij (of a synapse connecting a presynaptic neuron j to a postsynaptic neuron i) of excitatory neurons is updated every second based on the STDP function depicted in Fig. 1b: The change in W ij (ΔW ij ) depends on the temporal distance Δt ij between the time t j of presynaptic spike and the time t i of postsynaptic spike. We define the sequence of presynaptic spiking times to be t j m and that of postsynaptic spiking times to be t i n , where m and n = 1, 2, 3… are temporal indices for the action potentials in sequence. Incorporating the conduction time delays τ ij , the actual time difference between the mth spike of presynaptic neuron (j) and the nth spike of postsynaptic neuron (i) is   Throughout this paper, we used the parameter values of A + = 0.1, A − = 0.12, and τ = 20 ms. The choice of parameter A − > A + suppresses the potential over-excitation of the network. And, as we will explain shortly, the conduction time delay also plays a delicate role as for bringing inhibition to a bursting network.
Our numerical simulation was carried out by customizing the c ++ code, which was originally published in 32 .

Results
Network bursting activity and its suppression by a periodic sequence of single-pulse stimuli. The raster plot in the top-left frame of Fig. 1c shows a typical PB event series (of 3 seconds) that is autonomously generated by our model network. The histogram of the inter burst intervals of spontaneously generated PBs in general shows a rather broad distribution about 300 ms as shown in Fig. 1d. In fact, earlier we suggested that the time series can be viewed as a very noisy PB oscillation 13 , which is formed by two competing factors: inhibition via the depression-favored STDP rule and recurrent synaptic excitation. During each individual PB event, the cell-to-cell synaptic connections are to be weakened in general. There are two different elements in this weakening process. First of all, as mentioned previously the STDP function itself is biased towards the inhibition side: In the original model proposed by Izhikevich this bias was intentionally introduced to avoid network over-excitation 32 . Second, there is a conduction delay τ ij (which is randomly chosen from 1 to 20 ms) of the presynaptic firing of ith neuron towards its jth postsynaptic neurons. With the delays, the synaptic connection weights W ij are more likely to decrease than increase as each PB event comes to an end. This suppression effect is easy to understand if we imagine an unusual PB in which all cells fire simultaneously. For such a case, all synaptic weights will decrease because all post-synaptic neurons fire effectively before their presynaptic partners do (due to conduction delays). For autonomously generated PBs, cells do not fire simultaneously but their AP spikes are dispersed as shown in the left frame of Fig. 1e, more or less in a form of Gaussian function. Even these Gaussian PBs will cause an overall synaptic suppression for the same reasons. Then, the intermissions between any two successive PBs provide the time necessary for the recurrent synaptic excitation to regain its excitability. During these burst-free intermissions, single AP spikes prevail and the system recovers its excitability with the outnumbering excitatory synaptic connections over the inhibitory connections. The balance between the suppression and the recovery process sustains the spontaneous PB dynamics.
An extrinsic, stimulation pulse (for example, 60 pA amplitude, 1 ms duration given to a randomly selected group of 20 neurons) evokes an immediate PB, during which almost all neurons fire simultaneously. Each evoked PB has a firing rate profile which is quite skewed toward the onset of excitation. Then, the skewness becomes more conspicuous as the stimulation repeats as shown in the right frame of Fig. 1e. The 'front heavy' spike firing rate profile of the evoked PB accelerates the overall suppression of the synaptic connections and breaks the delicate balance between the suppression and the natural recovery process. After being exposed to a long series of external stimuli, only a small set of neurons that receive the stimulation directly and their immediate neighbors would be affected by the stimulation. By then, the perturbation cannot be relayed to the whole neuronal population. The increasing overall inhibition under a repeating stimulation is rather evident in three raster plots of Fig. 1c, in which recurrent burst activity following the perturbation (red vertical line) gets weakened and dies out gradually. In other words, the system is overwhelmed even by a sparse (3 s interval) pacing given only to 5% of the total neuronal population.
Novelty of Δt paired-pulse stimulation. In order to process and encode complex input information (eg., high-frequency spatiotemporal sensory signal), a neural network needs to utilize its whole network connectivity and be away from the network suppression like the one discussed in Fig. 1. For the PB generating Izhikevich neuronal model network, a high frequency (eg. 10~50 Hz) input information can be encoded without causing a complete network suppression, for example, as in the following Δt stimulation protocol, in which two selected groups of neurons are stimulated by a pair of pulses with a relative (small) time delay of Δt (see Fig. 2a, right frame) with a typical range of 50~120 ms.
Essentially, the Δt paired-pulse stimulation protocol reinforces the STDP rule simultaneously to two subpopulations of a large coupled neuronal network having many random recurrent connections. Having two evoked PBs in a short duration of time can result in a rather unexpected, interesting consequence. Two PBs together back-to-back with a small delay can make the system much less suppressed (when compared with the case of a single evoked PB) as illustrated in Fig. 2b: Almost all APs evoked by the first pulse are advanced in time than the APs generated by the second pulse for each paired-pulse stimulation; thus, there would be sufficiently many pre-and post-synaptic pair events that yield a synaptic potentiation over depression. That is so, even when the conduction time delays of 0~20 ms are taken into account. Now, the synaptic potentiation can compete with the depression to find a balance. Therefore, even with the increasing number of Δt stimuli each of the AP firing rate profiles of the two paired PBs does not go into a depression but gradually stabilizes into a non-trivial steady-state (the solid purple lines in Fig. 2b) which is not much different from that of the spontaneously generated burst shown in the first frame of Fig. 1e. In other words, the Δt stimulation (a simple spatiotemporal input pattern having a high-frequency component) could be successfully encoded into the whole network. After all, Δt paired-pulses can be far less suppressing than single pulses.
The purpose of the Δt stimulation is to train a PB generating network without completely suppressing its network connectivity; for that matter, a sufficient temporal overlap between the first evoked PB and the following second PB is crucial and we discuss that with Fig. 3. The heat maps in the top row of Fig. 3 quantify the AP firing activity while the system is externally driven by a periodic (period 3 s) sequence of Δt paired-pulse stimuli. For The significance of Δt in the paired-pulse stimulation protocol. The consequence of the perturbation can be quite different with a larger value of time delay between two stimulating pulses, for example, Δt = 120 ms (see Fig. 3b): Now, the system becomes severely depressed as the profiles of two evoked PBs get quickly separated in time and their APs loose STDP mediated interactions. The increasing separation is quite clear in the color activity map of Fig. 3b. Notably, for this case, 〈W ij 〉 decreases while R stim rises sharply as the number of stimuli increases. A large value of R stim means that there are too many simultaneous AP firings that would cause an overall synaptic depression. Again, for such a case the conduction delays (d = τ ij ) are a major factor determining the (suppressing) action of the STDP function shown in Fig. 1b. Eventually, all the synaptic weights, except for those of output connections emanating from the neurons that directly receive the extrinsic perturbation, go to zero. Figure 4 quantifies the evolution of probability distribution function of W ij for such a case. We view that  〈 〉 W 2 ij represents a nearly complete network suppression. In fact, any paired-pulse stimulation with  ∆t 120 ms is observed to undergo a complete suppression. If the value of Δt gets too small (eg. 20 ms), the two evoked PBs practically overlap each other to form a (roughly twice) densely packed single-hump PB that would give a more severe depression to the network connectivity than a single pulse stimulation would induce. After all, there should be an appropriate range of Δt values, with which one can train the network without causing a complete depression. For the selected set of system parameter values that we have used in this work, the 'non-depressive Δt range' is approximately 30~120 ms: Only for this range 〈W ij 〉 is greater than 4 and it decreases sharply to ~2 below 30 ms and beyond 120 ms. The degree of stimulation-driven depression is, of course, a function of the number of stimuli received by the system as shown in Fig. 3b: The depression sets in gradually and saturates around the 700th stimulation. More importantly, the degree of depression also depends on the group of selected neurons that receive the stimulation. For example, the result shown in Fig. 3c is obtained with the same Δt = 100 ms that was used for Fig. 3a, but with a different group of neurons being stimulated. For the case of Fig. 3c, the system gets depressed rather slowly at the beginning but accelerates quickly around 1000th stimulation to a complete depression. After all, these features (the dependence on the value of Δt, the group of neurons receiving the stimulation, and the number of stimuli) are exactly what one would want for encoding spatiotemporal information into the network.

Accessing the excitability of Δt-trained networks by a probing pulse stimulation. A useful
information about the dynamic state of a Δt-trained network can be extracted by delivering a probing pulse. A probing stimulation is a single pulse (60 pA amplitude, 1 ms duration) stimulation simultaneously given to a set of (randomly) selected 20 neurons. The raster plot of Fig. 5a (left frame) shows a typical response of the system to a probing pulse stimulation, which immediately follows a Δt training session. The almost periodic sequence of recurrent PBs (with a period ~0.5 s) can be viewed as a reminiscence of the noisy subthreshold PB oscillator that was discussed in ref. 13 . The recurrent PBs shown in Fig. 5a (left frame) are produced by a network which is trained with 300 Δt = 60 ms paired pulses, while the ones shown in Fig. 5a (right frame) are for 300 Δt = 100 ms paired pulses. Clearly, the recurrent PBs in the left frame of Fig. 5a are more dispersed than those in the right frame of Fig. 5a.
Importantly, the timing precision RBTP ≡ <T>/σ(T) (or RBTP = 1/coefficient of variation of T) of the recurrent bursts following a probing pulse stimulation depends on the value of Δt. Here, <T> and σ(T) refer to the mean and standard deviation of the time intervals between the probing pulse stimulation (marked by a red arrow  in Fig. 5a) and the first recurrent burst (marked by blue dot in Fig. 5a) measured over 100 different trials given at an interval of 10 seconds (while the values of W ij are fixed after the completion of each Δt training). Figure 5b plots RBTP vs. Δt for ten different initial network states (i. e., different initial W ij ). Clearly, the value of RBTP depends not only on the value of Δt but also on the initial network connectivity structure, the one just before the training pulses are begun to be delivered. In other words, the encoded memory state (accessed by the measure of RBTP) depends on the past history of the network before the training process starts. The dynamic network evolves autonomously in time, and practically the high dimensional system will never be in the same dynamics state when the stimulation is being delivered. In other words, a new input information gets encoded on the top of what is currently being stored in the network. Furthermore, as expected the value of RBTP also depends on the specific composition of the two groups of neurons receiving the training stimuli directly as well as on the value of Δt used for a priori training session (see Fig. 5c). There seems to be no simple functional relationship between the value of RBTP and Δt for the range of Δt, for which the system does not go into a complete inhibition.

Recurrent burst timing precision (RBTP) versus other measures of the network. The level of
individual neural spiking activity can be considered as a measure of the overall excitability that is relevant for the emergence of a population burst. The mean synaptic strength 〈W ij 〉 is a major factor affecting the spiking activity, and on average it is almost a linear function of RBTP: Shown in Fig. 6a are noisy, yet, almost linear functional relationships between 〈W ij 〉 and RBTP for 10 different initial states (marked by different colors) and 8 different Δt values (dots connected by a line). RBTP also depends on the uniformity of W ij : The larger the coefficient of variation C var of W ij is, the smaller RBTP is (see Fig. 6b). In other words, the more homogeneous the network connection weights are, the more reliable and accurate the recurrent PBs are when they are induced by a probing pulse stimulation. Also, RBTP is a linearly increasing function of the actual amount of synaptic transmission (traffic), which is defined by the sum of all weighted synaptic inputs triggered by APs during the total Δt training period of 900 s (see Fig. 6c). The total synaptic transmission can be further divided into three different categories, namely, the transmission from excitatory neurons to another excitatory neurons, that from excitatory neurons to inhibitory neurons, and that from inhibitory neurons to excitatory neurons. In all three cases, RBTP is found to be a linearly increasing function of synaptic transmission (see Fig. 6d). Another important factor that characterizes the connectivity of a functional network is the closeness centrality C C , which measures the degree of how near each neuron is to all other neurons or how well neural information spreads on the network. C c is computed by the following procedure. First of all, is computed for all connected pairs from ith neuron to jth neuron (only W ij > 0 is considered). Then, the weighted distance between mth neuron and nth neuron (dist mn ) is defined to be the shortest path P mn connecting the mth neuron to the nth neuron: where N m is the number of neighbors of the neuron indexed by m. The program CytoNCA 37 is used to calculate C C . Figure 6d plots RBTP vs. 〈C C 〉. As expected a well-connected network with a larger value of 〈C C 〉 gives a better precision in the timing of the recurrent burst. Thus, the 'overall excitability' of the PB generating neural system surely depends on the degree of network connectivity as well as the network homogeneity and 〈W ij 〉.
The value of RBTP also strongly depends on the location where within the network the probing pulse stimulation is delivered, as the network connection morphology matters for the value of RBTP. Naturally, the probing pulse delivered to the very group of neurons (labeled by group index #0 in Fig. 7a), which was used for the Δt training, yields the largest RBTP over 10 other trials made with other groups of neurons selected for the probing stimulation (see Fig. 7a). As explained earlier, during the time course of a Δt training, the efferent synaptic weights of the neurons that receive the training pulses would be generally enhanced (from 〈W ij 〉 group before training to 〈W ij 〉 group after training ), and Fig. 7b well confirms that: The mean synaptic weight of those neurons used for the Δt training has significantly increased after the training process. After all, the mean synaptic weight is the largest (~5) for the group of neurons that received the training pulses directly and RBTP is the largest when the same group of cells are stimulated by a probing pulse: The case is marked by a red dot and index #0 in Fig. 7c. In other words, with an appropriate Δt training some strong pathways are established connecting the group of neurons, which receive the Δt stimuli directly, to other neighboring neurons. The input information is impregnated into the network weight morphology forming a memory.

summary and Discussion
Essentially all STDP rules, developed so far, are for a pair of neurons: The rules consider the time difference(s) between two (or more) spikes, for example, one for the pre-synaptic neuron and the other for the post-synaptic neuron, and modifies the strength of the synapse connecting them as a nonlinear function of the time difference. Understanding the outcome of this rule applied to not just a pair but many neurons forming a well-connected network is far from trivial due to the intrinsically polysynaptic nature of the circuits having plastic and recurrent connections. Interesting issues can emerge 91) if such a system autonomously generate population bursts with a more or less well-defined time interval and 92) if it is perturbed by a periodic sequence of pulses. Self-sustained PB activity is sustained with a delicate balance of excitation and inhibition, as discussed in ref. 38 . When the PB generating network is stimulated at a high frequency, it will become easily depressed and fragmented. On the other hand, a temporally sparse stimulation will neither affect the autonomous network bursting activity in any meaningful way nor be able to form any input-specific synaptic structural adaptation in the network. The Δt paired-pulse stimulation protocol that we developed and implemented both in the earlier experiments 29 and the current model study is in a way a compromise between these two extremes and works in a regime causing significant changes to the network connectivity, yet, avoiding a complete depression. Most importantly, the Δt stimulation protocol includes a high-frequency component (i.e., the small value of Δt), which is relevant for the STDP rule, thereby enforcing a non-trivial interplay between the network plasticity and the external information.
In addition to the time delay Δt, which is associated the extrinsic paired-pulse perturbation, our model has also incorporated uniformly distributed conduction time delays τ ij : They were added to the time difference of spiking of the coupled neurons when the STDP function was evaluated. The time delays τ ij were, more precisely speaking, axonal conduction time delays as versus to dendritic time delays. According to a series of recent phase model (for weakly pulse-coupled oscillators) studies by Madadi Asl et al. 39,40 , the axonal and dendritic conduction delays both could be delicate parameters so that, for example, their different balances could bring quite a few Figure 7. The dependence of RBTP on the groups of neurons that receive a probing pulse (a) RBTP vs. group used for probing pulse stimulation (#0 red is for the same group of neurons that were used for Δt training beforehand, all others are just randomly selected groups of neurons), (b) The overall increase in the mean synaptic weight of neurons which were chosen to receive the Δt = 100 ms training pulses, and (c) RBTP vs. 〈W ij 〉 group after training . Note in (a) that RBTPs for 10 randomly chosen groups of neurons all are not as good as that of the first recurrent burst generated by a probing pulse delivered to the group of neurons that were used for the prior Δt training. different network connectivity patterns, including recurrent bidirectional couplings. In addition, their recurrent network model could also support multistability that were rendered visible with different initial distribution of the synaptic strengths. Knowing that in our model system the value of RBTP (after a Δt stimulation) has fluctuated significantly depending on the initial landscapes of W ij as shown in Fig. 5b, it is tempting to speculate that the fluctuation could be originating from the multistability of our (unperturbed) network model. So, in the future it will be interesting to check the possibility of multistability and explore the role of conduction time delays (eg. non-uniform distributions of axonal and dendritic time delays and their different balances) in a systematic way. Earlier, we claimed that the recurrent burst timing precision RBTP is a good readout measure of the overall excitability relevant for the emergence of a population burst 29 . In this work, we demonstrated that not only the mean synaptic weight 〈W ij 〉 but also the network connectivity, as measured by C var , and the homogeneity, as measured by 〈Traffic〉 and 〈C C 〉, are an important factor as for determining the overall excitability. In addition, we showed that the value of RBTP also relies on the morphology of the trained network: The detailed landscape of synaptic weights matters as for determining RBTP, as the probing pulse is given only to a small subset of neurons. After all, the Δt paired-pulse stimulation conveys an externally imposed spatiotemporal information to the network and forms a memory in the plastic landscape of synaptic weights. At this point, we should indicate that synaptic plasticity is often explored as a form of unsupervised adaption of neural circuits learning the structure of complex sensory inputs 41 . The evaluation of RBTP following the probing pulse stimulations is in a way for retrieving the Δt stimulation-encoded information. However, there is no simple, a priori set, functional relationship between RBTP and the value of Δt, as the initial state of the given network, the number of stimulations being given, and the subgroups of cells receiving the stimulations are also an important factor.

Conclusion
The modulation of RBTP in cultured networks of cortical neurons by Δt stimulation reported earlier in ref. 29 was successfully reproduced in computer simulations of a simple neuronal network model in which synaptic weights evolve according to a simple STDP rule as described in 32 . Subsequently, with a series of careful analyses we elucidated the significance and novelty of the Δt paired-pulse stimulation protocol as for an effective way of inducing input-specific structural changes as well as for modifying the overall bursting excitability of neuronal network. Finally, we indicate that the Δt stimulation that we developed may represent one of the most simplest forms of spatiotemporally structured sensory inputs that have a high frequency component. As such, in the future we will explore how the same burst-generating neural network would adapt its synaptic weights in response to more complex spatiotemporally patterned stimuli. Moreover, we can explore the same issues with a more sophisticated STDP rule that includes triple-spike, or even quadruple-spike, nonlinear interaction as in [42][43][44] for our future work.

Data Availability
All numerical data used in this work will be available upon request to the corresponding author.