Impact of viscoelastic coupling on the synchronization of symmetric and asymmetric self-sustained oscillators

We investigate a system of two viscoelastically coupled, modified Van der Pol oscillators to compare their synchronization properties due to elastic and viscoelastic coupling. We show that viscoelastic coupling leads to in-phase synchronization while elastic coupling favours anti-phase synchronization. To study the impact of symmetry and nonlinearity, the restoring forces in the Van der Pol oscillators are extended to include nonlinear and asymmetric components. If the asymmetry, or rather the nonlinearity, of the restoring forces exceeds a certain threshold, only in-phase synchronized motion is found to be stable. Another important finding is that chaotic solutions can only be found if the restoring forces are asymmetric and the coupling incorporates viscosity.


Introduction
In nature and engineering a variety of systems exist whose macroscopic behaviour critically depends on the synchronization of the oscillating components of which they are composed. To understand such systems, much effort has been made to analyse systems composed of coupled self-sustained oscillators [1][2][3][4][5][6][7]. Most of these studies assume a purely elastic coupling neglecting coupling forces with viscous properties. Other studies, although incorporating some viscous coupling into their mathematical models and experiments, deal with very complex or specialized systems [8][9][10] and do not tackle the general question of how viscoelastic coupling influences synchronization compared with purely elastic coupling.
An example of a biological system where viscoelastic coupling is crucial is collagen gel embedded and periodically contracting premature cardiomyocytes, which may synchronize and assemble into an engineered heart muscle [11,12]. During this process the cardiomyocytes synchronize while their surrounding medium changes its viscous properties. It is hypothesized that this synchronization is mainly driven by the mechanical interaction of the cardiomyocytes via the viscoelastic medium in which they are embedded. This system is an example where the viscoelastic coupling of self-sustained oscillating components plays an important role in their synchronization behaviour.
To study the impact of viscoelastic coupling we coupled two Van der Pol oscillators [13] with the Maxwell model [14]. Furthermore, we modified the classical Van der Pol equations in such a way that they are able to produce asymmetric oscillations in order to analyse how this affects their synchronization properties. In the first section of this article we introduce the mathematical model we analysed, explaining its parameters and naming the numerical tools we used. After that, in section 3, we present our results dealing with viscoelastic coupling and the impact of nonlinear or rather asymmetric oscillations. Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Model and methods
In this section we introduce and explain the system of differential equations we are going to analyse, and the numerical methods employed.

Model equations
The model equations basically consist of two parts: self-sustained oscillators and their coupling. As already mentioned, the oscillators we use are modified Van der Pol oscillators. The dynamics of a single Van der Pol oscillator is given by the following partially dimensionless 6 equation: Here ω is the natural frequency of the oscillator and μ describes the time scale on which the nonlinear damping acts. Integrating this equation in time yields a symmetric oscillation (symmetric in the x-ẋ plane). To generate non-symmetric oscillations we perturb the linear restoring force, given by a harmonic potential, by an exponential function: where ε is the perturbation parameter which is small in magnitude with respect to all other parameters and | | c describes the length scale of the asymmetry of the former harmonic potential. The sign of c describes which side of the harmonic potential is influenced more strongly. Figure 1 shows how the introduced perturbation influences the harmonic potential and what the resulting oscillation looks like.
Furthermore, we add a cubic term, as in the Van der Pol-Duffing oscillator, to the restoring force described in equation (2). In this way we are able to compare the effects of an asymmetric nonlinearity with a symmetric nonlinearity and thus reveal the pure influence of an asymmetric restoring force: Here  b 0 is the parameter which controls the strength of the symmetric nonlinear part of the restoring force. The impact of hydrodynamic interactions between Van der Pol-Duffing oscillators has been studied in other work [15]. There the parameter controlling the nonlinearity b was called α. One important finding was that α could be used to control which type of synchronization (inor anti-phase) was found.
The full system of differential equations consists of two oscillators given by equation (3) which are viscoelastically coupled with the Maxwell model. This model consists of a harmonic spring in series with a linear damper. The spring constant of the spring is w c 2 and the time scale on which the damping acts is g -1 . We call w c elastic coupling parameter and γ fluidity (inverse viscosity) or the viscous coupling parameter. Putting everything together yields the following system of differential equations: A sketch of the above system is shown in figure 2. For convenience, we parameterize the parameter space of the natural frequencies w w ( ) , 1 2 T by a centre of mass like parameter system w w (¯) , r T : Here w is the mean of the natural frequencies and w r is their relative deviation with respect to the mean. Furthermore, we are still free to choose the time scale in our system of differential equations and use this to fix w = 1. Table 1 summarizes all system parameters with their default values. It is interesting to note that the last line of equation (4), representing the viscous creep, introduces a distributed time delay coupling with an exponential integral kernel between the two oscillators under consideration. This can be shown by solving the linear but inhomogeneous differential equation which determines the dynamics of d by a Greenʼs function approach. Doing this yields: 0 is the initial condition which decays to zero over a time scale given by g -1 . Furthermore, g -1 can be considered as the amount of temporal memory possessed by the distributed time delay coupling. Thus, small viscous coupling parameters γ correspond to a long temporal memory and large γ to a short one.
Other studies have analysed the influence of distributed time delay coupling on amplitude death [16,17], amplitude and phase dynamics [18], and the synchronization of networks [19,20] using identical and nonidentical Stuart-Landau oscillators. Most of those studies discussed different kinds of time delay distribution and their influence on the mentioned topics. In this work the type of distribution will not change but the influence of the memory length on the synchronization behaviour of in-and anti-phasic synchronization states of two non-identical Van der Pol like oscillators will be analysed. Although we emphasize and focus on the meaning of γ in the context of viscoelastic coupling, our results might be interesting in the field of coupled oscillators with distributed time delay as well and we encourage the reader to keep this second possible interpretation of d and γ in mind while reading this article.

Numerical methods
For bifurcation analysis we used the software package AUTO [21,22] where one has to specify a number of numerical parameters. To ensure the reproducibility of our results, we will provide the values of the most critical parameters 7 which are: = -EPSL 10 7 , = -EPSU 10 7 and = -EPSS 10 6 . Apart from AUTO we used the dopri5 routine from the PYTHON [23] library SCIPY [24] for numerical integration. Furthermore, we calculated Lyapunov exponents [25] to distinguish quasi-periodic solutions from chaotic ones.

Results and discussion
As in systems with non-viscous but purely elastically coupled oscillators, our model equation (4) exhibits synchronization of periodic oscillations in terms of phase locking. For analysing the fundamental differences between viscous and elastic coupling, we restrict ourselves to periodic solutions with 1:1 phase locking. In this case, only two types of synchronization states remain: in-phase (1:1 phase locking) and anti-phase (−1:1 phase locking) synchronization. Since the two oscillators in our model have, in general, distinct natural frequencies, it is not possible to obtain perfect 1:1 or −1:1 phase locking in the sense of a constant phase lag of j D = 0 or j p D = , respectively.
Therefore, we define in-phase synchronization as  j p D < 0 2 or p j p < D <   ratio does not change much with respect to the viscous coupling parameter γ and it depends linearly on the elastic coupling parameter w c .

Synchronization with a symmetric linear restoring force
Each synchronization state corresponds to an asymptotically stable, periodic orbit in phase space. In the following we first consider the case of a symmetric and linear (c = 0, b = 0) restoring force and analyse the stability of the two periodic orbits, in-phase (1:1) and anti-phase (−1:1), with respect to the elastic and viscous coupling parameters w g and c , respectively. To do this we used AUTO to continue the saddle node bifurcations which are responsible for the change in stability of the orbits. The resulting phase diagram is shown in figure 5. The phase diagram shows a blue (dashed) and a purple (solid) saddle node bifurcation curve corresponding to the anti-phase and in-phase orbit, respectively. Above those curves the corresponding orbit exists. On the curves the stable orbits collide with their unstable counterpart and vanish due to a saddle node bifurcation. Below both curves only quasi-periodic solutions exist without any synchronization. The anti-phase and in-phase orbit are both located on the same two-dimensional torus. Quasiperiodic solutions are filling this torus.
Increasing the fluidity γ yields three distinct effects, which can be seen in the different parts (a), (b) and (c) of figure 5. If the two oscillators are coupled purely elastically (g = 0), the anti-phase orbit emerges before the in-phase orbit does when increasing the elastic coupling strength w c . As soon as the fluidity γ reaches a certain value, one observes a behaviour which is the other way round (a). At higher fluidity (b), a strictly increasing elastic coupling strength is necessary for the anti-phase orbit to emerge while the elastic coupling strength needed for the in-phase orbit to emerge stays nearly constant. Furthermore, there is a critical value of γ at which the system will never reach anti-phase synchronization regardless of the value of w c . If the fluidity γ is increased even further (c), a strictly increasing elastic coupling strength w c is needed for the in-phase orbit to emerge. This is because a fluidity approaching infinity would describe a medium which does not provide any diffusive coupling at all (e.g. liquid water).
These results indicate that systems composed of periodically oscillating components which are coupled viscoelastically may have a tendency to synchronize in-phase at low coupling strength. An example would be the system mentioned in the introduction, where premature cardiomyocytes synchronize in-phase while their surrounding medium gets stiffer and less fluid over time. This is in accordance with our results, since most of the paths (from small w c and large γ to large w c and small γ) the system can follow in the phase diagram (figure 5) lead to parameter values where the system exhibits in-phase synchronization first, while stable anti-phase oscillations may occur only for higher values of the elastic coupling parameter w c .

Synchronization with an asymmetric restoring force
We have analysed how the coupling parameters influence the synchronization of the two viscoelastically coupled Van der Pol oscillators (c = 0, b = 0). Now we want to find out what role the (a)symmetry (  | | c 0, b = 0) of their oscillations plays. We modified the dimensionless Van der Pol equation to produce asymmetric oscillations and introduced the parameter c which is the length scale of asymmetry in the potential generating the restoring force (see equation (2) and figure 1). To analyse how the resulting asymmetry of the oscillations influences the synchronization behaviour of the coupled oscillators we calculated a bifurcation diagram with respect to c, which is shown in figure 6 9 . For the symmetric case (c = 0) we see the two asymptotically stable, coexisting,  w c and viscous γ coupling parameters) showing the continuation lines of the saddle node bifurcations which give rise to the different synchronization states. Above each line, the corresponding synchronized asymptotically stable periodic orbits exist. 9 The L 2 -norm was chosen to emphasize the reflection symmetry at c = 0 and it is defined as follows: let T be the period of an orbit, N the number of dimensions of the system and u k the dynamic variables. Then periodic orbits corresponding to the in-and anti-phase synchronization states and their unstable counterparts. As the asymmetry increases ( > | | c 0), all orbits remain until a critical value is reached where they vanish due to a saddle node bifurcation (the anti-phase orbit vanishes before the in-phase orbit does). After that point, there is a region in which no synchronization of periodic oscillations occurs. As can be seen in figure 7, the system has a chaotic solution in this region. If the oscillations become even more asymmetric a saddle node bifurcation occurs which gives rise to a single asymptotically stable in-phase orbit which is stable for a wide range of | | c values. The oscillations of this synchronization state are shown in figure 8. Now we want to analyse how the existence of chaotic solutions depends on the viscous effects of the coupling. Figure 9 shows the two largest Lyapunov exponents l 1 and l 2 of the system in γ-c parameter space. One can see that the system exhibits chaos (l > 0 1 ) only for a specific c-range, which depends on the fluidity γ. We checked this result by continuing the two saddle node bifurcations marked in figure 6 and plotting their bifurcation curve on top of the Lyapunov exponent map. By increasing γ, the bifurcations come closer together until they annihilate each other, leaving the stable in-phase solution without any chaotic regime. It is important to note that there is no chaos to be found at g = 0 (l = 0 1 ). In this case, the second largest Lyapunov exponent is equal to zero as well. This can be explained by looking at the first derivative of the dynamic variable d in equation (4). Since g =  ȯ d 0 0a perturbation in d will neither grow nor decay. This is also true for a perturbation in the direction of the limit cycle corresponding to one of the two synchronization states. Because there are two dimensions in phase space in which a perturbation will be constant on average, there have to be two Lyapunov exponents equal to zero. Furthermore, we have calculated the Lyapunov exponents in w c -c parameter space with g = 0. There was no chaos to be found. So chaotic oscillations occur for the coupled system only with viscous coupling (g > 0) and asymmetric restoring forces ( > | | c 0).

Synchronization with a symmetric nonlinear restoring force
We have studied the synchronization dynamics of a system with asymmetric restoring force. The asymmetry was introduced by perturbing the former linear resorting force by an nonlinear exponential term, equation (2). Now we want to analyse if the dynamic behaviour found while studying the asymmetric system ¹ | | c 0 is due to the asymmetry or to the nonlinearity of the restoring force. For this purpose we study the system equation (4) with a symmetric but nonlinear (c = 0,  b 0) restoring force: two coupled Van der Pol-Duffing oscillators. To analyse how the symmetric but nonlinear restoring force influences the dynamics of the coupled oscillators we calculated a bifurcation diagram based on a Poincaré map and the three largest Lyapunov exponents, as shown in figure 10. In the case of a linear restoring force (b = 0) we see the two coexisting periodic orbits corresponding to the in-(purple) and anti-phase (blue) synchronization states. As the strength of the nonlinearity in the restoring force increases ( > b 0), the orbits remain until a critical value » b 0.2 is reached. At this point the anti-phase solution vanishes. The in-phase orbit remains until another critical point » b 1.6 is reached, at which the in-phase solution vanishes too. After this point no 1:1 synchronous state is found anymore. This was checked by calculating the Poincaré map and the Lyapunov exponents up to b = 20. Furthermore, by varying the initial condition at different values of b, no other coexisting attractors were found. Figure 11 shows the largest two Lyapunov exponents of the system in the γ-b-parameter space. It can be seen that the described behaviour does not depend qualitatively but quantitatively on the viscous coupling parameter γ. The second critical value of b, after which no synchronized solution exists (blue area), increases with γ. Thus a high fluidity seems to have a stabilizing effect. If γ becomes too large, this effect diminishes, since very high fluidity effectively decreases the coupling strength between the two oscillators.
In contrast to the results obtained while analysing the asymmetric nonlinear restoring force (section 3.2) the regime where neither an in-phase nor an anti-phase orbit exists is not chaotic, as the Lyapunov exponents ( figure 10(b)) indicate. Furthermore, there is no new in-phase orbit to be found as the parameter controlling the strength of the nonlinear part in the restoring force increases. Another difference to note is that the anti-phase  orbit vanishes long before the in-phase orbit does. In summary, it was found that the in-phase orbit is favoured by an increasing nonlinearity | | c or b in the restoring force; independent of symmetry. However, in system equation (4) only with a nonlinear asymmetric restoring force was chaotic dynamics observed.

Conclusion
The aim of this work was to analyse the effects of viscoelastic coupling on the synchronization behaviour of coupled self-sustained oscillators in comparison with purely elastic coupling. To perform this study we used modified Van der Pol oscillators that are able to exhibit asymmetric oscillations due to an asymmetry in the potential generating the linear restoring forces. With respect to the synchronization of two coupled oscillators, we restricted ourselves to 1:1 and −1:1 phase locking.
We found that in general the in-phase synchronization state is favoured by both more fluid coupling and oscillations generated with a nonlinear restoring force. This means in particular that if the fluidity exceeds a critical value or the nonlinearity in the restoring force becomes too large anti-phase synchronization ceases to exist, independent of the elastic coupling strength.
In the system with a linear restoring force it was found that, if the fluidity of the coupling is high enough, the in-phase synchronization state emerges before the anti-phase synchronization state when increasing the elastic coupling strength.
Another important finding was that viscous effects in the coupling and asymmetry of the restoring forces are necessary conditions for the emergence of chaos in the analysed system (in the investigated part of the parameter space). This means that, in comparison with purely elastic coupling, viscoelastic coupling yields not only quantitative but important qualitative differences.  ). It is important to note that for vanishing viscous coupling parameter γ an additional Lyapunov exponent is equal to zero since = d 0.