Synchronisation in networks of delay-coupled type-I excitable systems

Synchronisation in networks of delay-coupled type-I excitable systems

Andrew Keane 1Institut für Theoretische Physik, Technische Universität Berlin, 10623 Berlin, Germany 1    Thomas Dahms 1Institut für Theoretische Physik, Technische Universität Berlin, 10623 Berlin, Germany 1    Judith Lehnert 1Institut für Theoretische Physik, Technische Universität Berlin, 10623 Berlin, Germany 1    Sachin Aralasurali Suryanarayana 1Institut für Theoretische Physik, Technische Universität Berlin, 10623 Berlin, Germany 12Department of Physics, Indian Institute of Technology Bombay, Mumbai, India 2    Philipp Hövel 1Institut für Theoretische Physik, Technische Universität Berlin, 10623 Berlin, Germany 13Bernstein Center for Computational Neuroscience, Humboldt-Universität zu Berlin, Philippstraße 13, 10115 Berlin, Germany 34Center for Complex Network Research, Northeastern University, 110 Forsyth Street, Boston, MA 02115, USA 4    Eckehard Schöll 1Institut für Theoretische Physik, Technische Universität Berlin, 10623 Berlin, Germany 1
September 19, 2019

We use a generic model for type-I excitability (known as the SNIPER or SNIC model) to describe the local dynamics of nodes within a network in the presence of non-zero coupling delays. Utilising the method of the Master Stability Function, we investigate the stability of the zero-lag synchronised dynamics of the network nodes and its dependence on the two coupling parameters, namely the coupling strength and delay time. Unlike in the FitzHugh-Nagumo model (a model for type-II excitability), there are parameter ranges where the stability of synchronisation depends on the coupling strength and delay time. One important implication of these results is that there exist complex networks for which the adding of inhibitory links in a small-world fashion may not only lead to a loss of stable synchronisation, but may also restabilise synchronisation or introduce multiple transitions between synchronisation and desynchronisation. To underline the scope of our results, we show using the Stuart-Landau model that such multiple transitions do not only occur in excitable systems, but also in oscillatory ones.

1 Introduction

Studies on complex networks have gained much attention in recent times from various fields of research ALB02a ; NEW03 ; BOC06a . In particular, the dynamics on networks of nonlinear elements (or nodes) has become a very active area TIM02a ; STE06a ; ASH07 ; FLU10b ; LEH11 ; ENG11 ; KAN11 ; OME11 ; DAH12 ; HAG12 ; NEP12 . Here, we consider the dynamics of neurons as a network of nonlinear excitable nodes. While the model we consider is generic for excitability, the local dynamics being modelled in this context is the “firing” mechanism of neurons (i.e. the rise and fall of the electrical membrane potential). This mechanism is the basis for neural coding and information transfer or cell-to-cell communication PER68 . The type of network we primarily focus on in our investigations is the small-world (SW) network WAT98 , which has a short average path distance between nodes, as well as a large degree of clustering (in other words, many triangles in the network structure). These properties are found in many kinds of real-world structures, such as the collaboration between film actors, power grids, the World Wide Web and social relationships WAT98 ; ADA99 ; ALB02a ; NEW03 . In particular, large-scale cortical networks also show these properties SPO00 ; SPO06 . The brain has an architecture enabling both efficient global and local communication between neurons LAT01 , which is captured well by the SW model.

The interest in the synchronisation of neurons emerges from the role it plays in information processing and perception, as well as in epilepsy and similar diseases TRA82 ; ROE97 ; ENG01a ; PIK01 ; MEL07 . While there are several types of synchronisation, the focus here will be on zero-lag synchronisation. This means that the dynamics of all nodes in the network are not only identical but also temporally in-phase.

Inhibition plays an important role in the nervous system HAI06 . Here, when constructing a network we begin with a regular ring network of excitatory links and, as in Ref. LEH11 , we add long-range inhibitory links into the network structure. This creates a SW network of the form proposed in Ref. NEW99b . The introduction of inhibitory links and an increasing of their density in the network was shown in Ref. LEH11 to result in a loss of stable synchronisation.

In this paper, we consider delay-coupled networks of type-I excitability, in contrast to Ref. LEH11 where type-II was considered using FitzHugh-Nagumo local dynamics, and find qualitative differences for a particular range of coupling parameters. The most interesting implication of this is that for some network topologies with this type of excitability increasing the density of inhibitory links may not just lead to one transition from stable to unstable synchronisation, but may lead to multiple transitions. This means that increasing the density of inhibitory links in the network may stabilise synchronisation, rather than destabilise it.

In the following section, we introduce the model and the network topologies considered in the present study. In Sec. 3 we use the master stability function to investigate the stability of arbitrary synchronised networks with given coupling parameters (i.e. the coupling strength and the length of delay between coupled nodes). The master stability function is calculated in Sec. 4 for networks coupled within a range of small delay time and coupling strength, which delivers the main results of the paper - namely, the existence of synchronised states that have different stability conditions compared to coupling with larger delay times. The implications these results may have for specific complex networks are discussed in Sec. 5. In Sec. 6, we study networks of coupled Stuart-Landau oscillators and demonstrate that multiple transitions between synchronisation and desynchronisation can also appear in networks of oscillatory nodes and are not limited to excitable systems.

2 Model

In order to model excitability, the system must have a rest state, which corresponds to a stable fixed point. Small perturbations from the rest state allow for the creation of a large excursion in the phase space, which corresponds to the system entering an excited state before returning to the rest state. In the context of neurodynamics, this is the firing state of the neuron IZH00 .

Neurons can exhibit different excitability properties. In 1948, Hodgkin classified two types of neural excitability HOD48 . Both types can be ordered depending on the type of bifurcations that occurs when the bifurcation parameter is changed such that the system makes a transition from the excitable to the oscillatory regime:

Type-I neurons can generate action potentials of arbitrarily low frequency. This kind of behavior occurs near a saddle-node infinite period (SNIPER) bifurcation, also known as the SNIC bifurcation (saddle-node bifurcation on invariant cycle). The arbitrarily low frequency coincides with the period of the limit cycle going to infinity as the bifurcation parameter approaches a critical value, where the bifurcation occurs.

Type-II neurons are associated with a supercritical Hopf bifurcation. The frequencies of the action potentials lie within a certain non-zero range, while the amplitude of the limit cycle approaches zero with the bifurcation parameter.

As our model for type-I excitability we consider a generic normal-form of a SNIPER bifurcation HU93a ; HIZ07 . This model is mathematically represented by two first-order differential equations:


In polar coordinates Eq. (1) reads


is the bifurcation parameter. This bifurcation parameter influences the type of dynamics and determines where in the -plane the fixed points are located, as discussed below.

An important feature that must be taken into account in the mathematical model is that the effect one neuron has on another is not instantaneous. The signals sent between neurons have a finite propagation speed and take a certain time to arrive at the destination neuron, giving rise to a delay time .

The dynamics of a network of elements (labelled ) is written:


where is the local dynamics as described by Eq. (1) for each element . determines the matrix for the network structure, showing which elements are coupling together, and is the coupling function. is taken to be the identity matrix; this means that the -variable at time is coupled with the -variable at time , the -variable at time is coupled with the -variable at time , but there is no cross-coupling between the - and -variable. The coupling parameters, which are identical for all connections, are the coupling strength and delay time .

For synchronous dynamics, the constraints define a 2-dimensional synchronisation manifold (SM) within the -dimensional phase space and the coupled system is reduced to an effective single system with feedback:


with unity row sum , so that the nodes all receive the same level of input while they are synchronised. Any non-unity but constant row sum can be rescaled using the coupling strength .

With the bifurcation parameter there exists an unstable focus at the origin, as well as a saddle point and a stable node situated on the unit circle at and , respectively. At the saddle point and stable node collide, so that for a limit cycle exists along the unit circle. In the case of (excitable regime), a perturbation which pushes the system from the stable node beyond the saddle point can result in a single oscillations along the unit circle. Delayed coupling can induce a homoclinic bifurcation, such that a limit cycle is produced that bypasses the saddle point and stable node HIZ07 ; AUS09 . Here, we will consider the excitable regime with and investigate for which topologies the synchronised dynamics of complex networks is stable.

In this paper, two kinds of complex networks are considered: SW networks and Erdős-Rényi random networks ERD59 . We construct the SW networks as a variation to the method proposed in Ref. WAT98 , introduced in Refs. MON99 ; NEW99b : (i) Each of the nodes in a one-dimensional ring is given excitatory links to its nearest neighbours on each side. Note that in terms of the matrix , an excitatory link between the and node means that , while for an inhibitory link . (ii) For each of the links we add with a probability of another inhibitory link with strength -1 between two randomly chosen nodes (i.e. on average randomly distributed inhibitory links). (iii) Self-coupling and multiple links between the same two nodes are not allowed. (iv) Finally, the entries in each row of are normalised to ensure a unity row sum. If a row sum is equal to zero, then the network realisation is discarded.

The random networks are constructed with excitatory links between randomly chosen nodes, as well as randomly distributed inhibitory links, so that these may be compared to a SW network with the same parameters , and . Because the construction of both kinds of complex networks involves random processes, there will be many possible realisations of either kind of network with a given set of parameters.

In the following, we determine stability of the synchronised dynamics and compare the stability of these different kinds of networks.

3 Stability of Synchronisation

The Master Stability Function (MSF) is a tool that represents a measure of the stability of a synchronised state for given delay-coupling parameters in terms of the topology of an arbitrary network PEC98 . The master stability equation, given by


is found by linearising Eq. (4) around Eq. (5) and can be used to calculate the largest Lyapunov exponent , called the MSF. is the perturbation of away from the SM (i.e. ) and is the Jacobian matrix of Eq. (1) evaluated on the SM. Important for this approach is that, whereas one can calculate Lyapunov exponents for a specific network topology using the eigenvalues of the matrix , one considers here a continuous complex parameter representing the complex plane of possible eigenvalues scaled by the coupling strength (i.e. is a continuous parametrisation of , where are the eigenvalues of ). Thus, one can calculate the Lyapunov exponents for a region of the -plane which gives sub-regions of stability, where , and instability, where . It is then easy to compare the synchronous stability of various networks by simply observing whether any of their eigenvalues fall inside an unstable region of the -plane. If all the eigenvalues lie within stable regions, then perturbations away from the SM will decay exponentially.

Because of the unity row sum condition, always has an eigenvalue . This longitudinal eigenvalue (all others are called transversal) corresponds to perturbations within the SM and is always zero because we are looking at periodic dynamics. As such, it is only the transversal eigenvalues that are important for determining the stability of synchronisation.

Figure 1(a) shows the MSF of the SNIPER model with coupling parameters and . The white dots represent the eigenvalues of the matrix for a unidirectional ring of 11 nodes. All eigenvalues are within the stable region of the MSF, thus the synchronisation of all 11 nodes coupled with these parameters is stable.


[width=0.99height=0.82]Images/blank.jpeg (d)

Figure 1: (Colour online) Master stability function for coupling parameters and (a) , (b) , (c) , and (d) . and are the real and imaginary parts of the scaled eigenvalues of the coupling matrix, , respectively, and is the largest Lyapunov exponent. White dots in panel (a) represent the scaled eigenvalues of a unidirectional ring of 11 nodes (). .

According to Ref. FLU10b , for in the order of the system’s characteristic time scale (in this case, the period of the oscillations) and above, the MSF will always tend towards a rotational symmetry. Examples such as the one above in Fig. 1(a) confirm these general findings for the SNIPER model. When calculating MSFs for a large fixed while varying it becomes evident that the size of the stable region is scaled by , so that the stable region can be estimated very well by the circle (that is, a circle centred at the origin with a radius ). Large values of , however, despite influencing the Lyapunov exponents quantitatively, do not affect the shape of the stable region.

This rotational symmetry of the MSF was also found for type-II neurons, modelled as FitzHugh-Nagumo oscillators LEH11 . In this case and do not affect whether the eigenvalues fit into the stable region of the MSF, so that only the topology of a network is important for the stability of its synchronisation. Furthermore, because Gershgorin’s circle theorem GER31 guarantees that the eigenvalues of a network’s coupling matrix with no self-coupling and purely excitatory coupling (i.e. and ) lie within the unit circle on the complex plane, the synchronisation of such a network will always be stable. Finally, if additional inhibitory links are introduced with probability to construct a small-world network as described in Sec. 2, phase transitions from stable to unstable synchronisation are found with increasing probability of inhibition LEH11 . All these results also apply for the SNIPER neurons with sufficiently large .

As an example, Fig. 2 shows the transition from stable synchronisation, where the fraction of network realisations with unstable synchronisation is zero, to unstable synchronisation, where , as the probability of inhibition is increased. This example is for the coupling parameters and , but this transition to inhibition-induced desynchronization will look the same for all sufficiently large values of . A thorough investigation of the effects of additional inhibitory links when is small will be presented in Sec. 5.

Figure 2: (Colour online) Fraction of desynchronised networks vs. probability of inhibitory links for 100 realisations of SW networks of and for , , .

4 Small delay times

If the delay time is not large enough, the rotational symmetry of the MSF no longer holds. In fact, while reducing one can witness how the rotational symmetry begins to break down. This is depicted in Fig. 1. For a constant coupling strength of , the MSFs are numerically calculated for decreasing delay times. While at the MSF still has its circular form (Fig. 1(a)), when decreasing , the stable (i.e. dark blue/green) region of the MSF begins to show signs of deformation. By (see Fig. 1(b)) the stable region is clearly larger than the unit circle scaled by and has definitely lost its rotational symmetry. By (see Fig. 1(c)) the stable region has already split into two disconnected regions. Letting decrease further, the stable regions become increasingly smaller (see Fig. 1(d)). Note that the delay-induced limit cycle coexists alongside the stable fixed point and whether it is reached or not is therefore dependent on initial conditions. For less than about 4 (not shown here), the coupling is no longer able to induce the homoclinic bifurcation that creates the limit cycle at all (as discussed in Sec. 2), in other words, the neurons no longer oscillate. Instead, all solutions approach the stable fixed point.

It is now obvious that, in this regime of small delay, small changes in can have a great impact on the stability. Despite the seemingly sudden change in the MSFs in Fig. 1 between and , the evolution of the boundary of stability is in fact a continuous one (except at a critical value , which will be discussed below). This becomes clear after plotting the MSF versus the real part of the eigenvalue (with ) for varying . Taking this one slice of the eigenvalue plane gives a good indication of the growth and decay of the stable region in the MSF while changing . Figure 3 shows the MSF as a function of the real part () and the delay time , with fixed coupling strengths of (a) 0.25, (b) 0.30 and (c) 0.35. This produces an MSF on an -plane, which can be used for network topologies with a coupling matrix that yields real eigenvalues, e.g. symmetric matrices for undirected networks. In the following we restrict ourselves to undirected networks.


[width=0.9]Images/tau46e-1-12_sigma25e-2_alpha_tau_c.pdf (a){overpic}[width=0.9]Images/tau4-10_sigma3e-1_alpha_tau_c.pdf (b){overpic}[width=0.9]Images/tau45e-1-10_sigma35e-2_alpha_tau_c.pdf (c)

Figure 3: (Colour online) MSF for a fixed coupling strength (a) , (b) and (c) in the plane of the real part () and the delay time . The horizontal red lines show the position of the critical delay time . .

Based on the exemplary coupling strengths in Fig. 3 it is clear that the stability depends on both and . In all examples the vertical boundary line at , corresponding to the longitudinal eigenvalue (i.e. , where ), can be easily identified. It separates regimes of stable and unstable synchronisation. Another common characteristic is that the -dependent MSFs have a critical delay time at which the stable region splits into two separate, disconnected regions. For values above the stable region is found to the left of the longitudinal eigenvalue; whereas for values below there may be one stable region to the right of the longitudinal eigenvalue and one to the left. seems to be an important value, because it marks the most significant -dependent qualitative change in the MSF. Ultimately, the MSF can be divided into three regimes of : (i) two or more smaller separate regions of stability exist when ; (ii) one large region of stability exists when ; (iii) the MSF has one rotationally symmetric region of stability in the limit of (which holds already in good approximation if is of the order of the intrinsic oscillation period or larger).

Clearly, the topological criteria for stable synchronisation will be very different depending on whether the network is coupled with or . Observing the change in the limit cycles while varying offers the explanation for the qualitative change either side of in the MSF. For smaller values, the trajectory of the limit cycle is sensitive to changes in the coupling parameters, whereas, for larger values, the trajectory changes very gradually and approaches an elliptical form. The most predominant feature of the changing limit cycle is a large “kink” for values below . Figure 4(a) shows such a kink in the limit cycle induced by the coupling parameters and : upon bypassing the stable node the limit cycle deviates a little towards the unstable focus, located at the origin, before suddenly changing its trajectory back towards the saddle point. The kink has a large effect on the dynamics of synchronised nodes and is an important source of instability. Simulations have shown that it is at the kink that synchronised nodes begin to diverge from one another following any small perturbations. For larger the kink diminishes and disappears (see Fig. 4(b)).


[width=0.526]Images/tau535e-2_sigma35e-2_vectors_j00003400.pdf (a){overpic}[width=0.454]Images/tau6_sigma35e-2_vectors_j00002640_noy.pdf (b)

Figure 4: (Colour online) The position in phase space of one element representing the synchronised dynamics of a network (large green diamond) is shown relative to its delayed position (small green diamond). The black circle and the white circle represent the stable and saddle point, respectively. For small delay times such as (a) a kink may be the result of the relative positioning. For larger delay times such as (b) this does not occur. The movement of the phase point has a local component (dark green arrow) and a coupling component (light green arrow). The coupling component is directed towards the phase point at time (small green circle). , .

Let us explain the cause for the kink in greater detail. The large green diamond in Fig. 4(a) is the position of a node that represents the dynamics of a synchronised network in the phase space. The small green diamond is its delayed phase point, to which it is coupled, and represents the delayed dynamics of the same synchronised network (i.e. the small green diamond at time is the large green diamond at time ). A node coupled with a delayed node is attracted towards the direction of this delayed node (cf. Fig. 4(b)). As the node passes between the stable fixed point and the saddle point, it is slowed down, allowing its delayed node to “catch up” to it. For smaller this has the effect that the node is attracted in a clockwise direction towards its delayed node and gets slowed down (cf. Fig. 4(a)). This is the cause for the kink. Whereas for larger , the (small green) delayed node is always more than half a rotation behind, meaning that the (large green) node is attracted towards it in an anti-clockwise direction. This gives the limit cycle a form that is roughly circular as seen in Fig. 4(b).

As long as the kink is present in the limit cycle, a synchronised network will generally have very different topological conditions in order to be stable. Even when the limit cycle has a very predominant (sharp) kink, it will often be possible to find a class of networks that will still show stable synchronisation.

Note that in Fig. 3 one can see how the most stable part of the stable region (i.e. where is at its smallest) shifts away from as becomes smaller. Simulations have shown that this happens while the angle of the kink increases (making the kink “sharper”). A possible interpretation could be as follows. Each potential eigenvalue of the coupling matrix in a MSF is associated with one or more eigenvectors in the -dimensional space. These eigenvectors are all orthogonal to the vector that represents motion in the SM. This would suggest that the increasing angle of the kink corresponds to a shifting of the most stable direction of perturbation transverse to the SM.

5 Implications for complex networks with a small delay coupling time

An obvious observation is that networks with purely excitatory coupling, which are always stable for large delay times as mentioned at the end of Sec. 3, may not show stable synchronisation for . It was explained in Sec. 3 that increasing the probability of inhibitory links in the network can be a factor leading to unstable synchronisation. This occurs because a larger probability of inhibition can push a part of the eigenvalue spectrum of any network beyond the longitudinal eigenvalue at (because Gershgorin’s circle theorem no longer holds) and into the unstable region of the MSF. Now, for smaller , there may be a pocket of stability to the right of the longitudinal eigenvalue, so that increasing inhibition can make the otherwise unstable synchronisation of a network stable.

This means that the transitions between stable and unstable synchronisation as a function of the probability of inhibition discussed in Ref. LEH11 are no longer valid when is small. The transitions are now sensitive to the coupling parameters, not just the network parameters and . Due to the multiple regions of stability, eigenvalues may wander in and out of stable regions, while increasing the probability of inhibition . For large , increasing in the SW network model only results in one transition where the fraction of desynchronising networks (i.e. networks with unstable synchronisation) switches from 0 to 1 (cf. Fig. 2). For small , it is possible that jumps back to 0, before increasing again to 1. This will occur if there is a separate region of stability to the right side in the MSF that is large enough that all the eigenvalues lying over there can fit inside.

The observation of multiple transitions between stable and unstable synchronisation can be explained by looking at the eigenvalue spectra for SW networks for various values. As discussed above in relation to the MSF method, each network topology has a coupling matrix , the eigenvalues from which can be used to determine the stability of the network’s synchronised state. Because the “short-cuts” in the SW network are randomly introduced, each network with specific , and values may have many realisations. By calculating the eigenvalues for a large number of realisations of the SW network with certain parameters (i.e. given , and ), one can find the bounds for the eigenvalue spectra.

Figure 5 displays the superimposed eigenvalue spectra of 500 realisations for SW networks with parameters and (regular ring with excitatory coupling of nearest neighbours on either side) in dependence on the probability of inhibition . The longitudinal eigenvalues located at have been removed. One can see that the bounds for possible eigenvalues shift depending on the value of . One can also see how the spectra begin to increasingly resemble the semicircular distribution ALB02a of a random network for larger values, where the networks have lost their SW properties.

Figure 5: Spectra of eigenvalues of undirected small world networks with and and varying probability of inhibitory links. 500 realisations of the eigenvalue spectrum are plotted for each value of .

The eigenvalue spectra bridge the gap between observations of the MSF (i.e. the dependence of the stability of synchronisation on the eigenvalues) and what is seen in these transitions (i.e. the dependence of stability on the network topology). Increasing allows isolated eigenvalues to increase in value and, in terms of the MSF, shift their locations further to the right in the -plane. This is visualised in Fig. 6. Figures 6(b), (d) and (f) show the eigenvalue spectra for SW networks with elements and , 20, and 50, respectively.


[width=0.33]Images/trans_N200_k40_small.jpg (a){overpic}[width=0.33]Images/trans_N200_k20_small.jpg (c){overpic}[width=0.33]Images/trans_N200_k50_small.jpg (e){overpic}[width=0.33]Images/eigenvalues_N200_k40_small_p0-5e-1_tau65e-1_sigma3e-1.jpg (b){overpic}[width=0.33]Images/eigenvalues_N200_k20_small_p0-5e-1_tau65e-1_sigma3e-1.jpg (d){overpic}[width=0.33]Images/eigenvalues_N200_k50_small_p0-5e-1_tau65e-1_sigma3e-1.jpg (f)

Figure 6: (Colour online) Fraction of desynchronized networks in dependence of the probability of additional inhibitory links for 500 realisations of networks of with (a) , (c) , and (e) with and . Corresponding eigenvalue spectra for (b) , (d) , and (f) with realisations for each value. Here the green shaded regions represent the stable regions of the real part of the MSF. .

Figures 6(a), (c) and (e) show the corresponding fraction of desynchronised networks in dependence on the probability of additional inhibitory links . Consider, for instance, SW networks with parameters and . In Fig. 6(a) is shown for the exemplary coupling parameters and with a corresponding plot of eigenvalues in Fig. 6(b). Note that the grey (green) shaded regions represent the stable regions of the line on the -plane of the MSF for and (cf. Fig. 1(c)). For coupling parameters where this region of stability is not large enough, may briefly dip below without decreasing to , because only some realisations may have eigenvalues that fit inside the stable region. In Fig. 6(c) where , dips down to , because at most of the realisations show stable synchronisation (i.e. all the eigenvalues are inside the stable region). If the distance between the larger isolated eigenvalues matches the distance between stable regions (which is almost the case in Fig. 6(d)), then the transition curve only just touches the axis at some value of before increasing back to . Figure 6(e) is an example where a further transition is possible, because not only do all eigenvalues begin in stable regions at , but there is another regime of where all eigenvalues fit inside stable regions. Note that the length of the transition from to is actually a measure of the variance of the isolated eigenvalues for an ensemble of realisations. For example, the variance of the isolated eigenvalues decreases as is increased, so that the length of transition will be shorter in larger networks, and the transition becomes sharper.

When constructing a histogram of the eigenspectra for a particular value of – see Fig 7(a) – the larger isolated eigenvalues seen to the right, e.g. in Fig. 6(d), result in small peaks of eigenvalues.


[width=0.49]Images/eig_spectrum_N100_k10_2e-1_small.pdf (a){overpic}[width=0.49]Images/eig_spectrum_N100_k10_2e-1_random.pdf (b)

Figure 7: Histograms of eigenvalue spectra for realisations of (a) SW networks and (b) random networks with the same number of nodes, the same number of excitatory links and the same probability of inhibitory links. , , and number of bins is .

For each realisation of many SW networks there are isolated eigenvalues that are larger than most other eigenvalues in the spectrum. The small peaks of eigenvalues result from perturbations in isolated eigenvalues. The spectrum of a regular ring network (i.e. a SW network with ) can be found analytically using the graph’s symmetry operations FAR01 and is given by


where . Figure 8 shows this eigenspectrum for the case.

Figure 8: The eigenvalue spectrum for a regular network with and . is the value of the eigenvalue.

As is increased these eigenvalues will be slightly perturbed (in a random manner) by the changing network structure, of which many realisations are possible. In Fig. 8 one can see the isolated eigenvalues to the right-hand side of the plot that eventually evolve into the smaller side peaks of eigenvalues in histograms for larger (see Fig 7(a)). These peaks show the distribution of the eigenvalue under the influence of the random nature of the SW “short-cuts” creation process – which is why the smaller eigenvalue peaks can be approximated by a normal distribution 111Actually, because the multiplicity of the above mentioned isolated eigenvalues at is , the smaller eigenvalue peaks are two normal distributions that, at least for small values, overlap each other to a large extent. This is furthermore the reason why the desynchronisation transitions observed in Ref. LEH11 and Figs. 6(a), (c) and (e) look like the cumulative (integrated) distribution function. They reflect how increasing brings the area of the small eigenvalue peak accumulatively into the unstable region of the MSF.

To explain why the peaks wander with increasing , consider not normalising the rows of matrix , so that the row sums are not necessarily equal to . Then the location of the longitudinal eigenvalue decreases with , because it is equal to the average row sum of which is equal to . The locations of the other peaks increase slightly to maintain the eigenvalue sum of zero (a result of the trace of being zero, since there is no self-feedback coupling), although for large this has only a small effect. Thus, scaling the eigenvalues so that the longitudinal eigenvalue is always at means that the transversal eigenvalues are multiplied by , so that they appear to increase with , as seen in Figs. 6(b), (d) and (f).

The case of random Erdős-Rényi networks ERD59 is different. At there are no larger gaps between the eigenvalues, such as those for SW networks at . If a random network is fully connected (i.e. is one single component) then the distribution of its eigenvalues at are confined to a semi-circle, given by the semicircular distribution


where is the probability of excitatory links ALB02a . While there are initially no gaps in the eigenvalue distribution at , adding inhibitory links to the network (i.e. increasing ) does not lead to the creation of any gaps like those seen for SW networks. To illustrate this point, Fig. 7(b) shows the eigenvalue distribution for many realisations of Erdős-Rényi random networks with the same number of nodes, as well as the same number of excitatory links and the same probability of inhibitory links. As a result, multiple transitions between stable and unstable synchronisation will not occur.

One can, however, use Eq. (8) to find parameters and such that the eigenvalue distribution is confined to a stable region of the MSF at , as long as there is a stable region centred around (for example, see Fig. 1(c)). This allows for a desynchronisation transition with increasing , such as in Fig. 2.

6 Small-world networks of Stuart-Landau oscillators

The appearance of multiple transitions between synchronisation and desynchronisation is based on the interplay of the peculiar spectral properties of the small-world topology and the unconnected stable regions that show up in the MSF for the excitable SNIPER model for small delay times.

In this Section we show that such unconnected stable regions may also appear in the MSF of a purely oscillatory system. We consider the Stuart-Landau oscillator, which is the normal form of any oscillatory system near a supercritical Hopf bifurcation. The dynamics of the Stuart-Landau oscillator is governed by a complex variable , such that we have in Eq. (4). The local dynamics of the Stuart-Landau oscillator is given by


For , the uncoupled system exhibits self-sustained limit cycle oscillations with radius and frequency . The coupled system shows – depending on coupling parameters and topology – a variety of synchronised and desynchronised solutions, including amplitude death and cluster synchronisation CHO09 ; CHO11 .

Here we focus on complete (zero-lag) synchronization.


[width=0.46]sachin/fig1_1.pdf (a){overpic}[width=0.46]Images/trans_N200_k40_small_STU.jpg (c){overpic}[width=0.46]sachin/fig3_1a.pdf (b){overpic}[width=0.46]Images/eigenvalues_N200_k40_small_p0-5e-1_STU_tau15e-1pi_sigma8e-2 (d)

Figure 9: (Colour online) Network of Stuart-Landau oscillators according to Eqs. (4) and (9): (a) MSF in the plane for . (b) MSF in the plane for . (c) Fraction of desynchronized networks in dependence of the probability of additional inhibitory links for 500 realisations of networks of and , . (d) Corresponding eigenvalue spectra, where the green shaded areas represent stable regions. Other parameters: , , .

Fig. 9(a) shows the MSF for sets of parameters that allow for unconnected stable regions in the complex plane of the eigenvalues . As shown in Fig. 9(b) a connected stability region occurs for integer multiples of , while unconnected stability regions – which can lead to multiple phase transitions in SW networks – occur in between.

For the parameters used in Fig. 9(a), Fig. 9(c) shows the fraction of desynchronised networks for the modified small-world model as in Sec. 5. Figure 9(d) shows the eigenspectra for networks of elements, starting from a one-dimensional ring with excitatory coupling to neighbours to either side. Increasing the probability of additional inhibitory links broadens the spectra and leads to a second regime of stable synchronisation matching the transitions in Fig. 9(c).

7 Conclusion

We have investigated transitions between synchronisation and desynchronisation in complex networks of delay-coupled excitable elements of type I, induced by varying the balance between excitatory and inhibitory couplings in a small-world topology. In our analysis we have used the master stability function approach. For large delay times it seems that both type-I neurons, as modelled here, and type-II neurons, as modelled in Ref. LEH11 , must fulfill similar topological conditions in the network to allow for a stable synchronised state. This is different when considering small delay times. For a range of small coupling strengths and small delay times we have found novel multiple transitions between synchronisation and desynchronisation, when the fraction of inhibitory links is increased. Unlike the Erdős-Rényi random network, a small world model for complex networks with regular excitatory couplings and random inhibitory shortcuts has eigenvalue spectra with gaps between the larger eigenvalues, so that histograms of many realisations reveal isolated peaks of possible eigenvalues. It was shown that, because of this, synchronised small world networks can go through multiple transitions of stability in dependence on the probability of inhibitory short-cuts. Our results are valid beyond the special model of type-I excitability, as we have demonstrated using the Stuart-Landau oscillator.

This work was supported by the DFG in the framework of the SFB 910. PH acknowledges support by the BMBF (grant no. 01GQ1001B).


  • (1) R. Albert, A.L. Barabási, Rev. Mod. Phys. 74, 47 (2002)
  • (2) M.E.J. Newman, SIAM Review 45, 167 (2003)
  • (3) S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, D.U. Hwang, Physics Reports 424, 175 (2006)
  • (4) M. Timme, F. Wolf, T. Geisel, Phys. Rev. Lett. 89, 258701 (2002)
  • (5) A.J. Steele, M. Tinsley, K. Showalter, Chaos 16, 015110 (2006)
  • (6) P. Ashwin, G. Orosz, J. Wordsworth, S. Townley, SIAM J. Appl. Dyn. Syst. 6, 728 (2007)
  • (7) V. Flunkert, S. Yanchuk, T. Dahms, E. Schöll, Phys. Rev. Lett. 105, 254101 (2010)
  • (8) J. Lehnert, T. Dahms, P. Hövel, E. Schöll, Europhys. Lett. 96, 60013 (2011)
  • (9) A. Englert, S. Heiligenthal, W. Kinzel, I. Kanter, Phys. Rev. E 83, 046222 (2011)
  • (10) I. Kanter, M. Zigzag, A. Englert, F. Geissler, W. Kinzel, Europhys. Lett. 93, 60003 (2011)
  • (11) I. Omelchenko, Y.L. Maistrenko, P. Hövel, E. Schöll, Phys. Rev. Lett. 106, 234102 (2011)
  • (12) T. Dahms, J. Lehnert, E. Schöll, Phys. Rev. E 86, 016202 (2012)
  • (13) A. Hagerstrom, T.E. Murphy, R. Roy, P. Hövel, I. Omelchenko, E. Schöll, Nature Physics 8, 658 (2012), published online
  • (14) T. Nepusz, T. Vicsek, Nature Physics 8, 568 (2012)
  • (15) D.H. Perkel, T.H. Bullock, Neurosci Res Program Bull 6, 221 (1968)
  • (16) D.J. Watts, S.H. Strogatz, Nature 393, 440 (1998)
  • (17) L.A. Adamic, The Small World Web, Vol. 1696/1999 of Lecture Notes in Computer Science (Springer Berlin / Heidelberg, 1999), ISBN 978-3-540-66558-8
  • (18) O. Sporns, G. Tononi, G.M. Edelman, Cereb. Cortex 10, 127 (2000)
  • (19) O. Sporns, Biosystems 85, 55 (2006)
  • (20) V. Latora, M. Marchiori, Phys. Rev. Lett. 87, 198701 (2001)
  • (21) R.D. Traub, R.K. Wong, Science 216, 745 (1982)
  • (22) P. Roelfsema, A. Engel, P. K nig, W. Singer, Nature 385, 157 (1997)
  • (23) A. Engel, P. Fries, W. Singer, Nature Reviews Neuroscience 2, 704 (2001)
  • (24) A.S. Pikovsky, M.G. Rosenblum, J. Kurths, Synchronization, A Universal Concept in Nonlinear Sciences (Cambridge University Press, Cambridge, 2001)
  • (25) L. Melloni, C. Molina, M. Pena, D. Torres, W. Singer, E. Rodriguez, J. Neurosci. 27, 2858 (2007)
  • (26) B. Haider, A. Duque, A.R. Hasenstaub, D.A. McCormick, J. Neurosci. 26, 4535 (2006)
  • (27) M.E.J. Newman, D.J. Watts, Phys. Lett. A 263, 341 (1999)
  • (28) E.M. Izhikevich, Int. J. Bifurcation Chaos 10, 1171 (2000)
  • (29) A.L. Hodgkin, J. Physiol. 107, 165 (1948)
  • (30) G. Hu, T. Ditzinger, C.Z. Ning, H. Haken, Phys. Rev. Lett. 71, 807 (1993)
  • (31) J. Hizanidis, R. Aust, E. Schöll, Int. J. Bifur. Chaos 18, 1759 (2008)
  • (32) R. Aust, P. Hövel, J. Hizanidis, E. Schöll, Eur. Phys. J. ST 187, 77 (2010)
  • (33) P. Erdős, A. Rényi, Publ. Math. Debrecen 6, 290 (1959)
  • (34) R. Monasson, Eur. Phys. J. B 12, 555 (1999)
  • (35) L.M. Pecora, T.L. Carroll, Phys. Rev. Lett. 80, 2109 (1998)
  • (36) S.A. Gerschgorin, Izv. Akad. Nauk. SSSR 7, 749 (1931)
  • (37) I. Farkas, I. Derenyi, A.L. Barabási, T. Vicsek, Phys. Rev. E 64, 026704 (2001)
  • (38) C.U. Choe, T. Dahms, P. Hövel, E. Schöll, Phys. Rev. E 81, 025205(R) (2010)
  • (39) C.U. Choe, T. Dahms, P. Hövel, E. Schöll, Control of synchrony by delay coupling in complex networks, in Proceedings of the Eighth AIMS International Conference on Dynamical Systems, Differential Equations and Applications (American Institute of Mathematical Sciences, Springfield, MO, USA, 2011), pp. 292–301, DCDS Supplement Sept. 2011
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description