## Abstract

This paper studies the effects of distributed-delay coupling on the dynamics in a system of non-identical coupled Stuart–Landau oscillators. For uniform and gamma delay distribution kernels, the conditions for amplitude death are obtained in terms of average frequency, frequency detuning and the parameters of the coupling, including coupling strength and phase, as well as the mean time delay and the width of the delay distribution. To gain further insights into the dynamics inside amplitude death regions, the eigenvalues of the corresponding characteristic equations are computed numerically. Oscillatory dynamics of the system is also investigated, using amplitude and phase representation. Various branches of phase-locked solutions are identified, and their stability is analysed for different types of delay distributions.

## 1. Introduction

The dynamics of many complex physical, biological and engineering systems can be effectively modelled mathematically using ensembles of coupled oscillators [1,2]. A significant advantage of such an approach over other methodologies lies in the possibility to identify and study regimes with particular types of behaviour, such as emergence and stability of cooperation and (de)synchronization, different types of spatial patterns, travelling waves, etc. Such analyses often provide very important insights for practical applications. Our understanding and treatment of various brain pathologies and deficiencies, such as Parkinson's, Alzheimer's and epilepsy, relies heavily on the analysis of synchronization of oscillating neural populations [3–5]. Recent work on laser communication networks has extensively used coupled oscillator models to study the dynamics of in-phase or anti-phase and complete chaos synchronization [6–10]. Chimera states, where a network of oscillators splits into coexisting domains of coherent, phase-locked and incoherent, desynchronized behaviour, have also been observed [11–15]. These studies have shown that coupling between different elements can play a dual role in the dynamics: it can both lead to suppression of oscillations and also facilitate a certain degree of synchronization between different elements.

When analysing the dynamics of coupled systems, it is important to take into account the fact that, in many cases, the coupling between different elements is not instantaneous. Time delays associated with such coupling can have a major impact on the overall dynamics of the system. In practical examples, these time delays are often associated with delays in propagation or processing of signals, response time of mechanical actuators, translation and transcription time in genetic oscillators or driver reaction time in traffic dynamics [16–21]. In all these examples, explicit inclusion of time delays in the model has provided a more realistic and accurate representation of the system under consideration. From a mathematical perspective, differential equations with time delays are infinite-dimensional dynamical systems, which makes their analysis, as well as numerical simulations, quite involved [22,23].

In the case when the coupling between oscillators is sufficiently weak, it is possible to reduce the full system to a phase-only model [1]. A significant amount of research has been carried out over the years on the analysis of phase oscillators with different types of time-delayed coupling [2,24–32]. Such systems of phase oscillators demonstrate a rich variety of dynamical regimes, including chaos, synchronization, splay and chimera states, characterized by the coexistence of coherent and incoherent states. In many cases, coupling in the phase models leads to phase entrainment and the emergence of the so-called phase-locked solutions, where all oscillators start to oscillate with the same common frequency and have a constant shift between their phases. The existence and stability of such phase-locked solutions have been studied in a number of systems of oscillators with different types of local and non-local coupling [14,31].

For sufficiently strong coupling between oscillators, one can observe another interesting and important aspect of dynamics of coupled oscillator systems: the ability of external coupling to suppress otherwise stable periodic oscillations. This was first discovered by Bar-Eli in the context of chemical reactions with instantaneous coupling [33], and it has been shown for a system of two coupled Stuart–Landau oscillators that the ‘Bar-Eli effect’ can only occur provided the coupling strength and the frequency detuning of the two oscillators are both sufficiently large [34–36]. Later, it was discovered that, when the coupling is time-delayed, it is possible to achieve suppression of oscillations and stabilization of an unstable fixed point even for identical oscillators [37,38]. This phenomenon, named amplitude death [37], oscillator death, or ‘death by delay’ [39], has been demonstrated experimentally in nonlinear electronic circuits [40], in the dynamics of the slime mould *Physarum polycephalum* [41] and in thermo-optical oscillators linearly coupled by heat transfer [42]. Amplitude death has been subsequently studied for a number of different systems and couplings [31,43–52].

Despite successes in the analysis of systems of coupled oscillators with time-delayed coupling, one of the limitations of the majority of this research has been the restriction on the type of time-delayed coupling, which is usually taken to be in the form of one or several constant time delays. At the same time, in many realistic systems, the time delays themselves are not constant [53,54] and may either vary depending on the values of system variables (state-dependent delays) or just not be explicitly known. In order to account for such situations mathematically, one can use the formalism of distributed time delays, where the time delay is represented through an integral kernel describing a particular delay distribution [55–57]. Distributed time delay has been successfully used to describe situations when only an approximate value of time delay is known in engineering experiments [58,59], for modelling distributions of waiting times in epidemiological models [60], maturation periods in population and ecological models [61,62], as well as in models of traffic dynamics [63], neural systems [64], predator–prey and food webs [65].

In this paper, we investigate amplitude death and phase dynamics in a system of coupled Stuart–Landau oscillators with distributed-delay coupling. This system is a prototype for dynamics near a supercritical Hopf bifurcation, and, in this capacity, it captures the essential features of many realistic systems in such a regime. The corresponding mathematical model can be written in the form
1.1where , *ω*_{1,2} are the frequencies of the two oscillators, and are the strength and the phase of coupling, respectively, and *g*(⋅) is the kernel of the delay distribution. The kernel *g* is taken to be positive-definite and normalized to unity:
The case *g*(*u*)=*δ*(*u*) corresponds to instantaneous coupling (*z*_{2}−*z*_{1}), describing a situation when oscillators interact without any delay. For *g*(*u*)=*δ*(*u*−*τ*), the coupling takes the form of a discrete time delay: *z*_{2}(*t*−*τ*)−*z*_{1}(*t*). Atay [66] has analysed this system for the case of zero coupling phase and a uniform delay distribution kernel analytically for the case of identical oscillators and numerically for non-identical oscillators, and identified regimes of amplitude death in terms of coupling strength and mean time delay. Kyrychko *et al.* [67] have analysed amplitude death in model (1.1) with *ω*_{1}=*ω*_{2}=*ω*_{0} for uniform and gamma distributed-delay kernels and non-zero phase.

The outline of this paper is as follows. In the next section, we analyse the stability of the trivial equilibrium of system (1.1) and find regions of amplitude death depending on the parameters of the coupling and delay distribution. The eigenvalues of the corresponding characteristic equations determining the stability of the steady state are computed numerically to gain a better understanding of system dynamics inside amplitude death regimes. Section 3 contains an analysis of phase-locked solutions of system (1.1) for uniform and gamma distributions in the case of constant equal amplitudes of oscillations (when the system reduces to a coupled system of Kuramoto oscillators with distributed-delay coupling), as well as in the general case of arbitrary amplitudes. In each case, we identify branches of phase-locked solutions and numerically compute their stability. The paper concludes with a discussion of the results as well as possible further developments of this work.

## 2. Amplitude death

To study the possibility of amplitude death in the system (1.1), we linearize this system near the trivial steady state *z*_{1,2}=0. The corresponding characteristic equation is given by
2.1where *λ* is an eigenvalue of the Jacobian, and
2.2is the Laplace transform of the function *g*(*u*). Atay [66] has investigated the case *θ*=0 analytically for identical oscillators *ω*_{1}=*ω*_{2}=*ω*_{0}, and numerically for *ω*_{1}≠*ω*_{2} and a uniform delay distribution kernel. More recently, Kyrychko *et al.* [67] have studied analytically and numerically the case of *ω*_{1}=*ω*_{2}=*ω*_{0} and a non-vanishing coupling phase *θ*≠0 with uniform and gamma delay distributions.

### (a) Uniformly distributed kernel

To make further analytical progress, it is instructive to specify a particular choice of the delay kernel, which would also be relevant for applications. As a first example, we consider a uniformly distributed kernel 2.3This distribution has the mean time delay and the variance 2.4The uniformly distributed-delay kernel (2.3) has been successfully used in a number of different contexts, including ecological models [68,69] where it describes a lag between environmental changes and the reproductive response of species, models of traffic dynamics with delayed driver response [63], stem cell dynamics [70], time-delayed feedback control [53] and genetic regulation [71].

In the case of uniformly distributed kernel (2.3), it is quite easy to compute the Laplace transform of the distribution *g*(*u*) as
and this transforms the characteristic equation (2.1) into
2.5Because the roots of the characteristic equation (2.5) are complex-valued, the stability of the trivial steady state can change only if some of these eigenvalues cross the imaginary axis. To this end, we can look for characteristic roots in the form *λ*=i*ω*. Substituting this into the characteristic equation (2.5) and separating real and imaginary parts gives the following system of equations for (*K*,*τ*) as parametrized by the Hopf frequency *ω*:
2.6where
and we have introduced the frequency mismatch (detuning) *Δ* and the mean frequency as

Squaring and adding the two equations in (2.6) gives a single quartic equation for the coupling strength *K*:
2.7In a similar way, dividing the two equations in (2.6) gives the equation for the time delay *τ* at the Hopf bifurcation as
2.8

When the coupling phase vanishes (*θ*=0), the above equations for (*K*,*τ*) simplify to
2.9For *Δ*=0, we have two identical oscillators with the same frequency *ω*_{1}=*ω*_{2}=*ω*_{0}. Substituting this into equation (2.7) gives the equation for the Hopf frequency in the form
2.10and the equation for the critical time delay *τ* at the Hopf bifurcation reduces to

Using the trigonometric identity , we find
2.11Equations (2.10) and (2.11) have been recently studied by Kyrychko *et al.* [67], where the effects of the width of delay distribution *ρ*, as well as coupling strength *K* and the coupling phase *θ*, on the amplitude death were investigated. In the case of vanishing coupling phase (*θ*=0), these equations reduce even further to the system studied by Atay [66].

To illustrate the effects of varying the coupling strength *K* and the time delay *τ* on the (in)stability of the trivial steady state, we now compute the stability boundaries (2.7) and (2.8) as parametrized by the Hopf frequency *ω*. Besides the stability boundaries themselves, which enclose the amplitude death regions, we also compute the maximum real part of the eigenvalues using the traceDDE package in Matlab. In order to compute these eigenvalues, we introduce real variables *z*_{1r,i} and *z*_{2r,i}, where *z*_{1}=*z*_{1r}+*iz*_{1i} and *z*_{2}=*z*_{2r}+*iz*_{2i}, and rewrite the linearized system with the distributed kernel (2.3) as
2.12where
and **0**_{2} denotes a 2×2 zero matrix. When *ρ*=0, the last term in the system (2.12) turns into *KM***z**(*t*−*τ*), which describes the system with a single discrete time delay *τ*. System (2.12) is in a form that is amenable to the algorithms described in Breda *et al.* [72] and implemented in traceDDE.

First of all, we consider the case *ρ*=0 and *θ*=0, corresponding to a discrete time delay and vanishing coupling phase, as shown in figure 1. This figure shows that, as the frequency detuning increases, this leads to the emergence of new islands of amplitude death, and for sufficiently high detuning *Δ*, these islands merge into a single continuous region in the plane of coupling strength *K* and average time delay *τ*, where amplitude death can occur for an arbitrary value of *τ*, provided *K* lies in the appropriate range. One can note that, unlike the case of identical oscillators considered in Kyrychko *et al.* [67], in this situation, the values of *K* needed to achieve stabilization of the trivial steady state for any *τ* are in the lower part of the overall range of *K* values. Increase in the width of the uniform delay distribution *ρ* leads to a corresponding increase in the region of amplitude death, as illustrated in figure 2. When we consider the impact of coupling phase on amplitude death, as shown in figure 3, it becomes clear that the largest range of admissible *K* values is attained for *θ*=0, but overall the range of coupling phases, for which amplitude death is possible, increases with the frequency detuning *Δ*. It is worth noting, however, that, unlike the situation considered in Kyrychko *et al.* [67] where increase in the width of the distribution *ρ* led to the asymmetry of amplitude death regions with regard to the coupling phase, as *Δ* is increased, the region of amplitude death remains symmetric in *θ* (for small *ρ*).

### (b) Gamma distribution kernel

The second example we consider is that of a gamma distribution, which can be written as
2.13with *α*,*p*≥0, and *Γ*(*p*) being the Euler gamma function defined by *Γ*(0)=1 and *Γ*(*p*+1)=*pΓ*(*p*).

For integer powers *p*, this can be equivalently written as
2.14For *p*=1, this is simply an exponential distribution (also called a *weak delay kernel*) with the maximum contribution to the coupling coming from the present values of variables *z*_{1} and *z*_{2}. For *p*>1 (known as *strong delay kernel* in the case *p*=2), the biggest influence on the coupling at any moment of time *t* is from the values of *z*_{1,2} at *t*−(*p*−1)/*α*. The delay distribution (2.14) has the mean time delay
2.15and the variance
A gamma distributed-delay kernel (2.14) was originally analysed in models of population dynamics [73–76], and has subsequently been used to study machine tool vibrations [77], intracellular dynamics of HIV infection [78], traffic dynamics with delayed driver response [79] and control of objects over wireless communication networks [80].

When studying the stability of the trivial steady state of the system (1.1) with the delay distribution kernel (2.14), one could use the same strategy as the one described for a uniform distribution. The Laplace transform of the distribution kernel in this case is given by
Substituting this into the characteristic equation (2.1) yields a polynomial equation for *λ*:
2.16In appendix A, we show how the same characteristic equation can be derived using a linear chain trick without resorting to the Laplace transform.

Figure 4 illustrates the regions of amplitude death for a weak delay distribution kernel (2.14) with *p*=1, vanishing coupling phase and increasing frequency detuning *Δ*. Similar to the case of uniform delay distribution, as *Δ* increases, new islands of amplitude death appear, and eventually they merge into a single continuous region of amplitude death. Because, in this figure, the regions of amplitude death are plotted in terms of *α*, which according to (2.15) is the inverse average time delay *τ*, this implies that, for the case of a single connected region of amplitude death in the (*α*,*K*) plane, amplitude death can happen for an arbitrarily small value of the average time delay, provided that the frequency detuning is sufficiently large. In figure 5, we illustrate how amplitude death regions depend on the coupling phase. This figure suggests that, for the same average time delay, provided the coupling phase is sufficiently negative, it is possible to achieve amplitude death for an arbitrary value of the coupling strength *K* starting from some minimal value. The regions of amplitude death are strongly asymmetric in *θ*, and amplitude death is possible for higher positive values of *θ* for larger frequency detuning *Δ*.

When one considers a strong distribution kernel (2.14) with *p*=2 with vanishing coupling phase, there is a minimum value of the average time delay required to achieve amplitude death, as shown in figure 6. This figure also suggests that the actual value of the minimum average time delay increases with the increasing frequency detuning *Δ*. Figure 7 illustrates how the coupling phase affects amplitude death. One can note that, similar to the case of a weak distribution kernel, the regions of amplitude death are asymmetric in *θ*. However, there are two major differences from the case *p*=1, namely that the regions of amplitude death are closed in the (*K*,*θ*) parameter plane, and these regions shrink with increasing *Δ* rather than grow as was the case for *p*=1. The implication is that now it is not possible to find a coupling phase for which amplitude death could be achieved for an arbitrary value of the coupling strength *K*.

## 3. Phase-locked solutions

In the previous section, we analysed the situation when the distributed time delay in the coupling leads to the destruction of stable limit cycles and stabilization of the previously unstable fixed point. In order to understand the phase dynamics of the system, we introduce new real variables (*R*_{1}(*t*),*R*_{2}(*t*),*ϕ*_{1}(*t*),*ϕ*_{2}(*t*)) as follows:
Here *R*_{1,2}≥0 and *ϕ*_{1,2} denote the amplitudes and phases of the two oscillators, respectively. Substituting this representation into the system (1.1) yields the following system of equations for the amplitude and phase variables:
3.1A significant amount of work has been carried out on the analysis of this system for the case of instantaneous coupling *g*(*s*)=*δ*(*s*) or discrete time delay *g*(*s*)=*δ*(*s*−*τ*). In both of these cases, it has been shown that, besides amplitude death, the system can exhibit a number of other interesting solutions, such as phase-locked (also known as frequency-locked) and phase drift solutions [81–84]. In this section, we consider effects of distributed-delay coupling on the dynamics of such solutions, which have so far remained unexplored.

### (a) Constant amplitude dynamics

As a first step in the analysis of system (3.1), we assume that the amplitudes of both oscillators are equal to each other and constant: *R*_{1}(*t*)=*R*_{2}(*t*)=const for all times. In this case, neglecting amplitude dynamics, the system (3.1) reduces to a system of two *Kuramoto oscillators with distributed-delay coupling*:
3.2It is instructive to consider this system using the variables of the phase difference *ψ*=*ϕ*_{1}−*ϕ*_{2} and the mean phase *φ*=(*ϕ*_{1}+*ϕ*_{2})/2:
3.3The last equation suggests that, when *Δ*>2*K*, this system exhibits phase drift (i.e. unbounded growth of the difference between the phases of two oscillators) independently of the delay distribution kernel, and hence phase locking can occur only for *Δ*<2*K*. Following Schuster & Wagner [24], we look for solutions of the system (3.3) in the form
3.4where *Ω* is the ensemble frequency of oscillations, and *β* is the constant phase shift between the two oscillators. Substituting this into the system (3.3) yields the value of the phase shift as
3.5where we have introduced the notation
The new common frequency *Ω* satisfies the transcendental equation
3.6where the plus and minus sign correspond to the first and second value of *β* in (3.5), respectively. It follows from this equation that possible solutions for *Ω* can only lie in the admissible range
Stability of the phase-locked solution (3.4) is determined by the roots of the corresponding characteristic equation
3.7where
and the superscript L refers to this integral representing the Laplace transform of a modified kernel .

To make further analytical progress and derive results relevant for applications, one has to specify the form of the delay distribution kernel, which is taken to be the same as in the analysis of amplitude death in the previous section. For the uniform distribution (2.3), functions *F*_{c} and *F*_{s} can be computed as
and the case of a discrete time delay can be recovered by setting *ρ*=0. For the weak distribution kernel (2.14) with *p*=1, we have
and similarly, for the strong delay kernel (2.14) with *p*=2, we have

Figure 8 illustrates how the function *f*_{±}(*Ω*), whose roots determine the values of the collective frequency *Ω* of phase-locked branches, depends on frequency detuning * Δ* and the coupling strength

*K*. Computation of the stability of the branches shown in this figure indicates that, as the frequency detuning increases, branches of phase-locked solutions start to appear for higher values of the coupling strength.

A similar computation for the weak delay distribution kernel, shown in figure 9, suggests that in this case all branches of phase-locked solutions are unstable for any value of *Δ*. In the case of the strong distribution kernel, which is illustrated in figure 10, the branches of phase-locked solutions are stable for a very narrow range of *K* near the lowest tip of the branch for sufficiently small detuning (see insets). As the detuning increases, this causes the branches of phase-locked solutions to appear for higher values of the coupling strength in a manner similar to the case of uniform delay distribution kernel, and all these branches are unstable.

### (b) General phase-locked solutions

In the more general case, when *R*_{1} and *R*_{2} are not required to be equal to each other, one can still look for phase-locked solutions of the system (3.1) in the form
where are unknown constants, *Ω* is the new common frequency of oscillations and *β* is the phase shift between the two oscillators. These solutions can be found as the roots of the following system of equations:
3.8It is easy to check that the constant amplitude phase-locked solutions (3.4) are the only solutions of the full system (3.8) for the case of identical oscillators *ω*_{1}=*ω*_{2}, which implies *Δ*=0 and the phase shift *β* of either zero or *π* (describing in-phase or anti-phase solutions, respectively).

Although it is not possible to find solutions of the system (3.8) analytically, the phase-locked solutions can be computed numerically for each particular choice of delay kernel. Figure 11 shows the branches of phase-locked solutions (3.8) for the uniform delay distribution kernel. For *Δ*=0, we identify both the branch of equal amplitude phase-locked solution and also an unstable branch with unequal amplitudes. As the frequency detuning increases, stable branches of phase-locked solutions are observed for higher values of the coupling strength. For all values of detuning, as *K* tends to zero, the ensemble frequency of branches of phase-locked solutions approaches the values of *Ω*=*ω*_{1,2}, which should be expected from the fact that for *K*=0 the system (3.8) admits solutions (*R*_{1},*R*_{2},*Ω*)=(1,0,*ω*_{1}) and (*R*_{1},*R*_{2},*Ω*)=(0,1,*ω*_{2}).

Figure 12 illustrates the branches of phase-locked solutions for the weak distribution kernel. Similar to the case of uniform delay distribution kernel, for *Δ*=0, we find the branches of solutions with both equal and different amplitudes, and both of these branches are unstable. For higher values of frequency detuning, we have two unstable branches of phase-locked solutions located symmetrically around the average frequency . The distance between these branches in terms of their ensemble frequency is increasing with increasing *Δ*. For the strong distribution kernel, the situation is similar, as shown in figure 13. Once again, we observe two unstable branches of phase-locked solutions, whose frequencies are centred around , and the distance between the branches grows with *Δ*.

## 4. Discussion

In this paper, we have analysed amplitude death, as well as the existence and stability of phase-locked oscillatory solutions in a generic model of Stuart–Landau oscillators with delay-distributed coupling. Regions of amplitude death have been identified for uniform and gamma distributions depending on the intrinsic system parameters (average frequency and detuning), as well as coupling strength, phase and the distribution parameters. Computation of the boundaries of amplitude death regions, as well as the stability eigenvalues for the trivial steady state, suggests that the more disparate the oscillators are (i.e. the higher the frequency detuning), the larger are the regions of amplitude death in the case of uniform and weak delay distributions. However, for the strong delay distribution kernel, higher detuning corresponds to a smaller region of amplitude death. Special attention has been paid to the analysis of the effects of the coupling phase on amplitude death. While for uniform delay distribution, the maximum range of coupling strengths giving amplitude death is achieved for zero phase, in the case of gamma distribution, the non-zero coupling phase increases the range of such coupling strengths.

To understand the dynamics beyond amplitude death, we have analysed the appearance and stability of phase-locked solutions, characterized by both oscillators performing oscillations with the same common frequency and having a constant phase shift. When the amplitudes of both oscillators are constant and equal to each other, phase approximation results in a system of Kuramoto oscillators with distributed-delay coupling, in which case it is possible to find the range of admissible collective frequencies and the phase shift analytically. For a uniform delay distribution, the system of coupled Kuramoto oscillators exhibits both stable and unstable branches of phase-locked solutions, which appear for higher values of the coupling strength as the frequency detuning increases. For weak (exponential) delay distribution, all the branches of phase-locked solutions are unstable, and for the strong (gamma function) delay distribution kernel, they are stable for small detuning and sufficiently small coupling strength. We have also computed numerically branches of phase-locked solutions for uniform and gamma distributions for the full system. For a uniform delay distribution, there is coexistence of stable and unstable phase-locked branches, with branches being stable for higher values of *K* as the detuning increases. In the case of gamma distribution, there are two branches of phase-locked solutions, which are both unstable, and their frequency difference from the average frequency increases with increasing *Δ*.

So far, we have investigated phase-locked solutions in the system of Stuart–Landau oscillators with distributed-delay coupling. One possible extension of this work would be the analysis of a phase-flip transition [49,85], where in-phase and anti-phase phase-locked branches exchange stability. This phenomenon has been observed in systems with discrete time delays, and it has even been observed in the transient dynamics preceding amplitude death. Another possibility would be to consider whether coupled systems with delay-distributed coupling are able to exhibit other types of phase dynamics and synchronization. One such scenario, which is useful in laser applications, is the case when the sum of the phases of the two oscillators is constant [29,86]. Phase approximation in this case results in a delayed Adler equation, and it would be both theoretically and practically important to consider possible solutions of this model for different delay distributions. Other more complex locking scenarios have been found in quantum-dot lasers under optical injection [87,88].

## Funding statement

This work was partially supported by DFG in the framework of SFB 910: Control of self-organizing nonlinear systems: Theoretical methods and concepts of application.

## Acknowledgements

Y.N.K. and K.B.B. gratefully acknowledge the hospitality of the Institut für Theoretische Physik, TU Berlin, where part of this work was completed.

## Appendix A

A convenient way to derive the characteristic equation in the case of a gamma distributed-delay kernel is to use the *linear chain trick* [89], which allows one to replace the original equation by an equivalent system of (*p*+1) ordinary differential equations. To illustrate this, we consider a particular case of system (1.1) with a weak delay kernel given by (2.14) with *p*=1, which is equivalent to a low-pass filter [44]:
A1Introducing new variables
allows us to rewrite the system (1.1) as
A2where the distribution parameter *α* is related to the mean time delay as *α*=1/*τ*_{m}. The trivial equilibrium *z*_{1}=*z*_{2}=0 of the original system (1.1) corresponds to a steady state *z*_{1}=*z*_{2}=*Y* _{1}=*Y* _{2}=0 of the modified system (A2). The characteristic equation for the linearization of system (A2) near this trivial steady state is given by
A3which is the same as equation (2.16). For larger values of *p*, one would have to introduce a larger number of additional variables in a similar fashion (two additional variables for each increase of *p* by 1).

## Footnotes

One contribution of 15 to a Theme Issue ‘Dynamics, control and information in delay-coupled systems’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.