Effect of small-world topology on wave propagation on networks of excitable elements

We study excitation waves on a Newman-Watts small-world network model of coupled excitable elements. Depending on the global coupling strength, we find differing resilience to the added long-range links and different mechanisms of propagation failure. For high coupling strengths, we show agreement between the network and a reaction-diffusion model with additional mean-field term. Employing this approximation, we are able to estimate the critical density of long-range links for propagation failure.


Introduction
Excitable media are well studied model systems in a variety of applications ranging from chemical [1] to electronic systems [2] and lasers [3] and from heart-muscle tissue [4] to neural systems [5,6]. An excitable system rests in a stable steady state, but after a sufficiently strong perturbation performs a long excursion in phase space, i.e., emits a spike, before returning to the stable steady state again. Excitable media arise when excitable elements are coupled spatially.
The spatial coupling facilitates an abundance of dynamical behavior amongst which are Turing patterns [7], traveling waves [1], spots [8] and spiral waves [9], to name just a few. Especially wave-like spatio-temporal behavior is interesting from a neuroscience point of view. It is observed in living neural tissue and considered to play a role in neural information processing in different tasks [10,11]. Traveling waves are a generic phenomenon in cortical dynamics [12], and have successfully been used to describe features of cortical spreading depolarization [13,14,15].
In recent years dynamical systems coupled in complex network architectures have attracted a lot of attention [16,17,18,19,20]. Those systems can also occur in a wide variety of applications ranging from power grids [21] to biological networks [22]. In neuroscience, scenarios with excitable elements coupled in a chain-like one-dimensional topology have been suggested as a mechanism for the occurrence of these waves [11]. In related works, a one-dimensional network model of pyramidal cells and interneurons produces saltatory propagation, and excitatory connections play a crucial role [12].
The topology of a network is a key factor influencing the dynamic behavior of the system. Network topologies can range from well-ordered lattice systems, which resemble spatially extended systems, to random topologies where the notion of space loses its meaning. A very interesting intermediate form are the so-called small-world topologies [23,24], in which strong local connectivity is combined with a few longrange connections enabling short mean path lengths.
Small-world topologies are of interest for the description of anatomic and functional brain networks [25]. They are considered a powerful and versatile approach to the structure of those systems, and it is argued that within the cerebral architecture, they are preferred as a trade-off between network efficiency and wiring cost [26].
On the dynamical side, wave and front propagation have been studied in network topologies such as trees, where fronts can be pinned and waves cease to propagate [27,28]. Small-world topologies have been shown to make systems support sustained activity [29,30].
Regarding the influence of topology on the behavior of traveling waves, it has been shown that short-range connections mediate traveling waves in phase oscillator models which are used as simple models for cortical waves [31]. Also, additional longrange connections in a chain of locally coupled oscillators can be used to generate traveling waves of different wavelengths [32].
A change of topology may occur in pathological states as, e.g., multiple sclerosis, which affects the most expensive (long-range) links in the brain network [26]. Moreover, propagation failure of waves is an important aspect in different areas of physiology [33,34].
To the best of our knowledge, the combination of these relevant ingredients, i.e., generic excitable elements, small-world topology and propagation failure, has not been addressed in a satisfactory way, except for very few exemptions, e.g. [35]. We attempt to do so by using a generic model of excitability, the well-known FitzHugh-Nagumo model [36,37], combined with a Newman-Watts small-world architecture [17] as well as techniques from spatially continuous systems to study the behavior of excitation waves.
The structure of the paper is as follows: In section 2, we introduce the dynamics and the network model. We draw a connection between a ring network and a onedimensional continuous reaction-diffusion system. We discuss traveling wave solutions in both systems, in particular their spectral and stability properties, and we elaborate on the difference between the dynamics of the network and the reaction-diffusion system. In section 3, we consider a small-world topology. We highlight the different mechanisms leading to propagation failure in different dynamical regimes. In section 4, we modify the continuous reaction-diffusion model in order to incorporate the effect of the small-world topology. We give an estimate of the critical link density at which excitation waves cease to exist in the small-world network. In section 5, we summarize our findings.

Dynamics
As a generic model of excitable dynamics we use the FitzHugh-Nagumo model [36,37] on an undirected, unweighted network, where neighboring nodes are coupled by the difference in the activator concentrations. Thus, the model readṡ where u i is the activator, and v i the inhibitor at node i, and A ij is the adjacency matrix of the network, κ > 0 is the coupling strength, ε ( 1) and β are the timescale separation and the excitation threshold of the local dynamics, respectively. The local dynamics of Eq. (1) (i.e., without the coupling term) possesses a steady state at u * = −β, v * = −β + β 3 /3. For a value of |β| > 1, this steady state is stable. It undergoes a supercritical Hopf bifurcation at |β| = 1. For the network Eq. (1) a linear stability analysis shows that the eigenvalues of the linearization around the homogeneous steady state u * i = u * , v * i = v * are given by i are the eigenvalues of the local dynamics. Comparing terms of µ ± j , we note that the real part of the square root term is always smaller than the absolute value of the term before the square root. So, when the first term is negative, µ ± j is also negative. This is always the case if u * > 1 and thus if the steady state of the local dynamics is stable. On the other hand, if the first term is positive, µ + j > 0. A sufficient condition for this is u * < 1 and λ j = 0, which, as there always exists a λ j = 0, always happens if the steady state of the local dynamics is unstable.
Thus, we conclude that the homogeneous steady state inherits the stability of the local equations and thus when u * , v * is stable, so is the homogeneous steady state For ε 1, the local dynamics of Eq. (1) is a prototypical example for a slow-fast system. In that case, the periodic orbit which emerges from the Hopf bifurcations at |β| = 1 shows a canard explosion [38] immediately after the bifurcation point (|β| 1), evolving to the large-amplitude limit cycle of a relaxation oscillator.
Our interest is focused on the regime where |β| > 1 and no limit cycle (stable or unstable) exists in the local dynamics. Then the local dynamics Eq. (1) serves as a paradigmatic example for type-II excitability. Systems of type-II excitability do not possess a constant threshold that separates stimulations leading to a spike from those that do not. They are rather characterized by a "threshold trajectory", around which the system is very sensitive to the size of the stimulus. If we view the local dynamics of Eq. (1) as a singularly perturbed system, we note that by Fenichel theory [39], the normally hyperbolic parts of the u-nullcline (which give the critical manifold of the system in the limit ε → 0) are perturbed to attracting and repelling slow manifolds S a ε , S r ε , respectively. Any trajectory that passes near S r ε is repelled from it with v staying almost constant (because ε 1) until it hits S a ε . A trajectory that follows the repelling slow manifold S r ε for a 'considerable amount of time' is called a canard trajectory [38]. In the case of the FitzHugh-Nagumo dynamics it is a common choice to take such a canard trajectory as the threshold trajectory, e.g., one uses the trajectory that goes through the point where the repelling and attracting parts of the critical manifold meet (which is given by the u-nullcline), see Fig. 1 In the remainder of this work, we will fix the parameters at ε = 0.04 and β = 1.1, such that the model is in the excitable regime.

Network
The network topology for our model is chosen to be a ring topology consisting of N nodes, where each node is coupled to its k neighbors to the left and its k neighbors to the right, so that every node has degree 2k. We perturb the topology by adding a certain number n of long-range links, where both ends of these links are chosen randomly [17]. This is a small modification of the well-known Watts-Strogatz smallworld model [24], where the long-range links are replacing links of the regular ring network. For a certain range of n the Newman-Watts model shows small-world behavior as well, i.e., short average path length and high clustering coefficient, see Fig. 1(a).

Connection to one-dimensional reaction-diffusion systems
In general, a dynamic system on a ring network with coupling range k has a well defined continuum limit. Consider the equation for such a system (i=1,...,N) If we define a continuous spatial variable x and a distance h ≡ 1 √ κ between two adjacent nodes on the ring and assume that there is a function u(t, x) such that u i (t) ≡ u(t, i √ κ ), Eq. (2) can be expressed as Eq.
By rescaling κ → q(k)κ, this simplifies to where L is now given by L = N √ q(k)κ . Because in the second line of Eq. (1) no coupling term appears, Eq. Eq. (1) on a ring network simplifies in the limit (4) to Comparing the discrete ring network with this one-dimensional reaction-diffusion equation, we note that the three parameters that determine the coupling and the topology of the ring network N, k, and κ are translated into one parameter L in the continuum limit, while the number of parameters of the local dynamics does not change.
Traveling-wave solutions With the parameter values chosen here (β = 1.1, ε = 0.04), the local dynamics of Eq. (8) posesses exactly one stable steady state which becomes a stable homogeneous steady state in the spatially extended system. Apart from this solution, Eq. (8) is known to support traveling-wave solutions, stable 'fast waves' and unstable 'slow waves', see e.g. [40,41,42].
Varying L in Eq. (8) and monitoring the propagation velocity c of a traveling-wave solution that exhibits exactly one peak (a pulse), we obtain the dispersion relation of Eq. (8). The spectrum of the linear stability analysis around such a traveling-pulse solution is very closely connected to the essential spectrum around a wave-train in an infinitely extended spatial domain subject to the same differential equation. With the method proposed in [43], we calculate the essential spectrum of the linearization of Eq. (8) around the connected wave train. The essential spectrum is continuous and contains every point of the spectrum of a traveling wave solution of Eq. (8). These points additionally fulfil a certain boundary condition. The method we use to calculate the spectrum is elaborated in more detail in [43].
Starting with stable traveling-wave solutions and large L, the calculated spectrum tell us that a destabilization of the wave occurs at a critical value of L cr ≈ 30.756 by two complex conjugate eigenvalues with nonzero imaginary part. However this point does not yet mark the lowest possible propagation speed c. This is attained by the branch of unstable solutions after lowering L further to a value of L ≈ 19.6. After that, the branch of unstable solutions continues in form of a fold bifurcation in direction of rising L with the speed converging against c ≈ 0.493. On this branch, L can increase without bounds, and this branch is associated with the transition to an (unstable) solitary traveling pulse. Also on the stable branch, L can increase without bounds but the solutions do not lose stability and converge to a stable propagating solitary pulse.
After the destabilization at L cr , the spectrum changes in a complex way, leaving in the end an isolated eigenvalue as the only object in the right half-plane. The dispersion relation (propagation velocity c vs. domain size L) including the leading parts of the (essential) spectrum for selected points (c)-(k) is illustrated in Fig. 2(a). In the pictures of the spectrum, the essential spectrum is marked by continuous lines and the values of the essential spectrum that occur for periodic boundary conditions are marked with asterisks. The latter are the relevant ones for our model. The order of occurrence of the points (c)-(k) is the following, (c) All eigenvalues are in the left halfplane (except for the one that is always present at zero corresponding to the Goldstone mode of translation invariance), (d) two complex conjugate eigenvalues crossing the imaginary axis at L cr , (e) a second pair of complex conjugate eigenvalues crossing the imaginary axis, (f),(g) from the leading part of the essential spectrum, a circle comprising the Goldstone mode eigenvalue is forming and detaching, (h) the first two eigenvalues that have crossed the imaginary axis merge on the real axis and split in different directions, (i) one of the eigenvalues that has merged on the real axis crosses zero, (j) the second two eigenvalues that have crossed the imaginary axis cross the imaginary axis again in the opposite direction. (k) One eigenvalue is left in the right half-plane, the rest of the eigenvalues are in the left half-plane (again except for the Goldstone mode eigenvalue at zero). Thus after the first instability at L cr , a scenario with several secondary instabilities evolves.
The reaction-diffusion system Eq. (8) is a good approximation for the ring network Eq. (7) for large N and large κ.
On the ring network, the above sequence of changes in the spectrum would translate to the following sequence of bifurcations: A torus bifurcation through which the wave loses its stability at (d), another torus bifurcation at (e), a saddle-node bifurcation at (i), a torus bifurcation at (j). The minimum domain size L cr for which a traveling-wave solution is stable translates to a maximum coupling strength according to Eq. (4). We denote this maximum coupling strength by κ high = N 2 q(k)L 2 cr . In Eq. (8), we can make L arbitrarily large (corresponding to decreasing κ in Eq. (7)). The eigenvalue in the right half-plane will stay in the right half-plane and thus the solution will remain unstable. Also on the stable branch of wave solutions of Eq. (8), the domain length can increase arbitrarily and the solutions remain stable.
However, in Eq. (7) the regime of high L is equivalent to the regime of low coupling strength, and thus the discreteness of Eq. (7) becomes important. The profile versus time that a single node shows differs appreciably from that of a solution of Eq. (8) at a fixed spatial point.
Also, when the wave propagates in a discrete system, as Eq. (7), it cannot propagate at arbitrarily low coupling strengths [44,45]. The critical coupling strength κ at which the wave ceases to propagate has been approximately calculated in [45]. We will point out what happens to the dispersion relation at κ low . To the best of our knowledge this has not yet been investigated.
We calculate the dispersion relation from Eq. (7) by choosing a certain number of nodes and performing a numerical continuation for the resulting full system of coupled ordinary differential equations, using AUTO-07p [46]. Due to the closed ring topology, a traveling-wave solution manifests itself as a periodic solution. The propagation speed c can easily be calculated as c = , where T is the period of the solution. Using these parameters c and L instead of κ and T , we can compare the dispersion relation of Eq. (7) with that of Eq. (8) in the regime of high coupling strength κ or small (virtual) domain size L, respectively.
For high coupling strengths we find the same behavior as expected by the dispersion relation of Eq. (8). However, for decreasing coupling strength, the propagation speed of the wave solutions does not stay constant, as it does in the continuous system. For the unstable 'slow-waves' it increases whereas it decreases for the stable 'fast-waves' until both branches meet at the critical low coupling strength κ low in a saddle-node bifurcation. The value of κ low has been approximately calculated in [45]. However, the coalescence of the branch of stable propagating waves and that of the unstable ones in a saddle-node bifurcation at the destabilization point has to our knowledge not yet been reported. Because of this saddle-node bifurcation, the dispersion relation of wave solutions on the discrete ring Eq. (7) is given by a closed  The lower boundary of the coupling strength κ low is a genuine effect of the discreteness of Eq. (7) and does not depend on the number of nodes N (assuming that this is large enough). However it does depend on k. A comparison of the dispersion relations for different network sizes N and nearest neighbor numbers k is shown in fig. 3.

Setup
We are interested in the behavior of traveling waves on the ring, when the network topology is perturbed by adding links randomly as described in Sec. 2.2.
To this end, we numerically simulate a traveling-wave solution of Eq. (1) on a ring network without additional links. At a certain point of time we instantly add a number n of links at random positions to the network and monitor the resulting behavior, i.e., sustained wave propagation or decay to the homogeneous rest state.
We do this for varying coupling strength κ ∈ (κ low , κ high ) and for rings of different sizes N and different nearest neighbor numbers k. The number n of random links is changed from 1 to kN and for every combination (N, k, n), 200 realizations are investigated. For the numerics we use a Runge-Kutta-Fehlberg method with adaptive timestep for fast simulation and reliability.
The system either shows sustained activity in form of a perturbed traveling wave, or it decays into the stable homogeneous steady state. We simulate for 6000 time units, assuming that the wave is stable if no decay to the stable homogeneous steady state has occurred within this time ‡.

Numerical observations
Whether the wave decays to the stable homogeneous state or not, does not depend exclusively on either κ, n or even the particular realization of the small-world network.
For example there are different realizations with the same n that for the same coupling strength κ either support sustained wave activity or do not, see for example Fig. 5. On the other hand, we can also find one realization that for one value of κ supports ongoing activity and for a slightly different κ does not, see for example Fig. 4(a) and (b), where the particular networks are the same, leading to propagation suppression in 4(b) and to sustained activity in 4(a).
For a systematic investigation, we examine the fraction f sust (N, k, n, κ) of realizations that support an ongoing activity. We have performed one investigation for several different network sizes ranging from N = 150, to N = 2000, 3 values of k (k = 1, 2, 3), about 40 values of n (distributed logarithmically between 1 and N k) and 200-1000 values of the coupling strength κ (distributed logarithmically between κ low and κ high . ‡ Taking all solutions that actually decay within this timespan, the median time until decay is typically around 10 and always below 21. The 0.99-quantile is typically around 600 with one outlier (N = 1000, k = 4) where the 0.99-quantile is 3033. Thus it is safe to assume that almost all solutions that have "survived" the 6000 time units threshold will survive further. Generally we observe that the larger the number of additional links n, the more realizations do not support sustained wave activity. Finally, if n becomes too large, no realizations with ongoing activity can be found anymore. When κ is located well within the range [κ low , κ high ], we find that the transition takes place within a small range of n N k , see Fig. 6. However, the point of transition varies considerably depending on (i) the network size N , (ii) the nearest neighbor number k, and (iii) the coupling strength κ. Note that these parameters can also be chosen such that already for n = 1, f sust is considerably below 1.0 (e.g. N = 250, k = 2, κ = 7.0 in Fig. 6(b), leftmost point for those parameters).
In dependence on the coupling strength κ, we find two limiting regimes. When κ approaches either κ low or κ high , f sust is approaching zero. This behavior is expected from the dispersion relation discussed in Sec. 2, as this shows no possible travelingwave solutions beyond these values even on the unperturbed ring. However, as the propagation of the unperturbed wave differs in the two regimes, so does the sustained propagation or the propagation failure on the small-world network. As can be seen in Fig. 6(a), when the coupling strength is low, the network size N does not play a role and networks with the same k show the same transition point in n. When the coupling strength is high, as in Fig. 6(b), no effect like this is noticeable, but in Sec. 4.1 we will derive transformed coordinates in which the transition points fall together (approximately).
Low κ -discrete regime (Fig. 7) When approaching the lower limit value κ low , the fraction of realizations that support sustained wave activity f sust becomes smaller for all numbers of added links n > 0. At κ = κ low it is zero for all n > 0.
In this regime of low κ, the propagation of the excitation takes place in a saltatory fashion. By this, we mean that the excitation 'hops' from node to node. More precisely, the state of one node reaches maximum excitation before the next node In the transition that we observe we find three values κ low < κ 1 < κ 2 at which the effect which added links have changes suddenly. These values are independent of N , but they depend on k.
If κ < κ low , we find no stable traveling wave solutions whatsoever, even for n = 0. This is expected from the dispersion relation (see. Fig. 3).
If κ low < κ < κ 1 , one additional link will lead to propagation failure once the traveling wave reaches one of the nodes which this link joins. When the node that is about to become excited has one end of the additional link, the other end will point to a node that is in the rest state. This is because the propagation is saltatory, there is only one node that is excited at one instance of time. In this range of κ, this link is sufficient to prevent this node from becoming excited, and thus the propagation is quenched. As a consequence, f sust (n > 0) = 0 but f sust (n = 0) = 1 in this range of κ. This behavior is illustrated in Fig. 7(a). There the exemplary network has one additional link from node 246 to 406. Propagation suppression by coupling back to an unexcited node happens as soon as the wave reaches node 406.
If κ 1 < κ < κ 2 , one additional link may lead to propagation failure, but it does not necessarily (Fig. 7(b) In this range of κ, a traveling-wave solution can pass one end of the additional link without being suppressed. When the node at the first end of the additional link becomes excited, this excitation is also coupled to the node at the remote end. Because of the coupling scheme of the ring, this node is coupled to more nodes in the rest state than the successor node to the node at the first end of the link. As a consequence, the node at the first end will be able to trigger a full excitation in its successor node but not in the node at the remote end. Here a subthreshold excitation is generated, which does not propagate further. However, this k κ low κ 1 κ 2 1 ≈ 0.0324 ≈ 0.0339 ≈ 0.0359 2 ≈ 0.0233 ≈ 0.0235 ≈ 0.0481 3 ≈ 0.0169 ≈ 0.0170 ≈ 0.5890 Table 1: Approximate transition values κ low , κ 1 and κ 2 for different nearest neighbor numbers k in the discrete limit.
sub-threshold excitation leads to an increased inhibitor level which takes a certain time to decay back to the steady state value. If the remote end of the shortcut will be reached by the traveling wave before the inhibitor level has sufficiently decayed, the propagation will stop. This behavior leads to decrease of f sust for increasing n in the mentioned range of κ, because as n increases, the more likely it becomes that one of the additional links spans a short enough distance. An exemplary time series for this behavior is shown in Fig. 7(b). There the exemplary network has two additional links from node 60 to 176 and from 343 to 347. Propagation failure happens by the raised inhibitor level at node 347 due to previous sub-threshold excitation mediated by the additional link.
If κ > κ 2 , a traveling wave solution can still pass one end of the additional link without becoming suppressed. But now a full excitation, leading to a pair of traveling waves with opposite propagation directions will be generated at the remote end of the additional link. As traveling waves are only generated pairwise and the annihilation also takes place in pairs, the mechanism for propagation failure must work differently. One possible mechanism that we found and that can also be seen in Fig. 7(c) is that two or more additional links end very close to each other (or even on the same node). Then, as can be seen in fig. 7(c) it might (i) not be possible to excite a secondary wave pair here and also (ii) propagation can stop here due to the strong coupling back to the rest state as in the regime below κ 1 . Generally, as f sust is close to unity for small n in this regime, propagation failure seems to be mostly caused by more complex mechanisms which need several additional links. In the exemplary timeseries for Fig. 7(c) there are four additional links from node 90 to 127, 38 to 142, 142 to 221 and 43 to 326. The propagation stops at node 142, where two additional links end. A secondary wave pair is never excited at this node because the other additional link couples back this node to the rest state.
To summarize: For κ < κ low , no traveling wave solutions exist for any n. For κ low < κ < κ 1 no stable traveling wave solutions exist for n ≥ 1. For κ 1 < κ < κ 2 f sust decreases with increasing n, no secondary waves can be excited. For κ > κ 2 the excitation of secondary waves is possible, leading to a sudden rise of f sust for intermediate n.
The approximate observed values for κ low , κ 1 and κ 2 are given in Table 1. We note that the values for κ low and κ 1 (i) decrease with increasing coupling range k as does (ii) the distance between the two. This is expected as (i) increasing k has a similar effect as raising κ on the ring and (ii) as the long-range links have the same weight as the local links, they have less impact, if k becomes larger. κ 2 on the contrary increases with increasing k. This is also expected for the same reason as (ii). For k = 3, κ 2 is not even located in the regime of saltatory propagation anymore. High κ -continuum regime (Fig. 8) If we are in the regime of high coupling strength κ, f sust decreases with increasing κ until κ reaches a maximum value κ high above which no traveling wave solutions are found even for n = 0, see Fig. 8(a). The value of κ high is related to viewing the ring network as the discretization of a reaction-diffusion system. For that there exists a minimum domain length L cr which corresponds directly to κ high , see Sec. derivation and discussion. This is shown in Fig. 8(b), where the fraction of sustained wave activity f sust is plotted versus n/L for various N and k. The coupling strength κ in this plot is adapted for each (N, k) so that L = 200 is constant in that figure. Note that the transition points from sustained wave activity to propagation suppression coincide better, the lower k and the larger N . For higher k (k = 2, 3) the networks with smaller N need to have (much) smaller κ to have the same L and thus are not well located in the continuum regime anymore. Thus in order to show that the approximation works very well in that case, networks with N > 2000 would need to be simulated, which would have been numerically too expensive within the scope of this work.
Unlike in the discrete regime, for high values of κ, the propagation takes place in a continuum fashion. By this we mean that successive nodes become excited before the excitation at a particular node reaches its maximum. As a consequence, there is always more than one node with its state moving on the fast manifold. This can be seen in the exemplary phase portraits Figs. 4 and 5.
We do not find distinct values of κ at which the overall behavior of f sust changes drastically as in the discrete regime. We observe the excitation of secondary waves in the entire continuum regime. Again, as the generation of secondary waves occurs only pairwise, this mechanism can not lead to the decay of all activity directly by pairwise annihilation of counterpropagating waves.
This can be seen in the exemplary time series Fig. 5(b), where, at the time of the snapshots (grey horizontal line in the space-time plots), the annihilation of two wave pairs is taking place as well as propagation failure. The corresponding positions in the phase portraits of Fig. 5(b) are annotated by arrows. The propagation failure of a wave is different from the annihilation of a wave pair in that many nodes cross the middle part of the u-nullcline or rather the threshold trajectory (almost) simultaneously.
We find as the main mechanism for propagation failure again the distribution of too many ends of the additional links in a small region, thus coupling back nodes that are in the excited state too strongly to nodes having a low-activator concentration. This effect can be seen very nicely in the (u, v) diagrams in Figs. 4(b) and 5(b), where all nodes that constitute the original wave are pulled over the middle-part of the

u-nullcline.
We also observe another notable effect: For nearest neighbor numbers higher than k = 1, there is a fixed coupling strength (depending on k but not on N ) at which f sust starts to decrease already at lower numbers n of additional links. This can be seen very well in Fig. 9(c) at κ ≈ 1. So far we have no explanation for the mechanism behind this phenomenon.

High κ -continuum limit
In order to include the effect of the long-range links into the continuum limit description of Sec. 2.3, we split the adjacency matrix A ij in Eq. (1) into two parts. A ij ≡ R ij + S ij , with R ij being all links of the original ring network and S ij being the additional randomly added links. Thus the dynamics readṡ where in the last equality a rescaling κ →κ = κq(k) (see Sec. 2.3) has been used, and the tilde will be dropped in the following. The ring part of the coupling can be treated in the same way as in Sec. 2.3. For the small-world part of the coupling, we assume a large number of additional links and distribute the entries in S ij equally over all entries of the entire matrix S ij which leaves S ij = 2n N 2 a constant. For easier readability, we consider only the small-world part of the coupling term s i : In performing the transition to the continuum description, we replace the sum N j=1 by the integral L 0 dy and introduce the mean valueū The continuum limit including the additional long-range links reads x ∈ [0, L] and (u, v)(t, 0) = (u, v)(t, L) , with L = N √ q(k)κ and σ = 2n q(k)L . With this mean-field approximation, the four coupling parameters (N, k, κ and n) reduce to two parameters (L and σ). This kind of global feedback coupling has also been studied for the Rinzel-Keller model in [47].

Approximate boundary of wave propagation
If σ = 0, Eq. (9) is the same as Eq. (8). Employing the same methods as in Sec. 2.3, we examine the change of the dispersion relation c(L) if σ is increased, see Fig. 10.
If σ = 0, the destabilization happens by a torus bifurcation at L cr (σ = 0). If σ is increased, L cr increases as well, i.e., the parameter range of L for stable propagation becomes smaller ( Fig. 10(a)). At a certain value of σ, the mechanism of destabilization changes, when the torus bifurcation coincides with a saddle-node bifurcation (limit point) at σ ≈ 0.246. After that, the destabilization happens by a saddle-node bifurcation. However, the σ value of the locus of the destabilizing saddle-node bifurcation goes into saturation at a value of σ max ≈ 0.247. This means that L cr goes to infinity when σ approaches σ max from below, so that above σ max , no stable propagation is possible. We display these loci of destabilization as a curve σ 0 (L) in Fig. 10(a). This curve can be transformed to a curve n 0 (N, k, κ), yielding an approximation for the boundary in n above which no realizations of a small-world network will support stable traveling waves. The curve L cr (σ) in (a) is connected with the instability points of the dispersion relations shown in Fig. 10(b) as indicated by the vertical dotted lines.
In Fig. 9, f sust is shown in the entire range of κ and n for N = 1000 and k = 1, 2, 3 as well as for N = 200, k = 2. Also shown in these plots is the analytical curve n 0 (κ) for the respective values of N and k which would give the expected maximum long-range link number. The transition in n to quenched wave activity happens at lower values of n. This is expected, as a significant contribution by the coupling term arising through the long-range links can only occur if the difference in activator concentration at both ends of the link is large. This is only the case if the node at one end of the shortcut is in the excited state (wave peak). Thus the critical link density is only important in part of the network. Of course, in a random network this is more likely to occur in an (arbitrary) part than in the entire network. Also note that the approximation becomes worse for higher k. Moreover, there is always an optimum coupling strength κ where n can be highest without disturbing the propagation of the wave. This optimum κ is a result of the transition between the discrete and the continuum regime.

Conclusion
We have studied the propagation of a solitary pulse (or wave) on a ring network and the influence of small-world perturbations of the topology upon the propagation. Already on the unperturbed ring topology, there are two regimes. One regime corresponds to high coupling strength, in which the behavior of the system resembles that of a continuous reaction-diffusion system with mean-field coupling. The other regime is associated with low coupling strength, in which the discrete nature of the network is important and the behavior differs from that of a reaction-diffusion system.
In each regime, a too large number of long-range links leads to failure of wave propagation. However, the mechanisms which lead to the suppression of the traveling  wave differ in the two regimes.
We have identified three different subregimes of coupling strength in the discrete, weak-coupling regime, which are sharply separated from each other. In the first one (lowest coupling strength), one additional link, regardless of the distance it spans, is enough to prevent propagation. In the second one, one additional link can be sufficient to prevent propagation if the distance it spans is not too large. For coupling strengths above the second subregime, secondary wave pairs can be created through the longrange links. For the latter coupling strengths, the mechanism for the quenching of a traveling wave is similar to that in the continuum regime.
In the strong-coupling regime, the main mechanism appears to be a too large number of additional links in the excited part of the wave (high activator concentration). Those links collectively "pull" the excited part back over the threshold trajectory of the system and thus lead to propagation failure. We have successfully approximated this behavior in the continuum limit by including a mean-field coupling term in the equations of the continuous reaction-diffusion system.