Communication through Resonance in Spiking Neuronal Networks

The cortex processes stimuli through a distributed network of specialized brain areas. This processing requires mechanisms that can route neuronal activity across weakly connected cortical regions. Routing models proposed thus far are either limited to propagation of spiking activity across strongly connected networks or require distinct mechanisms that create local oscillations and establish their coherence between distant cortical areas. Here, we propose a novel mechanism which explains how synchronous spiking activity propagates across weakly connected brain areas supported by oscillations. In our model, oscillatory activity unleashes network resonance that amplifies feeble synchronous signals and promotes their propagation along weak connections (“communication through resonance”). The emergence of coherent oscillations is a natural consequence of synchronous activity propagation and therefore the assumption of different mechanisms that create oscillations and provide coherence is not necessary. Moreover, the phase-locking of oscillations is a side effect of communication rather than its requirement. Finally, we show how the state of ongoing activity could affect the communication through resonance and propose that modulations of the ongoing activity state could influence information processing in distributed cortical networks.


Introduction
The brain processes sensory stimuli by an organized flow of neuronal activity across a distributed network of specialized cortical areas. This flow requires mechanisms that route neuronal signals from one cortical area to another. However, the exact nature of this routing process remains poorly understood.
Experimental studies suggest that synchronization of spiking activity may play a pivotal role in the flow of neuronal activity, as synchronous neuronal firing can effectively drive downstream neurons [1][2][3][4]. To date, our understanding of synchrony-based neuronal routing has been dominated by two models which attribute the origin of synchrony to dissimilar mechanisms.
According to the first model, synchronous spiking activity is both created and routed through dense and/or strong convergentdivergent connections between subsequent layers of feedforward networks (FFNs). In this scenario, these connections are a source for shared and correlated input that provides sufficient synchronization for spiking activity to propagate across the FFN [5][6][7][8][9][10].
However, the requirements of either strong synapses or high connection probability pose serious constraints on the biological plausibility of these FFNs in the cortex, in which connectivity is in general sparse [11] and synapses are weak [3,9,12]. Even though, the sparse cortical connectivity could in theory host a large number of sparsely and weakly connected (diluted) FFNs, they would fail to generate enough synchronization to ensure propagation of spiking activity [7,13,14].
The second model suggests that population oscillations could soften the requirement of strong/dense connectivity by enhancing synchronization and neuronal excitability during the excitable phase of the oscillation [15,16]. A key requirement for this propagation mode is that oscillations, which are generated locally due to interactions between excitatory and inhibitory neurons, must maintain a consistent phase relationship (coherence) between the communicating networks (''communication through coher-ence''; [15,17,18]). However, the mechanisms underlying the generation and maintenance of such coherent oscillations between distant brain areas have remained elusive despite a number of theoretical proposals [19].
Here, we propose a novel mechanism by which oscillatory activity exploits the presence of resonance frequencies in networks of excitatory and inhibitory neurons (EI) to promote the propagation of synchronous activity across diluted FFNs (''communication through resonance''). The role of such network resonance is to amplify weak signals that would otherwise fail to propagate. According to our model, coherent oscillations emerge in the network during slow propagation of synchrony, while at the same time synchrony needs these oscillations to be propagated. Thus, spreading synchrony both generates oscillations and renders them coherent across different processing stages. This abolishes the requirement for separate mechanisms providing the local generation of oscillations and establishing their long-range coherence. Moreover, coherence between oscillations may be viewed as a consequence of propagation instead of being instrumental to establish communication through synchrony. Our results also suggest that the emergence of coherent oscillations is influenced by the dynamical state of the ongoing activity. We propose that changes in the ongoing activity state can have an influence on cortical processing by altering the communication between different brain areas.

Network model
The network models were multi-layered FFNs. Each layer consisted of two recurrently connected homogeneous neuronal populations. In Figure 1, Figure 2 and Figure 3 we used 2,000 excitatory (E) and 500 inhibitory (I) neurons. For the rest of the figures, we reduced the number of E neurons to 1,000 while keeping the number of interlayer projecting neurons fixed to 300. This reduction, which was done in order to improve simulation efficiency, did not affect the results in any qualitative manner.
The connectivity within each layer was random with the following connection probabilities: E EE~0 :05 and E IIẼ EI~EIE~0 :1, where E XY denotes the probability of connection from a neuron in population X to a neuron in the population Y . Connections between layers were strictly feedforward and excitatory, and restricted to a sub-population of 300 randomlychosen E neurons (in the rest of the paper referred to as P) in every layer. The interlayer connectivity was sparse with probability E PiPiz1~EPP~0 :1 (cf. Table 1).

Neuron model
Neurons were modeled as leaky integrate-and-fire neurons, with the following membrane potential sub-threshold dynamics: where V m is the neuron's membrane potential, I syn is the total synaptic input current, C m and G leak are the membrane capacitance and leak conductance respectively. When the V m reached a fixed threshold V th~{ 54 mV a spike was emitted and the membrane potential was reset to V reset~{ 70 mV: After the reset, the neuron's membrane potential remained clamped to V reset during a time period t ref~2 ms, mimicking the period of absolute refractoriness. All other parameters are detailed in Table 2.

Synapse model
Synaptic inputs consisted of transient conductance changes: where E syn is the synapse reversal potential. Conductance changes were modeled using exponential functions with t E~5 ms and t I~1 0ms. Synaptic delays were set to d EE~dII~2 ms, d EI~dIE~5 ms and d PiPiz1~dPP~5 ms in Figure 4a, Figure 5b-c and Figure  S4a. In the rest of the figures, delays were set to d EE~dII~1 ms, d EI~dIE~2 :5 ms and d PP~5 ms. Longer delays produced a stronger and more reliable propagation and therefore were chosen to illustrate the propagation across layers in Figure 4a. The choice of delays influenced the resonance properties of the network [20]. However, the general principle remained unaffected. Other parameters are detailed in Table 3.

External input
Each neuron was driven by 1,000 independent Poisson excitatory spike trains with an average rate of 1 Hz each (i.e., a total average input rate of 1 kHz), which mimicked uncorrelated background inputs coming from other brain areas. In Figure 6, E neurons received this external drive (referred to as E drive) with larger rates than 1 kHz as indicated in the figure.
The synchronous stimuli consisted of periodic trains of synchronous spikes (pulse packets) with different frequencies. Only P 1 neurons received these additional spikes. The individual pulse packets consisted of a fixed number (a) of spikes per neuron, distributed randomly around an arrival time t n . The time of each individual spike was drawn independently from a Gaussian probability distribution centered around t n and with s.d. (s). In Figure 2e, Figure 3b-c, Figure 4a, a~30 spikes and s~0 ms (i.e, perfectly synchronous). In the remaining cases a~20 spikes and s~3 ms.

Author Summary
The cortex is a highly modular structure with a large number of functionally specialized areas that communicate with each other through long-range cortical connections. It is has been suggested that communication between spiking neuronal networks (SNNs) requires synchronization of spiking activity which is either provided by the flow of neuronal activity across divergent/convergent connections, as suggested by computational models of SNNs, or by local oscillations in the gamma frequency band . However, such communication requires unphysiologically dense/strong connectivity, and the mechanisms required to synchronize separated local oscillators remain poorly understood. Here, we present a novel mechanism that alleviates these shortcomings and enables the propagation synchrony across weakly connected SNNs by locally amplifying feeble synchronization through resonance that naturally occurs in oscillating networks of excitatory and inhibitory neurons. We show that oscillatory stimuli at the network resonance frequencies generate a slowly propagating oscillation that is synchronized across the distributed networks. Moreover, communication with such oscillations depends on the dynamical state of the background activity in the SNN. Our results suggest that the emergence of synchronized oscillations can be viewed as a consequence of spiking activity propagation in weakly connected networks that is supported by resonance and modulated by the dynamics of the ongoing activity.
When the stimulus was a periodic train of pulse packets, we set the frequency of stimulation by adjusting the period (T) between arrival times t n (i.e., the center of the Gaussian p.d.f.). When sw0, the spikes were spread around t n , as indicated above, and therefore the time distance between the last spike from a given pulse packet and the first spike from the next was always variable for the same input frequency. The smallest interval that was used between arrival times was 10 ms (100 Hz) and the largest 100 ms (10 Hz). Additionally, in Figure 3b we used 1 Hz stimulation.
In simulations where the arrival times were jittered, the size of the jitter was drawn from a uniform distribution centered on the arrival time t n . The extent of the jitter window was chosen to be a function of the interval t n +T=f , where f~8, 4 or 2 in order to make the effect comparable across different frequencies.

Data analysis
To compute the auto-covariance functions A(t) (inset in Figure 2c and Figure 3c bottom right; only positive time lags are shown), time was divided into bins of Dt~5 ms and the population spike trains were transformed into spike count vectors y a (t), where a~E, I, P denotes the population. The autocovariance functions were then computed as follows: where t~({D, {DzDt, . . . , D), t~(D, DzDt, . . . , T s {D), D~150 ms, Dt~5 ms, n a indicates the population mean firing rate and the superscript b~o, a denotes ongoing (computed from a single T s~1 00 s simulation in absence of pulse packet stimulation) and activated (computed from T s~1 0 s of activity during stimulation starting 5 s after the stimulus onset and averaged across 20 trials), respectively. We used the population Fano factor (pFF) to classify the population spiking activity states as synchronous or asynchronous (dashed line in Figure 6a). We used the central value of normalized by the mean population firing rate: The signal-to-noise ratio (SNR) in Figure 6d was computed as follows: indicates the variance of the spiking activity of P 10 neurons as indicated above.
Pairwise correlations were computed using the Pearson correlation coefficient between the spike count vectors of pairs of neurons (y i (t) and y j (t)). where: Cov(y i ,y j )~S(y i (t){n i )(y j (t){n j )T~Sy i (t)y j (t)T{n i n j and S : T indicates time average and vectors y i (t) and y j (t) were computed using a time window of D t~200 ms. We used 10,000 pairs to compute the distributions shown in Figure 2c and where, b~a in Figure 4b, Figure S4a and Figure S5a and b~o in Figure 6b indicating the corresponding value of T s (cf. description of the auto-covariance function above).

Simulation tools
Network simulations were performed using the simulator NEST, interfaced with PyNest [22,23]. The differential equations were integrated using forth order Runga-Kutta with a time step of 0.1 ms. Simulation data was analyzed using the Python scientific libraries: SciPy and NumPy. The visualization of the results was done using the library Matplotlib [24]. The code to reproduce several results presented in this work ( Figure 1b, Figure 3a, Figure 4a, Figure 5b and Figure S6a) is available at https:// github.com/AlexBujan/ctr. Other results can be reproduced by modifying that code.

Network model: Diluted FFNs
We studied the propagation of synchronous spiking activity across diluted FFNs with sparse interlayer connectivity. In this model, each layer represented a small neocortical network with 2,000 excitatory (E) and 500 inhibitory (I) neurons. The connectivity within each layer was sparse and random. The connections between layers, which modeled long-range projections between different cortical networks, were strictly feedforward and excitatory. These interlayer projections were restricted to a  sub-population of 300 E neurons which we refer to as projecting neurons or P neurons throughout the manuscript (P i refers to all the projecting neurons in a layer with the subscript i indicating the position of the layer in the FFN; cf. Figure 1a). Interlayer connection probability E PP~0 :1 and hence, each P i neuron received, on average, N P E PP~3 0 connections from the previous layer (P i{1 ; cf. Table 1).
All layers were driven by external Poisson input spike trains and the synaptic weights were adjusted (cf. Table 3) to bring the network into an asynchronous-irregular (AI) activity regime [25,26], consistent with the statistics of cortical activity in awake behaving animals [27][28][29][30]. The mean firing rate of individual excitatory neurons showed a heavy tailed distribution with a mean of *1 Hz (+1:1 Hz; + s.d. across the population). The mean coefficient of variation of the inter-spike interval distribution (CV ISI ) was *0:95 (+0:19) and the distribution of pairwise correlations was centered around zero with a mean of * 1|10 {3 (+0:04) (cf. Figure 1b and Figure 2a-c). The activity of the I population was also irregular and asynchronous although with slightly higher mean firing rates (2+1:8 Hz). These results were computed from a single simulation of 100 s duration.
To study the propagation of synchrony, we stimulated all P neurons in the first layer (P 1 ) with synchronous events or pulse packets (cf. Methods). The synaptic strength of these input synapses was equivalent to the other EE synapses in the FFN (cf . Table 3).
First, we checked that the connectivity between layers was indeed too weak to support the propagation of single synchronous events. To this end, we generated an amplitude transfer map which we used to estimate the change in amplitude undergone by pulse packets as they travel across the FFN. This map, shown in Figure 2d, was generated using the ongoing membrane potential distribution (black trace) and depolarization transfer function (dark gray solid trace) of the P population. The measured membrane potential distribution (computed from 100 s of ongoing activity) is shown as the cumulative density function (c.d.f.) of the distance to threshold ({54 mV). When represented as such, the probability of being at a certain distance from threshold can be interpreted as the fraction of cells (here named k i , where i indicates layer index) that will spike if a depolarization (''jump'') equivalent to such distance is applied to all P cells. The membrane potential transfer function was calculated by measuring the averaged maximum depolarization across P neurons induced by a pulse of perfectly synchronous spikes (s~0) with different amplitudes a. The mapping between the two curves can be done by knowing the relationship between the activation level of the ith layer k i and the amplitude of the pulse packet received by the subsequent layer a iz1 which in this case is as follows: a iz1~ki N P E PP . Knowing this relationship, it is then possible to project a point from one curve to the other, thereby drawing an estimated trajectory of the pulse packet's amplitude across the chain. In the figure, an example of such a trajectory is illustrated with red dots and dotted lines. To make a convincing case, we started the trajectory with a fully activated first layer (k 1~1 ; upper red dot) and followed the pulse packet until it reached a stable point (intersection between the two curves). Such trajectories will always end at an intersection between the curves which in this case (E PP~0 :1) is found only at zero. This shows that any single pulse traveling across this FFN will eventually vanish, regardless of the initial value of a. Similarly, it can be shown that if the connectivity is raised to

Resonance in recurrent EI networks
After a perturbation caused by a synchronous pulse, the network's activity relaxed back to ongoing levels while displaying a stereotypical damped oscillation (Figure 2e). This dynamics, which was observed both at the spiking level (shown as conductances in P neurons in Figure 2e top) and the level of the membrane potential (Figure 2e bottom), indicated that the network had resonance frequencies. The presence of such resonance frequencies suggested that stimulating the network with a periodic train of pulse packets, within a specific frequency range, could induce a large response even for weak stimuli (e.g., pulse packets consisting of a few weakly synchronized spikes).
The existence of resonance behavior in EI networks has already been shown elsewhere [20]. Here, we analyzed the network response to a pulse packet stimulation in order to understand in more detail how resonant dynamics can emerge in these networks. During the transient damped oscillatory response, there was a brief time period of a few milliseconds (indicated approximately as a gray region in Figure 2e) during which P neurons were slightly    more depolarized (higher mean), more synchronous (decreased s.d.) and their inhibitory conductance was reduced. This suggested that the arrival of a second pulse packet inside this brief time window (e.g., around 45 ms after the arrival of the first pulse packet; shown as a green dot in Figure 2e) should result in a larger activation as compared to the first pulse. Conversely, the arrival of a second pulse outside of this window (magenta dot in Figure 2e) would only lead to a similar or even weaker activation.
To confirm this, we stimulated P neurons in an isolated layer with a sequence of 100 periodic pulse packets (identical to the ones described in the previous section; a~30 and s~0) and computed the mean firing rate within 20 ms after the arrival of each synchronous event (which was found to be an appropriate time window to capture the pulse packet induced modulation of the firing rate). We repeated the experiment using three different time intervals T: 35, 45 and 1,000 ms ( Figure 3b) and in each case the results were averaged across 100 trials. Pulse packets separated by T~45 ms, which matched the optimal window described above, resulted in an average spiking activity of 48 Hz (+1:7 Hz, + s.d. across trials; green bar in Figure 3b). By contrast, stimulation with pulse packets separated by T~35 ms could only induce a mean network response of 33:7 Hz (+14:8 Hz). This response was comparable to a stimulation in which pulse packets arrived at an interval of one second, long after the transient response to each individual event had died out (compare magenta and blue bars in Figure 3b). This result confirmed that a train of periodic pulses, with a period adjusted to match the optimal time window, was able to elicit a stronger response as opposed to a single pulse packet. Additionally, the fact that a higher input frequency resulted in a lesser activation suggested that this effect was not merely due to the temporal integration of the individual pulse packets.
To further understand the emergence of resonance in these networks, we analyzed the temporal evolution of the membrane potential distribution (mean and s.d. sampled 1 ms prior to the arrival of each pulse packet) during stimulation with a train of 100 pulse packets separated by 45 ms (Figure 3c). The results were averaged across 100 trials. A brief initial depolarization, caused by the first two pulse packets, was followed by a sustained hyperpolarization in both P and I neurons as more pulse packets were presented. The hyper-polarization reflected that a larger fraction of neurons was refractory (or close to the spike reset potential) due to the increase in firing rate (light gray bars in Figure 3c) and recurrent inhibition. The fact that most neurons were more hyperpolarized seemed to be at odds with the observation that the pulse packets were more effective in driving P neurons. Furthermore, a decrease of the s.d. (Figure 3c inset) indicated that P neurons were overall more synchronized, namely, that the hyper-polarization was shared across the entire population. Essentially, the increased responsiveness was a consequence of the fact that I neurons were effectively refractory at the time of the arrival of pulse packets, as indicated by the progressive reduction in their firing rates (Figure 3c dark gray bars). That is, although P neurons moved farther away from the spiking threshold, they received less inhibition at the time of the arrival of the pulse packets which resulted in stronger activation. This observation hinted to an important role of II connections in the emergence of resonance in these networks.
We investigated the contribution of the II loop to the generation of resonance by conducting simulations in which we progressively reduced the strength of the recurrent inhibitory connections ( Figure S1). We compensated the reduction in II input by adding an additional source of external inhibitory conductance in order to keep the firing rate of the I neurons (measured during the ongoing state) constant across conditions. Our results showed that although the II loop had a substantial effect on the resonance peak's amplitude and frequency, the network still had resonant properties in the absence of an II loop. This indicates that while EI connections are sufficient to create resonance, II dynamics play a facilitating role. In addition to the hyperpolarizing inhibition used in our model, other biologically plausible mechanisms, such as shunting inhibition or gap junctions, could also enhance resonance [31][32][33].
Although the overall activity of an isolated layer became more synchronized, with network oscillations that were locked to the stimulus, the overall activity of E neurons remained fairly irregular (CV ISI *0:95+0:20), and mean pairwise correlations were still relatively low (*0:06+0:11; compare (Figure 3d and Figure 2c). Hence, the activity of E neurons during stimulation was still consistent with biological data, which shows that cortical firing is highly irregular despite the presence of oscillations at the population level as measured by local field potentials [34,35]. Note however that the activity of the P neurons was more regular (they skipped fewer cycles) than the other E neurons. Such a level of regularity in the P population was needed in order to induce oscillations in the post-synaptic layer and was a consequence of the small number of P neurons together with the sparse inter-layer connectivity. Thus, the choice of a larger P population size and/or a higher connection probability could make propagation compatible with a more irregular firing in the projecting population (cf. below).
Additionally, we explored whether our network model operated in a linear regime in which case the tools of linear systems analysis could be applied to further understand the resonance [20]. To this end, we calculated the amplitude of the network's response when stimulated with synchronous pulses for different values of a. Our results indicated that the behavior of the simulated network was generally non-linear showing a saturation of the response amplitude with high a and a progressive shift in the resonance frequency ( Figure S2). However, we also found that within a restricted range of input amplitudes the network's response approached linearity (cf. straight lines in Figure S2e).

Communication of synchrony through network resonance
Next, we addressed the question whether the network resonance-induced amplification of stimulus responses, observed in isolated layers, could be sufficient to enable the transmission of synchrony in diluted FFNs, which did not support the propagation of individual pulse packets.
To this end, we stimulated a 5-layer FFN with three different frequencies, that were analogous to the ones introduced in the previous section (cf. Methods). The amplification, observed when the input frequency matched the resonance frequency of P 1 , proved to be sufficient to induce a successful transmission across the entire FFN (Figure 4a bottom). As expected, when the stimulus had a different frequency from the resonance frequency, or it consisted of a single pulse packet, the synchronous activity did not reach the last layer (Figure 4a top and middle). Since the transmission relies on the network resonance, we refer to this mode of synchronous activity propagation as ''communication through resonance'' (CTR).
After receiving a few input cycles at the resonance frequency, nearly all P 1 neurons started to fire near synchronously every time a new pulse was presented. At this point, even though a large number of synchronous spikes were produced in the first layer, the sparse interlayer connectivity (10 %) reduced this increased activation to a train of weak pulse packets with an average of a~28:21 spikes (+2:65 spikes) and s~6:32 ms (+0:5 ms), which prevented the propagation of synchronous volleys immediately after amplification had taken place in P 1 . Therefore, amplification through resonance was needed at every layer to propagate the activity across the FFN due to the diluted interlayer connectivity.
Next we investigated how the frequency of stimulation affected the propagation of synchrony across a 10-layer diluted FFN. Expectedly, we found a correlation between resonance-induced increase in synchrony in P 1 and the successful communication of synchronous events across the entire FFN (Figure 4b). To quantify the synchrony we calculated the variance of the P population spike train (V a P (v) where v indicates the stimulus frequency; cf. Methods). We then used V a P (v) to construct resonance curves as shown in Figure 4b bottom. A propagation was labeled as successful when V a P10 was significantly increased (wmeanz2|s:d:; white dots in Figure 4b bottom) with respect to the baseline value V o P10 . The spectral analysis of the spiking activity revealed that the increase in power in the last layer was always more pronounced at *20 Hz, which was approximately the resonance frequency of the network (cf. Figure 4b lowermiddle subpanel). Furthermore, CTR was not restricted to the FFN architecture discussed thus far. Our results showed that at least two alternative interlayer connectivity patterns also supported CTR: when receiving neurons were restricted to a specific subpopulation of E neurons but any E neuron could project to the next layer ( Figure S3a); when any E neuron could receive and send projections ( Figure S3b).
However, even when neuronal activity propagated to the last layer (white dots in Figure 4b), V a P10 was significantly lower than in V a P1 (compare red and blue curves in Figure 4b bottom). This result indicated that propagation was occasionally characterized by failures of synchronization of the last layers. Thus, the ratio V a P10 =V a P1 could be used as a proxy for the propagation reliability when activity was observed during long time periods (10 s). Generally, networks that produce a moderate amplification of the signal at the resonance frequencies would be more sensitive to noise fluctuations, which can transiently reduce the degree of synchrony and lead to frequent propagation failures. A larger amplification, which in our model was achieved by introducing longer delays within each layer, lead to a perfectly reliable propagation at the resonance frequencies ( Figure S4a).
For the parameters used here, the range of frequencies that led to a successful propagation approximately spanned from 22 to 26 Hz. The extent of this frequency range can be varied by an appropriate choice of network parameters (cf. Figure S4a; see [20] for a more detailed study on the effect of different parameters on resonance). The effect of different parameters on the resonance profile of the network can be estimated using the network's average response to a single pulse packet stimulation (cf. Figure 2e). When the input frequency is expressed as the time interval between pulse packets V a P1 (T) (dashed black trace in Figure 2e top), the resonance profile can be related to the average network response. Note that V a P1 (T) is again represented as a function of the input frequency V a P1 (v) in Figure 4b (blue trace). As can be seen in Figure 2e, the dominant peak in V a P1 (T) closely matches the trough of the average inhibitory conductance response (G in ; red curve). This suggests that the network's response to a single pulse packet stimulation can predict its resonance curve and thus can be used to understand how different changes in the network parameters may affect the resonance properties of the network. While different network parameters can alter its resonance curve, the activity propagation based on network resonance would remain essentially the same.
For this specific choice of parameters, PS a E1 (v) (Figure 4b top; cf. Methods) revealed that the resonance occurred mainly around two main stimulus frequencies: 23 Hz and 58 Hz (see also Figure  S5b bottom row). Note that similar resonance frequencies were found when neurons were stimulated with a sinusoidally modulated Poisson input, which indicates that the faster resonance frequency can not be explained by the existence of harmonics of the base frequency present in the periodic input pulse train ( Figure S6). Naturally, the smaller resonance frequency precisely matched the time window described in the previous section. The frequency of the second resonance peak can be explained using the network's average response as indicated earlier. To understand this effect, we can consider a simpler stimulus consisting of three pulses the frequency of which is systematically increased with respect to the main resonance frequency (23 Hz). Initially, the rise in frequency will cause the second and third pulses to arrive outside the optimal time window. However, as the frequency is further increased, a frequency will be reached for which the third pulse will fall inside the optimal window giving rise to an increase of the spiking response. Intuitively, this latter frequency should be approximately twice as large as the main resonance frequency, which is inconsistent with our results. This discrepancy can be understood when we notice that the second pulse, although not strong enough to activate P neurons, does accelerate their repolarization, thereby advancing the optimal time window within which the third pulse should arrive. That is, the subthreshold effect of these incommensurate pulses will speed up the network response resulting in the second resonance peak being faster than twice the main resonance frequency.

Deviations from periodicity and sustained input signals
Experimental evidence suggests that brain oscillations in the gamma range are not perfect periodic oscillators with a consistent phase [36][37][38][39]. Consequently, to be a biologically plausible mode of communication, CTR should be robust enough to facilitate the transmission of oscillatory spiking activity when the constraint of a constant phase has been relaxed.
To quantify the extent to which CTR could afford unstable phases within an oscillation, we probed 10-layer diluted FFNs with periodic trains of pulse packets whose arrival times were jittered. The jitter was drawn from a uniform distribution centered on the arrival time (t n ) of the pulse packet. The extent of the jittering window was chosen to be a function of the interval t n +T=f where f~8, 4 or 2. The results showed that CTR could still enable the transmission in the presence of moderate amounts of jitter ( Figure 5a). For this particular selection of network parameters, a jitter of T=8 ms did not alter the main characteristics of the amplification process and the activity propagated to the last layer (Figure 5a top right).
However, if the jitter was further increased the activity propagated to fewer layers and the propagation was more unreliable. Interestingly, for a jitter of T=2, which corresponds to completely aperiodic pulse packet train, we observed that activity propagation increased with increasing the stimulus frequency. However, in this case also the pulse packets propagated with a frequency of *20 Hz, close to that of the network resonance frequency ( Figure S4b). That is, each FFN layer acted like a bandpass filter, which suggested that a broad-band noise stimulus could also trigger the transmission since it can generate oscillations close to the resonance frequency.
Indeed, it is well known that the dynamics of EI networks can display oscillations at the population level when they are stimulated with strong unstructured external drive [26,40]. We hypothesized that in the FNN a constant rate Poisson input could bring the activity of the first layer into an oscillatory regime, thereby generating a train of weak pulse packets that provide rhythmic input to the subsequent layers.
We tested this hypothesis by replacing the oscillatory input to P 1 by an additional source of constant Poisson input to all E neurons in the first layer. When in the network shown in Figure 5b-c the E drive was increased from 1 to 1.8 kHz the activity became oscillatory with enough power to ignite the resonance in the second layer. Interestingly, the frequency of the oscillations in P 1 was comparable to the resonance frequency of the network. This is not surprising as both resonance and oscillations at higher input regimes are shaped by the same network time constants, e.g., synaptic delays and membrane time constants [20].
Thus, we show that both slightly phase-jittered oscillatory inputs at the resonance frequency and broad-band stimulation are compatible with CTR in diluted FFNs.

Effect of the network state
Thus far, we have assumed that ongoing activity in each individual layer of the FFN was AI with low firing rates. However, there is ample experimental evidence suggesting that cortical networks in vivo can display more synchronized ongoing activity regimes depending on the behavioral state of the animal [30,41]. We therefore explored how the propagation of pulse packets via CTR is influenced by the dynamical state of the spontaneous network activity.
The level of synchrony in recurrent EI networks can be modulated by adjusting the firing rate of the external excitatory input [9,26,42]. Here, we changed the dynamical state by increasing the E drive from 1 to 1.6 Hz. Lower rate E drive gave rise to very sparse and asynchronous firing patterns, which progressively became more synchronous as the E drive was increased (synchrony measured as population Fano factor; red line in Figure 6a). The spiking activity of individual neurons remained irregular (CV ISI &1) for the parameter space explored here (cf. blue line in Figure 6a). PS o E increased in the range between 10 and 100 Hz for larger values of E drive. This increase was more pronounced around the peaks, which progressively shifted towards faster frequencies as the external input became stronger (Figure 6b).
To study the effect of network synchrony on CTR, we stimulated P 1 in 10-layer FFNs with periodic trains of pulse packets for the different levels of E drive and computed the signalto-noise ratio (SNR) in P 10 (cf. Methods). Generally, more synchronized activity states enabled CTR within a broader range of input frequencies, however the largest SNR values in P 10 were found at the low input regimes (Figure 6d). Independent of the synchrony level, resonance frequency and subsequently CTR were always confined within a range of input frequencies that closely matched the frequency around the peaks of PS o E (compare Figure 6d and Figure 6b). Hence, the resonance frequencies also became faster at higher levels of E drive. This shift reflected the reduction of the time that neurons needed to recover from the effective refractory state (absolute refractory period and hyperpolarization time) due to the presence of larger amounts of excitation as E drive was increased. The main peak in PS a E10 when activity reached P 10 was invariably found at * 20 Hz. This value was slower than the mean peak measured in PS o E10 which was * 28 Hz (Figure 6c). The values of PS a E10 were larger than those of PS o E10 for all the frequencies analyzed here. Notably, this difference was more pronounced in the gamma range (30{100 Hz) as compared to lower frequencies (v20 Hz; cf. Figure 6c).
Interestingly, network synchrony improved the propagation for faster input oscillatory regimes (*50{60 Hz, Figure 6d). In summary, our results showed that the ongoing state had opposite effects on CTR depending on the input frequency range. For lower input frequencies, AI activity increased SNR, while for larger input frequencies SI could enable the propagation which was absent during AI.

Model validation and implications for population coding
A direct validation of the model will involve the induction of coherent oscillations between distant brain areas by stimulating excitatory neurons in the presynaptic area at the resonance frequency. The resonance profile of a neuronal population can be obtained by recording its activity during periodic stimulation of the neurons with different frequencies. Similar experiments, which made use of optogenetic tools, have already been performed to study the role of specific cell types in the generation of gamma oscillations [43]. According to our model, even weakly connected distant networks (verified, e.g., by anatomical or electrophysiological studies) with a similar resonance profile can engage in a coherent oscillation by stimulating the presynaptic population at the resonance frequency. In contrast, a stimulation protocol, which does not induce a strong oscillation in the stimulated area, will fail to form such a coherent activity with the distant population. Our model also predicts a progressive entrainment characterized by a gradual increase in the measured power over multiple stimulation cycles in the stimulated presynaptic network. A similar entrainment should be found in the postsynaptic network with a certain delay which should be a function of the connectivity strength (see discussion). Moreover, in CTR mode of propagation the oscillations emerge only after a delay and not directly at the onset of the stimulus. This feature of the model is consistent with the observation that c -band oscillations appear after 100 ms of the stimulus onset (e.g. [44]).
This would confirm that CTR is by definition a slow mode of communication and therefore it is not suited for the communication of signals which have to propagate across multiple areas within a short period of time. Note that, e.g., in the FFN shown in Figure 4a, synchronous activity reached the fifth layer only after approximately 10 stimulation cycles (*200 ms at 40 Hz). We further quantified this result by testing the number of cycles required in a given layer until a significant synchronization level was found in the subsequent layer. A significant degree of synchrony was reached when the instantaneous rate of P neurons, computed using 5 ms time bins, hit a threshold value equal to n o P plus five times its s.d.. The results, computed using 100 trials, are shown in Figure 7 as a function of stimulus frequency (represented as the inter-pulse interval). Our results showed that when stimulated within the main resonance frequency range (39-42 ms intervals) the average speed of propagation was approximately two cycles/layer with small variability. Small deviations from that resonance frequency range resulted in higher trial-to-trial variability of the propagation speed and increased mean while larger deviations resulted in propagation failure.
The results obtained with our example FFN are indicative of how much time it will take to encode a stimulus using CTR at each stage of a processing chain. Naturally, the amount of time will be proportional to the number stages that the activity has to traverse. However, synchrony-based coding using FFNs seems to be suited only for communicating binary signals, i.e., the asynchronous/synchronous activity of a given layer indicates the absence/presence of a particular stimulus (e.g., a specific orientation of a bar of light). By contrast, the encoding of graded signals would require a monotonic relationship between the input and the output of the FFN. We tested the capacity of a diluted FFN to communicate continuous signals using CTR. To this end, we applied periodic stimuli with different amplitude a and computed the amplitude response of the network. Our results showed that for these network parameters it was possible to find an input range within which the system's response changed monotonically. Moreover the response remained linear for a restricted range of inputs strength a (cf. gray lines in Figure S2e). Such a linear operating regime including even a modest degree of saturation, could allow for the communication of graded signals.
We note that our model supports communication of activity between areas that have similar resonance profiles. This automatically ensures selective communication and gives possibility of gating the propagation by small change in the resonance frequency of a network. The experiments proposed above could demonstrate, whether CTR is in principle compatible with the neuronal hardware and physiology, even though they will not necessarily rule out other proposed mechanisms like CTC [17].

Discussion
Here we propose a novel mechanism for propagation of synchronous spiking activity within weakly coupled FFNs based on the presence of resonance in EI networks. In our model, resonance is a network property that emerges due to the interactions between excitatory and inhibitory neurons in each FFN layer. Using numerical simulations of spiking neuronal networks, we show that a weak and sustained stimulus can be gradually amplified in every layer, thereby overcoming the limitations of synchrony transmission imposed by the diluted interlayer connectivity. We refer to this mode of synchronous activity propagation as ''communication through resonance'' (CTR).
Until recently, resonance was considered mostly at the level of single cells in both experimental [45][46][47] and theoretical studies [48,49]. Now, there is increasing experimental evidence showing that resonance also exists at the network level in inhibitory [43] as well as excitatory neuronal populations [50], and may play a crucial role in the generation of cortical rhythms. Theoretical studies have shown that resonance is a fundamental property of EI networks [20] and could be used to gate neuronal signals [51].
In our model, such EI network resonance is used to enable the propagation of synchronous spiking activity in diluted FFNs. In previous theoretical studies, propagation of neuronal activity was restricted to either densely and weakly connected FFNs, which promote the propagation of synchronous activity [7,9,42,52,53], or sparsely and strongly connected FFNs, which are capable of propagating asynchronous firing ( [54,55]; see [14] for a review). However, biological neuronal networks are typically neither densely connected nor have strong synapses [11] and therefore the mechanisms that govern the propagation of neuronal activity in dense/strong FFNs are not always applicable. Our results indicate that propagation is possible in diluted FFNs, when aided by network resonance, but is restricted to synchronous activity.
Oscillations in the gamma range (30{70 Hz), which are a key feature of task-related population activity in several brain areas [56,57], have emerged as a prominent mechanism that may facilitate propagation of synchronous spiking activity in weakly connected networks [15]. These oscillations can synchronize neuronal activity and provide appropriate temporal windows of excitability, which enable communication between different brain areas. Within these temporal windows, effective functional connections are generated where otherwise only weak structural links may exist [17,58]. This mode of propagation, however, requires communicating brain areas to oscillate with matched phase and frequency (i.e., their oscillations are coherent) such that synchronous activity from the sender can reach the receiver during its excitable phase and maximize its spiking response. It is commonly believed that coherent oscillations are generated by two independent mechanisms, one responsible for the local generation of oscillations [59] and another mechanism that can flexibly modulate the coherence between spatially distant oscillators [19]. However, the precise nature of the process responsible for achieving such long-range coherence still remains elusive.
Here, we argue that coherent oscillations arise due to the propagation of periodic synchronous spiking activity. In our model, weak rhythmic synchronization provided by the input initially fails to propagate further down the FFN due to the diluted connectivity. The crucial role of the oscillations is to amplify this weak synchronous stimulus by promoting resonance dynamics of the receiving network and enable its propagation across the FFN. This is in contrast to the idea that oscillations are generated independently at every layer and locally synchronize unstructured background input. Our results show that oscillations arise in the network as a consequence of the stimulus propagation, and at the same time the stimulus exploits these oscillations to propagate. Due to this propagation, oscillations in each layer are driven by the previous layer and are hence naturally coherent with a phase that is determined by the conduction delay between the layers [17]. From this perspective coherence becomes a side effect of the propagation dynamics. Thus, a separation of distinct mechanisms that create oscillations and provide coherence is not necessary, as both arise naturally as consequence of CTR. Indeed, recent experimental studies suggest that there is an unidirectional entrainment of coherent oscillations between areas [60][61][62], making the feedforward spread of coherent oscillatory activity, as explained by our model, biologically plausible.
We show that while CTR still works for moderate deviations from periodicity, it is most efficient for propagating periodic stimuli. Notably, the same FFN architecture can transform a sustained firing rate signal into a weak rhythmic stimulus that can then be propagated. Even though it can be argued that environmental stimuli are often not periodic, it has been recently suggested that sensory information could be actively converted into periodic signals by sensing organisms [63,64].
CTR requires amplification of activity in each layer and, as a consequence, the propagation is slow requiring several cycles to reach the target network. The numbers of cycles needed to transmit synchronous activity across the entire FFN is a function of the connectivity strength between the layers. As the synaptic weights become stronger, the number of cycles required to spread synchrony to the final layer of the FFN decreases and transmission becomes more reliable. Once the weights are sufficiently strong, synchrony flows through the network in one oscillation cycle, which is equivalent to the propagation of synchronous activity in dense/strong connected FFNs investigated by previous studies (cf. [14] for a review). Thus, CTR could generate FFNs with strong connections capable of propagating isolated synchronous events, when certain types of synaptic plasticity are recruited to strengthen the synapses between the different FFN layers. Indeed, coherent oscillations, like those generated by CTR, can provide an ideal dynamical environment to promote synaptic potentiation [65]. In this way, CTR could be regarded as an initial means to propagate activity before strong connections have been formed, while providing the ideal substrate for the generation of fast and reliable communication channels.
In the present study, we describe activity propagation in single FFNs. However, other more complicated network architectures in which multiple FFNs interact may also be possible. In such a scenario, the input could create a stronger response in one such FFN, while partially and weakly activating other FFNs with unmatched resonance frequencies, thereby generating a broadband increase in power around the resonance frequency of the activated FFN. Thus, such a scheme could indeed explain the increase in broadband gamma power of the LFP signal observed during behavioral tasks [66].
Signal gating is an intrinsic property of CTR, since in a given FFN only the stimuli that match its resonance frequency are able to propagate. Selective gating of signals through network resonance has been suggested by previous theoretical studies [51,67]. Interestingly, the resonance frequency of the network can be dynamically modulated offering the possibility to gate signals differently in time. In our study, we show that modifying the level of external excitation shifts the resonance frequency of the FFN. Additionally, other mechanisms such as neuromodulator mediated changes of the effective connectivity within each layer can have similar effects on the resonance properties of the network. Another alternative gating mechanism is the use of gating signals [14]. Gating activity in dense/strong FFNs requires highly precise and strong gating signals [53]. However, the fact that in CTR the initial phase of the propagation in a given layer is characterized by low amplitude synchrony, which is still insufficient to elicit responses in the next layer, makes CTR suited for a gating mechanism that utilizes relatively imprecise and weak gating signals. Thus, overall CTR constitutes a flexible process that could implement complex spatio-temporal routing of neuronal signals.
As we show here, the dynamical properties of the background activity affect the quality (SNR) of the neuronal signals that are communicated using CTR. More specifically, SNR at low frequency stimulation (*20{30 Hz) was maximized when background activity state was asynchronous-irregular. This result is in line with experimental evidence which found oscillations in the gamma range to be associated with cortical desynchronization [68,69]. In contrast, the propagation of *50 Hz stimuli was successful only when ongoing activity was in a synchronousirregular state. These findings hint at a hypothetical scenario in which slow periodic modulations of the background dynamics could rhythmically improve or even gate signals that propagate using fast oscillations. The fact that the nesting of slow and fast cortical oscillations (e.g., beta-gamma) is commonly found in experiments (see [70] for a review) could be indicative of such a collaborative effort between different cortical rhythms. These findings open up the possibility that top-down signals may provide the change of background activity state required for coherent feedforward oscillations to be generated. Importantly, CTR is not restricted to the specific neuron and network model used in this work. The resonance mechanism, which is the essence of the model, is a general property of recurrently connected populations of excitatory and inhibitory neurons [20] and therefore it is widely applicable. Notably, a specific range of propagating frequencies can be achieved by a proper selection of network parameters. In summary, we have shown that communication of neuronal signals across weakly connected networks can be achieved by combining oscillatory activity with resonance dynamics.  Figure S3 Alternative FFN architectures that support CTR. (a) All E i{1 neurons were equally likely to establish longrange connections with the next layer but only P i neurons were allowed to receive projections from the previous layer. E PP~0 :1 which implied effectively increasing the total number of connections with respect to the architecture used in the main text from 9,000 to 30,000 synapses. The plots bellow the schematic drawing of the FFN are analogous to the ones used in Figure 4b in the main text. (b) In contrast to (a), all E i neurons were allowed to receive afferents from the downstream layer and project to the next layer. The total number of connections was kept as in (a) which implied a low E PP~0 :03.