Figures
Abstract
Zero-lag synchronization between distant cortical areas has been observed in a diversity of experimental data sets and between many different regions of the brain. Several computational mechanisms have been proposed to account for such isochronous synchronization in the presence of long conduction delays: Of these, the phenomenon of “dynamical relaying” – a mechanism that relies on a specific network motif – has proven to be the most robust with respect to parameter mismatch and system noise. Surprisingly, despite a contrary belief in the community, the common driving motif is an unreliable means of establishing zero-lag synchrony. Although dynamical relaying has been validated in empirical and computational studies, the deeper dynamical mechanisms and comparison to dynamics on other motifs is lacking. By systematically comparing synchronization on a variety of small motifs, we establish that the presence of a single reciprocally connected pair – a “resonance pair” – plays a crucial role in disambiguating those motifs that foster zero-lag synchrony in the presence of conduction delays (such as dynamical relaying) from those that do not (such as the common driving triad). Remarkably, minor structural changes to the common driving motif that incorporate a reciprocal pair recover robust zero-lag synchrony. The findings are observed in computational models of spiking neurons, populations of spiking neurons and neural mass models, and arise whether the oscillatory systems are periodic, chaotic, noise-free or driven by stochastic inputs. The influence of the resonance pair is also robust to parameter mismatch and asymmetrical time delays amongst the elements of the motif. We call this manner of facilitating zero-lag synchrony resonance-induced synchronization, outline the conditions for its occurrence, and propose that it may be a general mechanism to promote zero-lag synchrony in the brain.
Author Summary
Understanding large-scale neuronal dynamics – and how they relate to the cortical anatomy – is one of the key areas of neuroscience research. Despite a wealth of recent research, the key principles of this relationship have yet to be established. Here we employ computational modeling to study neuronal dynamics on small subgraphs – or motifs – across a hierarchy of spatial scales. We establish a novel organizing principle that we term a “resonance pair” (two mutually coupled nodes), which promotes stable, zero-lag synchrony amongst motif nodes. The bidirectional coupling between a resonance pair acts to mutually adjust their dynamics onto a common and relatively stable synchronized regime, which then propagates and stabilizes the synchronization of other nodes within the motif. Remarkably, we find that this effect can propagate along chains of coupled nodes and hence holds the potential to promote stable zero-lag synchrony in larger sub-networks of cortical systems. Our findings hence suggest a potential unifying account of the existence of zero-lag synchrony, an important phenomenon that may underlie crucial cognitive processes in the brain. Moreover, such pairs of mutually coupled oscillators are found in a wide variety of physical and biological systems suggesting a new, broadly relevant and unifying principle.
Citation: Gollo LL, Mirasso C, Sporns O, Breakspear M (2014) Mechanisms of Zero-Lag Synchronization in Cortical Motifs. PLoS Comput Biol 10(4): e1003548. https://doi.org/10.1371/journal.pcbi.1003548
Editor: Boris S. Gutkin, École Normale Supérieure, College de France, CNRS, France
Received: April 16, 2013; Accepted: February 20, 2014; Published: April 24, 2014
Copyright: © 2014 Gollo et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: This work was supported by the Australian Research Council TS0668690, BrainNRG JSMF220020255 and NMHRC PG1037196. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
The study of large-scale brain dynamics, and the cortical networks on which they unfold, is a very active research area, providing new insights into the mechanisms of functional integration and complementing the traditional focus on functional specialization in the brain [1], [2]. Whilst progress towards understanding the underlying network structure has been impressive [3], [4], the emergent network dynamics and the constraints exerted on these dynamics by the network structure remain poorly understood [5]. The problem is certainly not straightforward, as the dynamics between just a pair of neural regions already depends critically on the nature of the local dynamics and the nature of the coupling between them [6]: Although non-trivial, a complete description of nonlinear dynamics between a pair of nodes is nonetheless typically possible [7]. However, aggregating such duplets into larger arrays and introducing noise and time delays leads to further challenges and prohibits an exact description of the precise functional repertoire, motivating recourse to the broader objective of finding unifying and simplifying principles [8]. Structural and functional motifs – small subnetworks of larger complex systems – represent such a principle [9]. As depicted in Fig. 1 a, they characterise an intermediate scale of organization between individual nodes and large-scale networks that may play a crucial role as elementary building blocks of many biological systems [10]. Motif distribution in cortical networks has also been shown to be highly non-random, with a small set of motifs that appear to be significantly enriched in brain networks [9]. The relative occurrence of 3-node motifs in three different anatomical networks of the Macaque brain and cat cortex (Figs. 1 b–e) is shown in Figs. 1 f–i. These motifs may play distinct roles in supporting various computational processes. In this report we examine the principles of neuronal dynamics that emerge on small motifs and consider their putative role in neuronal function.
(a) The thirteen different motifs of size 3. (b–e) Connectivity matrices, and (f–i) Structural motif counts for each cortical network. Data (from the CoCoMac database [67], [68]) and algorithms are available at the brain connectivity toolbox website [69].
The mechanisms supporting zero-lag synchrony between spatially remote cortical regions can be considered paradigmatic of those mediating between structure and function. Since first reported in cat visual cortex [11], zero-lag synchrony has been widely documented in empirical data and ascribed a range of crucial neuronal functions, from perceptual integration to the execution of coordinated motor behaviours [12]–[16]. In particular, zero-lag synchrony between populations of neurons (quantified through synchrony between the local field potentials) may play a crucial role in aligning packets of spikes into critical windows to maximize the reliability of information transmission at the neuronal level [17], and to bring mis-aligned spikes into the time window of spike-time-dependent plasticity [18]. The situation is particularly pertinent in sensory systems, where precise differences in the timing of inputs, between left and right cortex for example, may carry crucial information about the spatial location of the perceptual source [19]. However, the empirical occurrence of zero-lag synchronization is at apparent odds with the observation that two mutually coupled oscillators interacting through a time-delayed connection do not, in general, exhibit zero-lag synchrony [20]. Indeed, in many models of neuronal systems the presence of a reciprocal delay has been found to introduce a ‘frustration’ into the system such that zero-lag synchrony is unstable and out-of-phase synchrony is instead the preferred dynamic relationship [21]. In fact, this phenomenon occurs quite generally in systems of oscillators with time-delayed coupling [21], [22].
Complex dynamics in spatially embedded systems arise in a broad variety of physical and biological contexts. Arrays of coupled semiconductor lasers are a prominent example. Because of their extraordinary internal speed, even small time delays due to the finite speed of light are usually nonnegligible in arrays of coupled lasers [23]. Detailed analysis of delay-coupled laser systems has suggested that an intermediate and reciprocally coupled relay node in a motif of three nodes could represent a general mechanism for promoting zero-lag synchrony in delay-coupled systems [24]. In previous work, it was also shown that such motif arrangements also represent a candidate mechanism for zero-lag synchrony in delay-coupled neuronal systems [25]. This is encouraging because there exist several candidate neuronal circuits in the mammalian brain which are characterized by reciprocal coupling between an intermediate delay node, including corticothalamic loops and the hippocampus [26], [27]. There also exist strong reciprocal connections in the visual system, such as the heavily myelinated connections between primary visual cortex and the frontal eye fields. Indeed, the corresponding motif occurs disproportionally in mammalian cortex (Fig. 1), hence being embedded in many cortical subsystems [9].
The presence of a node that drives two common-driven nodes that reach zero-lag synchrony between them due to the driver's influence is intuitively appealing and finds anatomical support, for example, by shared input through bifurcating axons [13]. Certainly, a common-driving input of sufficient intensity can generate virtually perfect spike-time correlation, as long as the time delay to both driven nodes is identical. However, this scenario is not robust if the time delays lose symmetry or the coupling is not sufficiently strong. The common-driving setup is nonetheless a key prototype that offers insights into the synchronization between the driven nodes and the roles of the dynamics of the nodes [28]–[32]. Here we consider dynamics on the 3-node motifs that occur abundantly in large-scale networks of the brain (Fig. 1), adding connections to the prototypical common-driving motif. We confirm that common driving – a coupling arrangement that is widely invoked in the literature – is an ineffective means of inducing zero-lag synchrony in the presence of weak coupling (a neurophysiologically plausible regime). However, the additional incorporation of a single reciprocally coupled connection between the driver and an edge node – which leads to synchrony between that pair – is found to be a novel and efficient way of promoting zero-lag synchrony amongst other nodes in these small motifs. We demonstrate that this effect – which we term resonance-induced synchrony – arises consistently in candidate computational models at the neuronal, population and mesoscopic spatial scales and is robust to mismatches in system parameters and even time delays. Remarkably, we show that the resonance effects of a synchronized pair are not necessarily localized, but may instead propagate throughout the network. We hence propose resonance-induced synchrony as a general and unifying mechanism of facilitating zero-lag synchrony in the brain.
Results
We studied zero-lag synchronization – quantified as the average zero-lag cross-correlation between two nodes A and B () – in a variety of different motifs involving a common driving node. We considered the dynamics of nodes expressing different neuronal systems across a hierarchy of scales. At the microscopic scale, each node was modeled to represent a single spiking Hodgkin-Huxley neuron; at the circuit scale, each node was taken to represent a population of 400 excitatory and 100 inhibitory randomly connected neurons described by the Izhikevich model; and at the mesoscopic scale each node was modeled as a neural mass model with chaotic activity. This last model permits systematic parameter exploration that is not possible with populations of spiking neurons. In all the three modeling levels, coupling between nodes was via excitatory chemical synapses (see Methods for details on models and integration scheme). For the sake of simplicity, we initially assumed homogeneous delays in the motifs, i.e., all connections between nodes had the same time delay. We later explored the robustness of the results when relaxing these assumptions in the section “Mismatch in the conduction delays”.
The notation we adopt for the motifs of three nodes follows the notation of Sporns and Kötter (2004) [9] who denoted all 13 possible connected subgraphs (motifs) composed of three nodes, denoted from M1 to M13 (Fig. 1). The genuine common-driving motif (illustrated in Fig. 2) is designated M3. Node 2 is the common driver whereas nodes 1 and 3 are the common-driven ones. In particular, we pay special attention to the cross-correlation between nodes 1 and 3. With the exception of illustrative time traces and their corresponding analysis, the results represent an average over 40 independent runs, unless otherwise stated, with different random initial conditions and with error bars given by the corresponding standard deviation. We characterize the synchronization in other motifs that represent structural variations of the M3 motif: the addition of one or more connections (M6, M8, M9, M13), or the addition of connections and nodes (e.g., M3+1). In particular, M9 is the prototypical dynamical-relaying motif [24], which has been previously shown to promote zero-lag synchronization in a variety of systems [33]–[38], including neuronal systems [25]–[27].
Dynamics of common driving motif (M3) versus common driving motifs with resonant sources (M6, M9 and M3+1) in motifs of excitatory delayed-coupled neurons with delay . First and second columns (panels a, b, e, f, i, j, m, n, q, r) correspond to individual spiking time traces of neurons, whereas the third and forth columns (panels c, d, g, h, k, l, o, p, s, t) correspond to the cross-correlation functions of the corresponding single time series and average over 40 trials respectively. Descending rows show motifs M3, M6, M9 and M3+1, respectively.
Common-driving motifs without and with resonance pairs
We first focus on the four motifs depicted in Fig. 2. The simple common driving motif (M3), in which node 2 drives the dynamics of nodes 1 and 3 was contrasted with three other motifs (M6, M9 and M3+1), which represent structural variations of M3. Because motif M3 lacks any feedback or cyclical structure, the conduction delay plays no role in the dynamics or in the synchronization between nodes 1 and 3: Hence the outer nodes passively receive the driver's input. Onto this “backbone”, motif M6 has a single feedback connection added, forming a reciprocal connection between nodes 1 and 2. Motif M9 has reciprocal connections between node 2 and nodes 1 and 3. Motif M3+1 possesses an extra node (4) reciprocally connected with node 2.
Motifs of Hodgkin-Huxley neurons.
For the smallest-scale system we consider, each node comprises a single excitatory Hodgkin-Huxley neuron [39] weakly inter-connected with a conduction delay of 6 ms. Each neuron receives independent Poisson trains of spikes, representing background stochastic input. Stimulated by such external input, neurons exhibit continuous spiking behavior with average inter-spike interval of approximately 15 ms and are hence suprathreshold, regardless of the input from the other motif neurons. When the neurons are coupled according to the M3 motif, as shown at the top row of Fig. 2, spikes from the center neuron 2 only sporadically trigger simultaneous spikes of neurons 1 and 3 (following the common 6 ms delay). Panels a and b illustrate an exemplar time trace of the neurons. As a consequence of the absence of regular coincident spikes in the outer neurons, the maximum of the cross-correlation between nodes 1 and 3 is small, evident in both the single trial (Fig. 2 c), and average (Fig. 2 d) results. Note in particular that the cross-correlation function between the central neuron and an outer neuron has a modest peak corresponding to the 6 ms time delay (blue trace in third and fourth columns). The smaller peak at zero lag (red trace) reflects this common time delayed peak from the center to each of the two outer nodes. On the other hand, when the neurons are coupled according to the structural variations of the M3 motif (namely M6, M9 and M3+1), as shown in the second to the fifth rows of Fig. 2, spikes from neuron 2 reliably trigger simultaneous spikes in neurons 1 and 3. This is evident in the exemplar time series as well as the single and average cross-correlation functions. This is quite a striking change, given that all other parameters of the model remain unchanged from M3.
The structural variations introduced in motifs M6, M9 and M3+1 over the common driving M3 share an essential feature: The driver node 2 is mutually connected and synchronized with at least one other node. We denote this mutual connection resonance pair: Its presence dramatically alters the dynamics and synchronization properties of the driven nodes. Supplementary Fig. S1 compares the dynamics of two oscillators with different types of time-delayed coupling. Synchronization between these pairs appears exclusively when they are mutually coupled (in this particular case the synchronization is in anti-phase at the slow rhythm). Therefore, the resonance pair, identified by the red stars in the motifs, is the source of resonance-induced synchronization, leading to zero-lag synchronization between the outer nodes.
The emergence of zero-lag synchrony in motifs of coupled Hodgkin-Huxley neurons shows a strong dependence on the time delay, consistent with prior work [40]. This delay effect is crucial to the dynamics of motifs containing a reciprocal coupling, but not for the common driving M3. To compare the dynamics of motifs M3 and M6 it is instructive to analyze both the cross-correlations between pairs of nodes and the regularity of the inter-spike intervals (ISIs). To measure the irregularity of the inter-spike intervals, we use the incoherence R, defined as the coefficient of variation (CV) of the ISI, ([41], [42]), where STD stands for the standard deviation. Large values of R indicate more irregular patterns of ISIs. As shown in Fig. 3 a, the incoherence of each node in M3 is independent of the delay and is larger for the driver node. In contrast, the incoherence of each node is similar in M6 and shows a strong effect on the time delay (top panels of Figs. 3 b–d), increasing and decreasing with a period of approximately half of the average ISI.
(a) Incoherence in motif M3 does not depend on time delay. Colors indicate the different nodes. (b–d) Top panels show incoherence for M6, where colors represent different nodes, and bottom panels show crosscorrelations for M6 (blue) and M3 (black). Continuous lines indicate the cross-correlation coefficients at zero time lag, and dashed lines indicate the maximum cross-correlation coefficients for all time lags. Panels b, c and d represent pairs of nodes: 1–2, 1–3, and 2–3 respectively. Phase, anti-phase synchrony, and asynchrony can be found in motif M6 depending on the time delay (see exemplar time traces in supplementary Fig. S2). Results are averaged over 40 trials.
The bottom panels of Figs. 3 b–d compare the cross-correlation between pair of nodes at zero-lag (continuous lines) against the maximum across all time lags (dashed lines) for motifs M3 (black) and M6 (blue). For motif M3, the maximum cross-correlation does not depend on the time delay. The input from node 2 solely arrives at nodes 1 and 3 after different latency times, but this delay does not impact on the dynamics of the driven nodes. In contrast, for motif M6 the maximum cross-correlations (blue dashed lines) vary, with peaks that coincide with the minima of incoherence (top panels). The synchronization in the Hodgkin-Huxley model, which has only one oscillatory frequency, appears either in phase or in antiphase. For neighboring nodes (1–2 or 2–3), phase synchronization occurs when a peak of the maximum cross-correlation coincides with a peak of the cross-correlation at zero lag. Anti-phase synchronization occurs when a peak of the maximum cross-correlation coincides with a minimum of the cross-correlation at zero lag. Supplementary Fig. S2 illustrates example time traces of Hodgkin-Huxley neurons in M6 for anti-phase synchronization (), no synchronization (
), and phase synchronization (
).
The delay in a resonance pair can either enhance or reduce the synchronization. In motif M6 the driven nodes (1 and 3) synchronize whenever the resonance pair (1 and 2) synchronizes, whether this is in-phase or anti-phase synchrony. This corresponds to a drop in the incoherence R of the driving node. Hence, the synchronization between nodes 1 and 3 depends on the time delay in the resonance pairs because the synchronization between the resonance pair (1 and 2) – and thus the incoherence – also depends on the time delay. It therefore appears that synchronization between the reciprocally connected nodes leads to a more regular (less incoherent) output from the master node which then facilitates synchronization between this node and the other slave node.
Motifs of neuronal populations.
To investigate whether these results translate to neuronal activity at the next spatial scale, we exploited the computational parsimony of the neural model of Izhikevich [43], [44] to study populations of spiking neurons, each node comprised a population of (400) excitatory and (100) inhibitory randomly interconnected neurons [27], with each neuron in these populations receiving an independent Poisson spike train. Neurons of the same populations were synaptically coupled without conduction delay and with a latency of for (exclusively) excitatory inter-population connections. We focused on the dynamics of the ensemble mean membrane potential
of all neurons within each population. As shown in the time traces of Fig. 4, the activity of each population consists of two time scales, a higher frequency (
) brief network spikes and a lower frequency fluctuation (
) on which these transients typically recur. Notably, the dominant (low frequency) time scale – which does not occur in the single neuron system – is much longer than the conduction delays. Despite discrepancies in the time scales and nature of the dynamics, the zero-lag synchronization reported in Fig. 4 largely resembles that shown in Fig. 2. Dynamics on the common driving motif M3 between the central and outer nodes show a moderate time delayed cross-correlation (
) at approximately
and a corresponding weak to moderate zero-lag synchrony between the outer nodes (
). However zero-lag synchrony is substantially stronger on the motifs possessing at least one resonance pair (M6, M9, M3+1). Notably, the anti-phase relation between node 2 with respect to nodes 1 and 3 appears solely at the faster time scale, comparable to the delay period.
Panels (a–t) as per Fig. 2, for populations of 500 (400 excitatory and 100 inhibitory) spiking neurons and delay .
Similar to motifs of Hodgkin-Huxley neurons, synchronization of populations of spiking neurons also shows a dependence on the delay time between nodes in the presence of a resonance pair. Supplementary Fig. S3 shows that the incoherence (here using the inter-burst interval instead of the inter-spike interval) and cross-correlations in motif M3 do not depend on the delay (panel a). However, they vary considerably for motif M6. Supplementary Figs. S3 b–d show that large incoherence values for motif M6 correspond to the transition between the regimes of phase and anti-phase synchronization for neural populations. Supplementary Fig. S4 illustrates two cases of synchronization with one dominant oscillatory frequency (one in phase with , and another in anti-phase with
), and one case (
) of synchronization with two oscillatory frequencies that is in phase for the slow rhythm and in anti-phase for the fast rhythm. For any time delay, the synchronization between 1 and 3 is enhanced for motif M6 when compared to M3, and the synchronization between nodes 1 and 3 largely resembles the maximum synchronization between first neighbors (nodes 1 and 2, or 2 and 3).
Motifs of neural mass models.
To further study the robustness of the relationship between motifs and synchronization with respect to the underlying dynamical systems, we next utilized a neural mass model, which represents a reduced model of cortical dynamics. A neural mass model is a parsimonious representation of the dynamics of a very high-dimensional system, and replaces thousands of equations for each population of neurons with a small number (here only three) of nonlinear equations per node. These represent the dynamical behavior of the essential summary system statistics (here mean firing rate) and hence a reduced representation of spontaneous cortical dynamics. Here we employ a population representation of conductance-based model neurons [5], [6], [45], [46], as has been previously used to elucidate important features of large-scale brain dynamics [47], [48]. This system also breaks from the previous two scales studied above in that irregularity is dynamically generated (through endogenous chaotic dynamics within each mode) rather than introduced through external stochastic spikes.
The dynamics of these delayed-coupled neural masses shows chaotic oscillations fast dynamics () superimposed on slower return times of about 110 ms [6]. As shown in Fig. 5, the dynamics in this system clearly replicate those observed above, namely that zero-lag synchrony between nodes 1 and 3 was strongly and exclusively expressed in the motifs with resonant sources (M6, M9, M3+1). It can be seen that within the resonance pairs, node 2 is in anti-phase synchrony with nodes 1 and 3. Notably, however, the anti-phase relation typically occurs at a much slower time scale (110 ms) than the coupling delay (10 ms).
Panels (a–t) as per Fig. 2, but for neural mass models with coupling strength c = 0.01, and delay .
Despite the dissimilarities between the neuronal systems at different scales, synchronization and incoherence of the neural mass model also exhibits a dependence on the time delay in the presence of a resonance pair (see supplementary Fig. S5). To better understand synchronization dynamics in this system – which has multiple internal time scales – it is necessary to study the combination of time delays and coupling strength: For weak coupling strength (c = 0.01), phase synchronization is not reached. However, as illustrated in supplementary Fig. S6, rich dynamics can arise, including anti-phase synchronization at the slow () or fast (
) time scales, or an asynchronous state (
). In contract, when the coupling is stronger (Supplementary Figs. S7, S8, left), phase synchronization emerges for very short time delays. In the cases of stronger coupling, for example c = 0.05 or c = 0.15 (Figs. S7, S8, right), zero-lag synchronization between nodes 1 and 3 is also more stable for long delays.
Mismatch in the conduction delays
Biological systems are naturally diverse, and therefore, any relevant behavior should not be highly dependent on the fine-tuning of the delay – and particularly its symmetry. We next tested the generality of the zero-lag synchronization between nodes 1 and 3 with respect to delay mismatch in the different motifs containing the resonance pair. The connections preserved the conduction delay of except for a single feedback connection to the driver node 2 in motifs M6′, M9′ and M3+1′ in which we introduced a variable conduction delay in one direction (
), as illustrated in Fig. 6 a. The three motifs exhibited zero-lag synchronization that was substantially larger than that of motifs M3 (black line) or even M3 plus a unidirectional input (yellow line) across a large region of the parameter space (Figs. 6 b–d). In the motifs of Hodgkin-Huxley neurons (Fig. 6 b), the behaviors of all three motifs are similar for
. In contrast, for
zero-lag synchrony decays in a similar way for motifs M6′ and M3+1′, whereas synchronization in motif M9′ is virtually independent of
for up to fivefold
(not shown in the plot). Supplementary Fig. S9 shows the analyses of the dynamics of motif M6′ in more detail: It shows that synchronization arises in M6′ only when the delay mismatch
yields synchronization with the same phase relation as
, which – in the case of
– is anti-phase synchronization between neighboring neurons (see Fig. 3). The motifs of neural mass models show a systematic consistency of synchronization across
for a biologically plausible range of delays (Fig. 6 c). However, a behavior similar to that observed in motifs of Hodgkin-Huxley neurons occurs for greater delay mismatches (Fig. 6 d). Such differences in the time scales are consistent with the different time scales of these systems: The Hodgkin-Huxley neurons oscillate with periods of about 15 ms, whereas the neural masses oscillate with periods of about 110 ms.
The top schemes (a) illustrate the motifs of neurons considered. Motifs M6′, M9′ and M3+1′ have one connection with delay , and all the other connections have delays of 6 ms. The bottom panels show the zero-lag crosscorrelation between nodes 1 and 3 in motifs of Hodgkin-Huxley neurons (b) and in motifs of neural mass models with c = 0.01 (c) averaged over 40 trials for varying
. Panel (d) shows the same as (c) but across a broader range of
. Plot colors correspond to motifs as per panel a.
Characterizing the dynamics of the motifs
From herein, we focus on motifs of neural masses, exploiting their relative computational parsimony to gain deeper insight into the mechanisms of the resonance pair. In particular we studied the robustness of our findings with respect to the most salient parameters of the system, namely the coupling strength and the delay. As shown in Fig. 7, the strength of the synchronization in the motifs with a resonance pair, but not M3, show an increase as a function of coupling strength (panels a, b). Although an expected feature of the model [49], the emergence of synchrony even at very weak coupling () is somewhat surprising for a biological system. There are, however, some regions of complex dynamics (evidenced as large error bars) in which there is not a unique solution, thereby entailing significant trial-to-trial variability. At relatively weak coupling (c = 0.01), zero-lag synchronization between nodes 1 and 3 holds across a broad regime of physiologically plausible time delays (Fig. 7 c). Analysis of longer coupling delays (supplementary Fig. S10) reveals an influence on synchronization that resembles the system of Hodgkin-Huxley neurons (Fig. 3), albeit weaker and at a much longer time scale.
Top: Common-driving motifs, labeled as per Sporns and Kötter (2004) [9], see Fig. c1 a. Bidirectional connections (red stars) indicate active resonance pairs. Top row panels compare common driving (M3) to common driving with resonance pairs (M6, M9 and M3+1) for varying coupling (panels a and b) and varying delay (panel c). Bottom row panels compare common driving (M3) to a ring of mutually coupled nodes (M13), and to common driving plus a bidirectional connection between 1 and 3 (M8) as a function of coupling (panels d and e) and time delay (panel f). Curve colors correspond to the motifs depicted on the top of the figure.
These analyses suggest a partition of the common-driving motifs into three distinct families: (i) The simple common-driving motif (M3) where synchronization at zero lag is not achieved in the weak-coupling regime, independent of the time delay; (ii) A ring of three mutually coupled systems (M13) or a common-driving motif that also contain direct coupling between the driven nodes (M8) require a relative strong coupling and negligible delay in order to promote synchronization (Figs. 7 d–f), because of the existence of frustration; and (iii) Common-driving motifs enhanced by active resonance pairs (e.g., M6, M9, M3+1) which exhibit zero-lag synchronization even for very small couplings, irrespective of the time delay (up to ). It is clear in these analyses that the increase in zero-lag synchrony in motifs with a resonance pair is not due to the additional coupling introduced by the backward connection, but rather through the placement of the additional edge. For example, the motifs with the greatest number of edges (M8 and M13) are amongst the most difficult to achieve zero-lag synchrony with an increase in coupling. Closing the outer nodes with two additional edges (going from M9 to M13) leads to a substantial decrease in zero-lag synchrony.
Propagation of the effect of the resonance pair
The preceding analyses show that the effect of the resonance pair can influence the common driving motif even when it is placed outside the motif itself (e.g. M3+1). Here we further investigate the propagation of the resonance pair effect by considering larger structures in which the resonance pair is distant from the driver node (2). This procedure is schematically shown in Fig. 8 a, and illustrated for a particular network of N = 7 nodes in Fig. 8 b. We are particularly interested to understand if the effects of the resonance pair are strictly local, and, additionally, on how the polysynaptic distance to the resonance pair influences the dynamics and synchronization.
(a) A resonance pair (nodes N and N-1) arbitrarily distant from a pair of commonly driven neural masses (1 and 3). (b) A seven-node chain configuration with a common-driving motif at the edge. (c) Zero-lag cross-correlation functions between nodes 1 and 3 for different chain sizes as illustrated in panel (a) are shown in solid lines. Thin dashed line represents the chain of panel (a) without the feedback connection from node N-1 to node N. (d) Zero-lag cross-correlation functions between every other node in the chain depicted in panel (b) are shown in solid lines. Thin dashed line represents the chain of panel (b) without the feedback connection from node 6 to node 7.
We observe that zero-lag synchronization between the driven nodes 1 and 3 is virtually independent of the distance along a polysynaptic chain from the resonance pair (Fig. 8 c). For a fixed motif length (N = 7), we also characterized the zero-lag synchronization of different pairs of nodes that did not interact directly, but interacted indirectly through a common neighboring mediator (see Fig. 8 d). Apart from pairs 5–7, all such pairs correspond to a strict flux of information flow, mandated by the direction of the coupling. Thereby, the synchronization decreased with the distance from node 7, unless the system was set with a specific coupling (see arrow in Fig. 8 d) that gives rise to global synchronization. This corresponds to identical synchronization between nodes 2, 5 and 7, which are anti-phase synchronized to nodes 1, 3, 4 and 6 occurring at this particular coupling strength.
Finally, to highlight the influence of the resonance pair in the dynamics, we removed the feedback connection to node N (results shown as thin dotted lines in Figs. 8 c and d). By means of this control simulation, we find that: (i) Zero-lag synchronization between 1–3 is consistently reduced (Fig. 8 c); and (ii) Zero-lag synchronization between 5–7 (Fig. 8 d) completely disappears in the absence of a resonance pair.
Characterizing active resonance pairs
We have denoted an active resonance pair as two mutually connected nodes that synchronize in the presence of appropriate time delays and coupling strength. This effect propagates through the motifs because the driven nodes show a strong tendency to synchronize with the driver node (hence promoting zero-lag synchronization between driven nodes). That is, the emergence of synchronization between the resonance pair then stabilizes synchrony amongst unidirectionally coupled nodes. The same phenomenon underlies the propagation down a polysynaptic chain (Fig. 8). Interestingly, the impact of the resonance pair extends beyond this propagation, giving rise to other dynamical effects for coupling delays in which anti-phase synchrony between neighbors prevails. Geometrical frustration is an example: In some motif configurations, anti-phase synchrony between pairs of mutually connected nodes (potential resonance pairs) is simply not a stable solution. In the case of motif M13 (illustrated in Fig. 7), for example, anti-phase synchronization between any pair is frustrated because the third node cannot be simultaneously synchronized in anti-phase with respect to the other two neighbor nodes. This situation illustrates that frustration can disturb potential resonance pairs. Large mismatches in the delays of the mutual connection between the pair can also disturb the effects of a resonance pair. As depicted in Fig. 6, both motifs M6′ and M3+1′ are similarly susceptible to mismatches in the reciprocal latencies.
Transient behavior and the stability of synchronization in resonance-pair motifs
Connectivity also plays a role on the onset of synchronization. We studied the temporal onset of zero-lag synchronization in neural mass models for different motifs by (1) examining the transient dynamics following random initial conditions, and (2) studying the response to a transient perturbation. An example is shown in Fig. 9 a, in which dynamics on M6 begin from random initial conditions, then approach synchronization between masses 1 and 3. The dynamics are then perturbed by a brief current from 800 to 1000 ms – that is distinct for each driven node – before rapidly regaining synchrony after a few hundreds of milliseconds. It is noteworthy that the approach to zero-lag synchrony in both scenarios is approximately exponential, with an exponent that can be used as a numerical estimate of the stability of the synchronous state (Fig. 9 b). In contrast, edge nodes on motif M3 remain unsynchronized. The dependence of the exponent
with the coupling strength for the 1200 ms following offset of the transient perturbation is shown in Fig. 9 c. Motifs with resonance pairs (M6 and M9) showed a negative exponent, consistent with stable synchrony, whereas the exponent associated with motif M3 was positive throughout. Interesting, the coupling strength associated with the strongest synchrony (most negative exponent) occurred for a relatively weak coupling strength of c = 0.01.
(a) Example of time-trace synchronization following random initial conditions (starting at time = 0) and consequent to a brief perturbing current (green bar) at time = 1000 ms in motif M6 with c = 0.01. (b) averaged over 400 trials with c = 0.01 in motifs M3 compared to M6. (c) Exponent
estimated from
averaged over 400 trials on the interval between 1200 and 2400 ms for varying coupling strengths in motifs M3, M6 and M9. Delay
.
Fine-tuning and synchrony in the absence of resonance pairs
Synchronization hence arises quickly in the presence of a resonance pair. Is it possible to adjust the dynamics of the driver node without such reciprocal coupling to induce synchronization? We next studied this possibility by fine-tuning the input current ( to the driver node (2) in motif M3, whilst keeping all other parameters fixed. As shown in Fig. 10 a, introducing a slight mismatch in the input current can indeed lead to large changes in the zero-lag synchronization between nodes 1 and 3. Crucially, careful fine-tuning of this current mismatch can lead to a near complete synchronization in motif M3 (A), or at least lead to a strong enhancement of synchronization (B and C). As depicted in Fig. 10 b, the maximum synchronization (A) occurs when the input current causes the driver node to exhibit the same oscillatory frequency as the driven edge nodes. The other local maxima occur when the driver node oscillates with a frequency that is an integer multiple of the driven nodes (2∶1 in B and 3∶1 in C). In contrast to this need for fine-tuning in motif M3, the resonance pair guarantees that node 2 oscillates with the same frequency as the driven nodes, with strong synchronization hence arising regardless of the coupling strength, as shown in Fig. 10 c for motif M3+1.
(a) Crosscorrelation averaged over 40 trials, (b) dominant oscillatory frequencies of neural masses 1 (green) and 2 (magenta) as a function of the mismatch on the input current over node 2. (c) Dominant oscillatory frequencies of neural masses 1 (green) and 2 (magenta) for varying coupling strength.
Beyond resonance pairs
The effects of a resonance pair can enhance the synchronization locally and even propagate in a polysynaptic way to influence distant dynamics. Reciprocally connected nodes can also interact in a way that disturbs the synchronization if they introduce frustration as in motifs M8 and M13, as shown in Fig. 7. To more deeply understand the role of reciprocally connected nodes and loops, we studied resonance motifs that go beyond the resonance pairs. Starting with a common driving motif M3, we added chains of bi- or uni-directionally coupled nodes of varying sizes as shown in Fig. 11 a. Adding one node reciprocally connected to node 2 recovers the resonance pair, which is clearly a more effective way of synchronizing the driven nodes than adding one extra unidirectionally connected node (the blue dashed line of Fig. 11 b). The addition of two reciprocally connected extra nodes in a closed loop (resonance triplet) had an effect that was analogous to the resonance pair, and again far more effective than the counterpart of two extra unidirectionally connected nodes in a loop (green dashed line). The addition of three or more reciprocally connected extra nodes in closed chain had a similar effect to the resonance pair. However, the influence of the unidirectionally coupled loops gradually approaches that of their reciprocally connected counterparts, which have already attained the ceiling effect (magenta dashed line). Hence, the interaction of unidirectionally connected nodes in a loop gradually enhances the synchronization of the driven nodes as the size of the loop increases. Therefore, even in the absence of reciprocally connected nodes, synchronization between 1 and 3 can be enhanced by a loop of at least three extra nodes connected to the driver node. Interestingly, the addition of a single resonance pair is the most efficient means of achieving zero-lag synchronization compared to loops of any size.
(a) Loops of reciprocally connected versus unidirectional connected loops. (b) Zero-lag cross-correlation between neural masses 1 and 3 with neural mass 2 connected to bidirectional or unidirectional chains of varying length. Blue dashed line highlights the effect of the resonance pair, and green (magenta) dashed line highlights the effect of the resonance triplet (quad). Red (yellow) curve represents the cross-correlation averaged over 40 trials for reciprocally (unidirectionally) connected loops. The coupling strength is 0.01, and delay 10 ms.
Effects of the common driving input at higher orders
Our final analysis concerns the synchronization properties of commonly driven nodes with higher polysynaptic orders (Figs. 12 a–c). In particular, we study the synchronization of the symmetrically located nodes n−n′ for the different connectivity states of the driver node A. Figure 12 a illustrates the case in which node A was part of a resonance pair together with node B; Fig. 12 b illustrates the case in which node A received a unidirectional input from node B; Fig. 12 c illustrates the case in which node A did not receive input from any neighboring regions. It can be seen in Figs. 12 d–g that only the motifs with the resonance pair (red line) yielded high correlation between nodes n and n′ (for n = 1,2,3,4). Interestingly, when the coupling strength is fixed (c = 0.024) and the number of elements further increased (Fig. 12 h), the cross-correlation coefficient remained quite high for the chain containing the resonance pair. A similar behavior occurred for the maximum cross-correlation coefficient (for all time delays) between node A and node n (Fig. 12 i): Again, the resonance pair was required for the propagation of synchronous activity.
Common driving to first (1,1′), second (2,2′) and n-th (n,n′) order for the resonance-induced pair (a), a unidirectional input (b), and simple common driving (c). (d) to (g): Zero-lag crosscorrelation for the different types of common driving from the first to the forth order versus the coupling strength. (h) Zero-lag crosscorrelations between pairs of nodes (n,n′) as a function of the distance from the driver node A. (i) Maximum (non-zero-lag) crosscorrelations as a function of the distance from the driver node A. Red, yellow and black curves represent the crosscorrelation averaged over 40 trials for the system depicted in (a), (b) and (c) respectively.
Discussion
Zero-lag synchronization between distant neuronal populations confers a number of important computational advantages, and finds broad empirical support. Here we report that common driving of passive nodes by a central “master” (motif M3), a scenario that is broadly assumed to underlie zero-lag synchrony, fails completely in the weak-coupling regime and is sensitive to parameter mismatch. However, the addition of one or more mutually coupled pairs fosters the emergence of zero-lag synchrony in the outer nodes of triplet motifs, and beyond. We find that this effect is robust to many of the particular details of the system, the spatial scale and parameter asymmetry, and can propagate through a multi-synaptic relay chain. In stark contrast, the further addition of a reciprocal connection between the driven nodes introduces frustration for delays that favor out-of-phase synchrony and fails to promote zero-lag synchronization. The disruptive effect of adding new edges that close the motif reinforces the observation that it is the topology (not the total amount of coupling) that determines the zero-lag synchrony. This is also evident by the fact that an increase in the coupling over two orders of magnitude in the unidirectional motif (M3) is less effective than adding a single feedback connection (where the effective coupling within that pair is simply doubled).
We have denoted this reciprocal pair a resonance pair because it can induce zero-lag synchronization between outer nodes. We find that an entire family of three- and four-node motifs exhibits zero-lag synchronization in the presence of such a resonance pair. Perhaps the archetypal motif in this family is M9 (see Fig. 1) also known as the dynamical relaying motif [24]–[27], [33]–[38], [49]. This motif contains two active resonance pairs (Fig. 1). Here we find that one feedback connection to the driver node can be removed (i.e., transforming the motif into M6) without compromising the synchronization between the outer nodes (confirming a recent observation in electronic circuits [50]). Similarly, the addition of one extra node mutually connected to the driver node, M3+1 (thereby comprising a resonance pair) causes robust zero-lag synchronization of the driven nodes where M3 alone fails. This indicates that a necessary condition for nodes 1 and 3 to synchronize is that the resonance-pair nodes also synchronize, regardless of their exact phase relationship. The synchronization of the resonance pair appears in turn to enhance its propensity to synchronize the driven nodes because when the driving node is synchronized its internal incoherence diminishes: This change in the regularity of the master node in turn enslaves the unilaterally driven node onto the synchronization manifold (Fig. 9). Thereby, we propose that the mechanism that promotes zero-lag synchronization in the dynamical relaying motif is indeed the resonance pair, in common to all other motifs in the broader family we examined.
We observed the effect of the resonance pair in a variety of different models (Hodgkin-Huxley neurons, populations of Izhikevich neurons, and neural mass models) and scales: motifs of neurons and motifs of cortical regions. The results are also robust with respect to the delay, the coupling strength, the oscillatory frequency band, and arise in autonomous, chaotic systems as well as noise-driven excitable dynamics. It seems reasonable to propose that resonance-induced synchronization will prove important for other neuronal systems, such as dendritic oscillations in single-neuron dynamics [51], and indeed other physical and biological systems of any domain characterized by weak interactions. Although the responses of neural populations to noisy inputs have been well studied [52], it remains to be seen if our results prove robust to further physiological details, including embedding stronger synaptic inputs into the noisy background [53] and stronger balanced background inhibitory and excitatory inputs [30]. We also note that although our study focused mainly on interactions with time delay, the resonance-induced synchronization can also occur in systems with no time delay (Figs. 7 b and c, and supplementary Figs. S3, S5, S7, S8 and S10).
Despite the robustness of the present effect in different classes of models and dynamical regimes, the universality and extent of the phenomenon remains to be clarified. Phase-resetting curves (PRCs) can be useful to predict whether phase or out-of-phase synchronization will arise [54]: This is a crucial factor in the dynamics because frustration does not occur in the case of in-phase synchronization. While usually studied in systems without delay, PRCs can also be used in systems in the presence of conduction delays [25], [55]. Analysis of the PRC can also be employed for formal stability analysis of synchronization of motif dynamics [25]. A second caveat, at least in the model of population of spiking neurons, is the type of dynamics studied – namely that in the dynamical regime studied here, neurons spike at least once per population cycle. An alternative approach would be to analyze synchronization in motifs of populations of spiking neurons in a sparsely synchronized regime [56] – that is when individual neurons spike less often than the background ensemble cycle. Further analysis is hence required to elucidate the extent to which our results translate to other physical and biological systems, perhaps focusing on canonical models that are more amenable to mathematical analysis such as the Kuramoto system.
Computational studies of anatomically derived brain networks have shown that motifs M9 and M6 are the first and second most abundant of all three-node motifs in the macaque visual cortex [9] and are among the most frequent motifs in other cortical networks (Fig. 1). Moreover, they appear to be clustered around the core “rich club” backbone of the structural connectome [57]. The presence of a resonant-pair in these motifs, and the robust zero-lag synchrony that they confer, may provide a dynamical advantage for these pairs. However, given the additional wiring cost, it is not clear why motif M9 is more common than M6. A possible explanation we provide derives from our observation that synchronization on motif M9 is robust to longer delays in one branch of the resonance pair in comparison to M6 (Fig. 6). Hence, the gain in robustness might overcome the cost of maintaining this extra feedback connection.
The influence of a resonance pair is not limited to local synchronization dynamics but also, through propagation, to larger networks, decaying only slowly with the polysynaptic distance (see Figs. 8 and 12). In a sufficiently sparse network like the brain, the number of neurons grows roughly exponentially with the inter-node distance. The coexistence of the slow decay (long correlation length) of the influence of the resonance pair, with rapid growth in the number of affected elements as a function of synaptic distance suggests that the zero-lag synchronization arising locally through a resonance pair has the capability to impact globally on network dynamics. Reframed in terms of a branching process, the slow decay of zero-lag synchronization and rapid growth of neuronal connectivity could lead to critical or supercritical propagation of zero-lag synchrony, consistent with prior theoretical considerations [58], and also suggesting a means for analytic extension of the present results.
The notion of motifs as fundamental building blocks of complex networks has yielded considerable prior success [9], [10], [53]. Degree distribution, the relative density of reciprocal synapses, convergence, divergence, and chains of synapses have been shown to play a crucial role in shaping the dynamics and synchronization properties of large networks [59]–[62]. In contrast to these studies, which focus on the global statistical features of large-scale networks, we have focused on particular features of small motifs. Future work, aimed at immersing these small motifs into larger networks, and focusing on the role of reciprocal nodes on the global synchronization properties of such networks, would be of significant interest. Our work confirms that the interplay between structural, functional and effective connectivity, while likely complex [63], may nonetheless be reliant upon a small number of unifying principles.
Methods
We simulated neuronal motif dynamics at different scales, and for different dynamical scenarios. First, representing the microscopic scale, each node was taken as a single neuron. For this endeavor we utilized the Hodgkin-Huxley model. Second, at the circuit scale, we took each node as a large population of spiking neurons. Third, at the mesoscopic scale, we considered a simplified coarse-grained version in which each population was taken as a neural mass model.
Hodgkin-Huxley neurons
Each node was modeled by the well-known Hodgkin-Huxley equations [39]. The dynamics of the membrane potential depends on sodium, potassium, leaky, and synaptic (intra-motif and external) current components,(1)where
is the membrane capacitance. The maximal conductances of the channels occur for completely open channels, with conductances given by
,
, and
, and
,
, and
stand for the corresponding reversal potentials. Generally, the voltage-gated ionic channels are not fully opened. The probability of finding them open depends on the gating variables. The
channel depends on the combined effect of gating variables
and
, whereas
depends on
. They evolve according to the equations,
(2)
(3)
(4)Hodgkin and Huxley set the empirical functions
and
to fit the experimental data of the squid giant axon,
(5)
(6)
(7)
(8)
(9)
(10)The synaptic current due to the interactions between neurons of the motifs are given by,
(11)where
,
, and
stands for the Dirac delta function. The summation over
stands for the spikes of the presynaptic neurons (all excitatory).
is the time at which the
spike occurred. We varied the conduction delay
. In agreement with the literature [40], the delay
can shape the synchronization (Fig. 3). The external current incoming to each neuron is,
(12)where
,
,
runs over excitatory spikes, and
corresponds to the spike times, modeled by an independent Poisson process for each neuron with rate
. As shown in supplementary Fig. S11, nearly identical results can be also obtained by assuming the external current term as synaptic contribution and including it as an extra term in equation (4) with
, and
. The equations were integrated by the Runge-Kutta method of fourth order, with time steps of
. Initial transient dynamics were discarded.
Populations of Izhikevich neurons
For this large-scale circuit model, each node represented populations of 500 randomly connected neurons described by the Izhikevich model [43]. 400 neurons were excitatory and 100 neurons were inhibitory. The neurons were described by the following equations:(13)where
represents the membrane potential,
represents the recovery variable, accounting for the
and
ionic currents, and
is the total synaptic current. The neurons have a threshold at
. Once this value is reached,
is reset to
and
to
. Following [44], we added dispersion to these four parameters (
,
,
and
) to account for neuronal heterogeneity. Excitatory neurons have
, and
, where
is a random number drawn from a uniform distribution in the interval [0,1]. Inhibitory neurons have
, and
.
Each neuron receives input from 80 neurons of the same population and from 25 excitatory neurons of each afferent population. The synaptic current is given by(14)where the dynamics of the excitatory and inhibitory synapses are described by
(15)
in the equations above stands for the Dirac delta function. The summation over
(
) stands for the spikes of the presynaptic excitatory (inhibitory) neurons.
(
) is the time at which the
excitatory (or
inhibitory) spike occurred. Conduction delays
, associated with excitatory long-range connections, varied. We modeled short-range (intra-node) connections with negligible delays. Synapses were modeled by exponential decay functions [64], with time constants
for excitatory and
for inhibitory synapses. Each neuron was subject to an external driving given by independent Poisson spike trains at a rate of
, which was also included in the sum over excitatory postsynaptic contributions (
index) of the equations above. With these parameters, individual neurons fire spontaneously, although not periodically.
The equations were integrated using a fixed-step first-order Euler method with time steps of 0.05 ms, starting with random initial conditions. To avoid spurious synchronization at the onset of simulations, neural populations were activated with random noise in 600 ms sequential windows (with a 500 ms overlap). The first transients of 1 s were discarded before further analysis.
Neural mass models
The preceding large-scale circuit model is a high dimensional system. Whilst the dynamics are instructive, the large number of parameters and equations preclude an intuitive perspective of the system. We therefore additionally studied a reduced system [65], which represents the large cortical scale that permits characterization of the system dynamics with respect to the most salient parameters. In contrast to the previous models, the coupling is not through discrete pulses, but by means of smooth sigmoidal rate functions, which embody population-wide neuronal responses to synaptic inputs in the presence of parameter and state dispersion [66]. This also allows us to study the robustness of the resonance-induced synchronization in relationship to the precise details – and dynamical regime – of the models.
Each node represents the mean dynamics of an ensemble of neurons, with spontaneous dynamics arising from the interaction between excitatory and the inhibitory sub-populations. The model is derived from the biophysical Morris-Lecar model [45], extended to a neural mass model with passive diffusive chemical [46], then synaptic interactions [6] and subsequently extended to large networks to model whole brain activity [5]. We utilize this most recent approach developed by Honey et al. [5], [47] systematically varying the features of the connectivity: architecture, coupling strength, and delay.
This neural mass model comprises three state variables: The mean membrane potential of the excitatory pyramidal neurons, ; the mean membrane potential of the inhibitory interneurons,
; and the average number of open potassium ion channels,
. Our main focus is on the dynamics of the pyramidal neurons. Their average membrane potential V depends on the passive leak conductance, and on the conductance of voltage-gated channels of sodium, potassium and calcium ions. The flow of current across the local pyramidal cell membranes, assumed as capacitors, governs its dynamics. In turn, the local activity of the inhibitory interneurons is course-grained modeled; its dynamics is modulated by the activity of the pyramidal cell. For each ensemble
, the equations for the dynamics of the mean membrane potential of the neurons are given by
(16)
(17)
The fraction of channels open are the neural-activation function, whose shape reflects a sigmoidal-saturating grow with
(18)The third differential equation of each node
stands for the fraction of open potassium channels:
(19)The neuronal firing rates (
, and
) averaged over the ensemble are assumed to obey Gaussian distributions, thereby giving rise to the sigmoidal activation functions [66],
(20)
(21)
Our simulations employ the previously published parameter values: ,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
were set to physiological values taken from [6]. These are associated with aperiodic fluctuations arising without external noise, but rather due to homoclinic chaos [6]. Equation 16 includes the other important parameters in our analysis: the presynaptic neighboring (afferent) regions of region
; c, the coupling strength between cortical regions;
, the synaptic delay between cortical regions. The model was simulated in Matlab (Math Works) using the function dde23.
Supporting Information
Figure S1.
Dynamics of pairs of neural mass models. (a), (d) and (g) show the time traces of the average membrane potential of the excitatory pyramidal neurons; (b), (e) and (h) show the auto-correlation function of node A; (c), (f) and (i) show the cross-correlation function between nodes A and B; respectively for a pair bidirectionally connected, unidirectionally connected, and disconnected nodes. Parameters are c = 0.01, and .
https://doi.org/10.1371/journal.pcbi.1003548.s001
(TIFF)
Figure S2.
Example dynamics of Hodgkin-Huxley neurons coupled on motif M6 for different time delays. From left to right, panels show anti-phase synchronization (), no synchronization (
), and phase synchronization (
).
https://doi.org/10.1371/journal.pcbi.1003548.s002
(TIFF)
Figure S3.
Synchronization dynamics and incoherence in populations of Izhikevich neurons. Panels (a–d) as per Fig. 3 but for populations of spiking neurons. Phase, anti-phase synchrony, and a state of phase synchrony at the slow rhythm and anti-phase synchrony at the fast rhythm can be found in motif M6 depending on the time delay (see exemplar time traces in supplementary Fig. S4).
https://doi.org/10.1371/journal.pcbi.1003548.s003
(TIFF)
Figure S4.
Example dynamics of populations of Izhikevich neurons coupled as motif M6 for different time delays. From left to right, panels show phase synchronization (), phase synchronization at the slow rhythm and anti-phase synchronization at the fast rhythm (
), and anti-phase synchronization (
).
https://doi.org/10.1371/journal.pcbi.1003548.s004
(TIFF)
Figure S5.
Synchronization dynamics and incoherence in weakly coupled neural mass models. Panels (a–d) as per Fig. 3 but for neural mass models with coupling strength c = 0.01. Anti-phase synchrony at the slow or at the fast rhythms, and a state of low synchrony can be found in motif M6 depending on the time delay (see exemplar time traces in supplementary Fig. S6).
https://doi.org/10.1371/journal.pcbi.1003548.s005
(TIFF)
Figure S6.
Example dynamics of neural mass models coupled on motif M6. From left to right, panels show anti-phase synchrony at the slow rhythm (no time delay), weak synchrony (), and anti-phase synchrony at the fast timescale (
). The coupling strength is weak, c = 0.01.
https://doi.org/10.1371/journal.pcbi.1003548.s006
(TIFF)
Figure S7.
Cross-correlation for strongly coupled neural mass models. Top panels show the crosscorrelations between nodes 1 and 2 (a–b), and nodes 1 and 3 (c–d) as a function of the delay for coupling strength c = 0.05. Bottom panels panels show the cross-correlations between nodes 1 and 2 (e–f), and nodes 1 and 3 (g–h) as a function of the delay for coupling strength c = 0.15. First and third columns correspond to a zoom of second and forth columns respectively. Black (blue) lines represent results for motif M3 (M6), and continuous (dashed) lines represent the crosscorrelation at zero lag (maximum for all time lags). Phase synchrony, and complex synchronous states can be found in motif M6 depending on the time delay (see exemplar time traces in supplementary Fig. S8). Results are averaged over 40 trials.
https://doi.org/10.1371/journal.pcbi.1003548.s007
(TIFF)
Figure S8.
Example dynamics of neural mass models strongly coupled as motif M6 for different time delays. From left to right, panels show phase synchrony (), and a transition from a state of anti-phase synchrony at the slow rhythm to a state of out-of-phase synchrony (
). The coupling strength is c = 0.15.
https://doi.org/10.1371/journal.pcbi.1003548.s008
(TIFF)
Figure S9.
Incoherence and synchronization dependence on the time-delay mismatch in Hodgkin-Huxley neurons. (a–c) Top panels show incoherence: Colors represent different nodes. Bottom panels show cross-correlations for motif M6′. Continuous lines indicate the cross-correlation coefficients at zero time lag, and dashed lines indicate the maximum cross-correlation coefficients across all time lags. Panels a, b and c represent pairs of nodes: 1–2, 1–3, and 2–3 respectively. Results are averaged over 40 trials.
https://doi.org/10.1371/journal.pcbi.1003548.s009
(TIFF)
Figure S10.
Zero-lag synchronization dependence on the delay in motifs of neural mass models. The curves are color coded as in Fig. 7. In agreement with [40] the synchronization depends on the coupling delay for long delays. The coupling strength is c = 0.01. Crosscorrelation is averaged over 40 trials.
https://doi.org/10.1371/journal.pcbi.1003548.s010
(TIFF)
Figure S11.
Kernel test in motifs of Hodgkin-Huxley neurons. Nearly identical cross-correlation functions are obtained when the external driving is considered identical to the spikes within the motifs (see Fig. 2, panels d and h). Plot corresponds to an average over 40 trials.
https://doi.org/10.1371/journal.pcbi.1003548.s011
(TIFF)
Acknowledgments
We thank Christopher J. Honey for technical assistance and helpful discussions, and Andrew Zalesky for discussions and a careful reading of the manuscript.
Author Contributions
Conceived and designed the experiments: LLG CM OS MB. Performed the experiments: LLG. Analyzed the data: LLG CM OS MB. Contributed reagents/materials/analysis tools: LLG CM OS MB. Wrote the paper: LLG CM OS MB.
References
- 1. Bullmore E, Sporns O (2009) Complex brain networks: graph theoretical analysis of structural and functional systems. Nature Reviews Neuroscience 10: 186–198.
- 2. Friston KJ (2011) Functional and effective connectivity: a review. Brain Connectivity 1: 13–36.
- 3. Hagmann P, Cammoun L, Gigandet X, Meuli R, Honey CJ, et al. (2008) Mapping the structural core of human cerebral cortex. PLoS Biol 6: e159.
- 4.
Sporns O (2010) Networks of the Brain, volume 1. MIT Press, 375 pp.
- 5. Honey C, Kötter R, Breakspear M, Sporns O (2007) Network structure of cerebral cortex shapes functional connectivity on multiple time scales. Proc Natl Acad Sci 104: 10240–10245.
- 6. Breakspear M, Terry JR, Friston KJ (2003) Modulation of excitatory synaptic coupling facilitates synchronization and complex dynamics in a biophysical model of neuronal dynamics. Network: Computation in Neural Systems 14: 703–732.
- 7. Ashwin P, Buescu J, Stewart I (1999) From attractor to chaotic saddle: a tale of transverse instability. Nonlinearity 9: 703.
- 8. Breakspear M, McIntosh AR (2011) Networks, noise and models: reconceptualizing the brain as a complex, distributed system. Neuroimage 58: 293.
- 9. Sporns O, Kötter R (2004) Motifs in brain networks. PLoS Biology 2: e369.
- 10. Milo R, Shen-Orr S, Itzkovitz S, Kashtan N, Chklovskii D, et al. (2002) Simple building blocks of complex networks. Science 298: 824–827.
- 11. Gray CM, König P, Engel AK, Singer W (1989) Oscillatory responses in cat visual cortex exhibit inter-columnar synchronization which reflects global stimulus properties. Nature 338: 334–337.
- 12. Roelfsema P, Engel A, Konig P, Singer W (1997) Visuomotor integration is associated with zero time-lag synchronization among cortical areas. Nature 385: 157–161.
- 13. Singer W (1999) Neuronal synchrony: A versatile code for the definition of relations. Neuron 24: 49–65.
- 14. Varela F, Lachaux J, Rodriguez E, Martinerie J (2001) The brainweb: phase synchronization and large-scale integration. Nat Rev Neurosci 2: 229–239.
- 15. Fries P (2005) A mechanism for cognitive dynamics: neuronal communication through neuronal coherence. Trends in Cognitive Sciences 9: 474–480.
- 16. Uhlhaas P, Pipa G, Lima B, Melloni L, Neuenschwander S, et al. (2009) Neural synchrony in cortical networks: history, concept and current status. Front Integr Neurosci 3: 17.
- 17. Fries P, Nikolić D, Singer W (2007) The gamma cycle. Trends Neurosci 30: 309–316.
- 18. Masquelier T, Hugues E, Deco G, Thorpe SJ (2009) Oscillations, phase-of-firing coding, and spike timing-dependent plasticity: an efficient learning scheme. Journal of Neuroscience 29: 13484–13493.
- 19. Nolte G, Bai O, Wheaton L, Mari Z, Vorbach S, et al. (2004) Identifying true brain interaction from EEG data using the imaginary part of coherency. Clinical Neurophysiology 115: 2292–2307.
- 20. Heil T, Fischer I, Elsässer W, Mulet J, Mirasso CR (2001) Chaos synchronization and spontaneous symmetry-breaking in symmetrically delay-coupled semiconductor lasers. Phys Rev Lett 86: 795–798.
- 21. Breakspear M, Heitmann S, Daffertshofer A (2010) Generative models of cortical oscillations: neurobiological implications of the kuramoto model. Frontiers in human neuroscience 4: 1–14.
- 22. Acebrón JA, Bonilla LL, Vicente CJP, Ritort F, Spigler R (2005) The kuramoto model: A simple paradigm for synchronization phenomena. Reviews of modern physics 77: 137.
- 23. Soriano MC, García-Ojalvo J, Mirasso CR, Fischer I (2013) Complex photonics: Dynamics and applications of delay-coupled semiconductors lasers. Reviews of Modern Physics 85: 421.
- 24. Fischer I, Vicente R, Buldú JM, Peil M, Mirasso CR, et al. (2006) Zero-lag long-range synchronization via dynamical relaying. Physical Review Letters 97: 123902.
- 25. Vicente R, Gollo LL, Mirasso CR, Fischer I, Pipa G (2008) Dynamical relaying can yield zero time lag neuronal synchrony despite long conduction delays. Proc Natl Acad Sci 105: 17157–17162.
- 26. Gollo LL, Mirasso C, Villa AEP (2010) Dynamic control for synchronization of separated cortical areas through thalamic relay. NeuroImage 52: 947–955.
- 27. Gollo LL, Mirasso CR, Atienza M, Crespo-Garcia M, Cantero JL (2011) Theta band zero-lag long-range cortical synchronization via hippocampal dynamical relaying. PLoS ONE 6: e17756.
- 28. de la Rocha J, Doiron B, Shea-Brown E, Josić K, Reyes A (2007) Correlation between neural spike trains increases with firing rate. Nature 448: 802–806.
- 29. Shea-Brown E, Josić K, de la Rocha J, Doiron B (2008) Correlation and synchrony transfer in integrate-and-fire neurons: Basic properties and consequences for coding. Phys Rev Lett 100: 108102.
- 30. Renart A, De La Rocha J, Bartho P, Hollender L, Parga N, et al. (2010) The asynchronous state in cortical circuits. Science 327: 587–90.
- 31. Barreiro AK, Shea-Brown E, Thilo EL (2010) Time scales of spike-train correlation for neural oscillators with common drive. Phys Rev E 81: 011916.
- 32. Tchumatchenko T, Malyshev A, Geisel T, Volgushev M, Wolf F (2010) Correlations and synchrony in threshold neuron models. Phys Rev Lett 104: 058102.
- 33. Gomes I, Buldú JM, Mirasso C, García-Ojalvo J (2006) Synchronization by dynamical relaying in electronic circuit arrays. Chaos 16: 043113.
- 34. Wagemakers A, Buldú JM, Sanjuan MAF (2007) Isochronous synchronization in mutually coupled chaotic circuits. Chaos: An Interdisciplinary Journal of Nonlinear Science 17: 023128.
- 35. Chitra R, Kuriakose V (2008) Phase effects on synchronization by dynamical relaying in delaycoupled systems. Chaos: An Interdisciplinary Journal of Nonlinear Science 18: 023129.
- 36. de Sousa Vieira M (2010) Properties of zero-lag long-range synchronization via dynamical relaying. Chaos: An Interdisciplinary Journal of Nonlinear Science 20: 013131.
- 37. Banerjee R, Ghosh D, Padmanaban E, Ramaswamy R, Pecora LM, et al. (2012) Enhancing synchrony in chaotic oscillators by dynamic relaying. Phys Rev E 85: 027201.
- 38. Jian-Ping L, Lian-Chun Y, Mei-Chen Y, Yong C (2012) Zero-lag synchronization in spatiotemporal chaotic systems with long range delay couplings. Chinese Physics Letters 29: 050501.
- 39. Hodgkin A, Huxley A (1952) A quantitative description of the membrane current and its application to conduction and excitation in nerve. Journal of Physiology 117: 500–544.
- 40. Dhamala M, Jirsa VK, Ding M (2004) Enhancement of neural synchrony by time delay. Phys Rev Lett 92: 074104.
- 41. Pikovsky AS, Kurths J (1997) Coherence resonance in a noise-driven excitable system. Physical Review Letters 78: 775.
- 42. Medeiros BN, Copelli M (2013) Synaptic symmetry increases coherence in a pair of excitable electronic neurons. PloS one 8: e82051.
- 43. Izhikevich EM (2003) Simple model of spiking neurons. IEEE Transactions on Neural Networks 14: 1569–72.
- 44. Izhikevich EM, Gally JA, Edelman GM (2004) Spike-timing dynamics of neuronal groups. Cerebral Cortex 14: 933–944.
- 45. Morris C, Lecar H (1981) Voltage oscillations in the barnacle giant muscle fiber. Biophys J 35: 193–213.
- 46. Larter R, Speelman B, Worth RM (1999) A coupled ordinary differential equation lattice model for the simulation of epileptic seizures. Chaos: An Interdisciplinary Journal of Nonlinear Science 9: 795–804.
- 47. Honey CJ, Sporns O, Cammoun L, Gigandet X, Thiran JP, et al. (2009) Predicting human restingstate functional connectivity from structural connectivity. Proc Natl Acad Sci 106: 2035–2040.
- 48. Alstott J, BreakspearM, Hagmann P, Cammoun L, Sporns O (2009) Modeling the impact of lesions in the human brain. PLoS computational biology 5: e1000408.
- 49. Vicente R, Fischer I, Mirasso CR (2008) Synchronization properties of three delay-coupled semiconductor lasers. Phys Rev E 78: 066202.
- 50. Suresh R, Srinivasan K, Senthilkumar D, Mohamed IR, Murali K, et al. (2013) Zero-lag synchronization in coupled time-delayed piecewise linear electronic circuits. The European Physical Journal Special Topics 222: 729–744.
- 51. Remme MWH, Lengyel M, Gutkin BS (2009) The role of ongoing dendritic oscillations in singleneuron dynamics. PLoS Comput Biol 5: e1000493.
- 52. Fourcaud N, Brunel N (2002) Dynamics of the firing probability of noisy integrate-and-fire neurons. Neural computation 14: 2057–2110.
- 53. Song S, Sjöström PJ, Reigl M, Nelson S, Chklovskii DB (2005) Highly nonrandom features of synaptic connectivity in local cortical circuits. PLoS biology 3: e68.
- 54. Van Vreeswijk C, Abbott L, Ermentrout GB (1994) When inhibition not excitation synchronizes neural firing. Journal of computational neuroscience 1: 313–321.
- 55. Goel P, Ermentrout B (2002) Synchrony, stability, and firing patterns in pulse-coupled oscillators. Physica D: Nonlinear Phenomena 163: 191–216.
- 56. Brunel N, Hakim V (2008) Sparsely synchronized neuronal oscillations. Chaos: An Interdisciplinary Journal of Nonlinear Science 18: 015113–015113.
- 57. Harriger L, van den Heuvel MP, Sporns O (2012) Rich club organization of macaque cerebral cortex and its role in network communication. PloS one 7: e46497.
- 58. Kumar A, Rotter S, Aertsen A (2010) Spiking activity propagation in neuronal networks: reconciling different perspectives on neural coding. Nature Reviews Neuroscience 11: 615–627.
- 59. Roxin A (2011) The role of degree distribution in shaping the dynamics in networks of sparsely connected spiking neurons. Frontiers in Computational Neuroscience 5: 1–15.
- 60. Zhao L, Bryce Beverlin I, Netoff T, Nykamp DQ (2011) Synchronization from second order network connectivity statistics. Frontiers in Computational Neuroscience 5: 1–16.
- 61. Trousdale J, Hu Y, Shea-Brown E, Josić K (2012) Impact of network structure and cellular response on spike time correlations. PLoS computational biology 8: e1002408.
- 62. Hu Y, Trousdale J, Josić K, Shea-Brown E (2013) Motif statistics and spike correlations in neuronal networks. Journal of Statistical Mechanics: Theory and Experiment 2013: P03012.
- 63. Adachi Y, Osada T, Sporns O, Watanabe T, Matsui T, et al. (2012) Functional connectivity between anatomically unconnected areas is shaped by collective network-level effects in the macaque cortex. Cerebral Cortex 22: 1586–1592.
- 64.
Dayan P, Abbott LF (2001) Theoretical neuroscience, volume 31. Cambridge, Massachusetts: MIT press.
- 65. Breakspear M, Terry JR (2002) Nonlinear interdependence in neural systems: motivation, theory, and relevance. The International journal of neuroscience 112: 1263–1284.
- 66. Marreiros AC, Daunizeau J, Kiebel SJ, Friston KJ (2008) Population dynamics: variance and the sigmoid activation function. Neuroimage 42: 147–157.
- 67. Stephan KE, Kamper L, Bozkurt A, Burns GA, Young MP, et al. (2001) Advanced database methodology for the collation of connectivity data on the macaque brain (cocomac). Philosophical Transactions of the Royal Society of London Series B: Biological Sciences 356: 1159–1186.
- 68. Kötter R (2004) Online retrieval, processing, and visualization of primate connectivity data from the cocomac database. Neuroinformatics 2: 127–144.
- 69. Rubinov M, Sporns O (2010) Complex network measures of brain connectivity: uses and interpretations. Neuroimage 52: 1059–1069.