Synchronization Fronts in a Spatially Extended System of Hybrid Rayleigh-van der Pol Oscillators

Numerous biological systems exhibit transitions to synchronised oscillations via a population-density-dependant mechanism known as quorum sensing. Here we propose a model system, based on spatially distributed limit-cycle oscillators, that allows us to capture the dynamics of synchronization fronts by taking the continuum limit as the number of coupled oscillators tends to infinity. We explore analytically a family of wavefront-type solutions to the system in terms of model parameters and corroborate its validity via a numerical finite-difference method algorithm. Finally, we review our results in light of some experimental systems in synthetic biology based on a synchronised quorum of genetic clocks, coupled via a diffusive auto-inducer signalling molecule analogue to our Hooke-like coupling scheme within the proposed hybrid Rayleigh-van der Pol model.

Keywords— Coupled Oscillators, Genetic Clocks, Quorum Sensing, Synchronization Fronts, Limit-cycle, Travelling Waves, Autoinducer, Repressilator, Synthetic Biology.

Synchronization in a large population of interacting oscillatory elements is a remarkable phenomenon that has attracted attention in multiple research areas [1]. One of the first historical reports regarding this phenomenon was made by a Dutch physician named Engelbert Kaempfer after a journey to Siam during the second half of the 17th century. Kaempfer describes an event which happens if a swarm of fireflies occupies a tree: the insects “hide their lights all at once, and a moment after make it appear again with the utmost regularity and exactness" [2]. This mechanism of synchronous emission of light pulses by a population of fireflies has been studied later on by biologists [3] and mathematically through the study of simplified models of coupled oscillators [4].

A parallel biological phenomenon, but within the microscopic domain, has been investigated recently through the use of repressilator based synthetic gene regulatory networks of synchronized clocks coupled via quorum sensing [5][6][7]. These systems made of bacterial colonies are engineered to produce synchronised oscillations through a fluorescent protein reporter within a microfluidic device.

In both phenomena, composed by a group of interacting biological clocks and their environment, frequency synchronization is a logical approach since we are interested in states for which the temporal rates of change of the phases are equal, rather than the phases themselves.

One of the main models to study those kinds of systems that exhibit synchronization was first introduced by Winfree [8], and later refined by Kuramoto and others [9][10][11]. (A comprehensive review and perspectives of the field can be found in [12]). In the existing literature systems of fully coupled oscillators have mostly been considered due to its analytical tractability and simplicity . On the other hand, systems of spatially extended oscillators with local couplings have received less attention even though they are more realistic and frequently observed in nature [13][14][15][16][17].

Here we propose a system based on spatially distributed limit-cycle oscillators, which have been extensively used as a model of physical and biological clocks, based on a hybrid Rayleigh-van der Pol (RvdP) oscillator. The main idea is to capture the dynamics of synchronization fronts by taking the continuum limit as the number of coupled oscillators tends to infinity. This approach makes sense in many biological systems such as colonies of bacteria in which a network-like coupling is generated via a quorum sensing mechanism. Usually when the population density exceeds a certain threshold then the corresponding local levels of signal can induce an spatially extended synchronised response in gene expression throughout the entire population [18].

We firstly consider the second order homogeneous van der Pol differential equation, describing an oscillator with non-linear damping,

where $u\in\mathbb{R}$ is the dynamical variable, $\omega_{0}$ is the linear frequency and $\varepsilon>0$ is a non-linear scalar parameter of the system. An important aspect of this equation is the relation between the negative damping proportional to $-\dot{u}$, and the non-linear damping, $u^{2}\dot{u}$, which leads to the existence of self-sustained limit cycle oscillations.

A simple way to understand this equation is given by applying a transformation into a rotating reference frame [19], via

with a complex amplitude that changes slowly as

where $u$ and $v=\dot{u}$ denote the position and the conjugate moment respectively. Within this rotating frame one can neglect the terms corresponding to fast oscillations in (2) if $\varepsilon<<1$ holds. This allows us to write equation (2) as a generalised Stuart-Landau type equation

which describes the oscillator dynamics. Looking for stationary states, i.e. when $\dot{\alpha}(t)=0$, this equation admits a limit cycle defined by $|\alpha(t)|^{2}=1$.

In our case we introduce a modified version of equation (1), namely

which is of the form $\ddot{u}=f(u,\dot{u})$.

First of all we can observe that equation (5) is a combination of the previous van der Pol equation (2) and the classic Rayleigh equation defined by

In that sense the proposed equation (5) contains the cubic term proportional to $\dot{u}^{3}$ which is contained in the Rayleigh equation (6) and not in the van der Pol equation (5). Concurrently it contains the $u^{2}\dot{u}$ term which the Rayleigh equation (6) doesn’t have, but the van der Pol equation (5) does. This hybrid Rayleigh-van der Pol equation has been studied in an approximate way within the context of bipedal motion modelling [20].

Here we analyse equation (5) starting from the case in which $\delta=1/\omega_{0}^{2}$ and considering the spatial location of each oscillator explicitly. To begin with, we consider a discrete system assuming that each element $u_{i}$ is coupled to its nearest neighbours $u_{i-1}$ and $u_{i+1}$ via a coupling constant $\mu$ (see figure 1). According to this scheme we obtain

Adding and subtracting the term $2u_{i}$

we observe that the finite difference approximation of the second spatial derivative naturally appears.

In other words the term $u_{i+1}+u_{i-1}-2u_{i}$ is the discretisation of the partial derivative

so the ODE becomes a PDE in the folling way,

and finally for one spatial dimension we obtain

which can be rewritten as

Now we show that a travelling wave of the form
$u=sin(\xi)$, $\xi\equiv\omega_{0}t+kx$ can be a solution of (11) for an appropriate choice of the parameter $k$.

So we compute the time derivatives

and the spatial ones as,

Then if we substitute in (11) we get

to which we obtain two solutions: the trivial one if $\mu=0$ and another one with wave number $k=\sqrt{2}$ ,

which is a solution to our system independently of the value of $\varepsilon$.

Next we can analyse the case where we take into consideration a 3-dimensional version of the previous coupling scheme of the system. By using a lattice coupling instead of the 1-dimensional represented in figure (1), the equation (7) stands as

as done previously in the 1-dimensional case we can add and subtract the $2\mu u_{i,j,k}$ term which gives us

so as seen previously in equation (9) the ODE becomes a PDE as

being $\nabla^{2}$ the Laplacian operator defined as $\nabla^{2}=\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}+\frac{\partial^{2}}{\partial z^{2}}$ .

Similarly to equation (14), we find that a solution in the form $u=sin(\xi),\xi=\omega_{0}t+kx$ in equation (18) leads to $6\mu sin(\xi)+\mu(-3k^{2}sin(\xi))=0$, so we recover the exact same solutions $\mu=0$ and $k=\sqrt{6/3}=\sqrt{2}$ as in equation (15).

Moreover, we can also analyse the spatially homogeneous case in equation (11) by dropping the second spatial derivative in its last term, which leads to

If we propose again an analytic solution in the form of a travelling wave but without the spatial term by $u=sin(\Omega t)$ we get

which corresponds to the case where $\delta=1/\Omega^{2}$ and $\Omega=\sqrt{\omega_{0}-2\mu}$ independently of the value of the parameter $\varepsilon$.

Finally, with the aim of finding an asymptotic frequency we linearize equation (5) and we take $\varepsilon<<1$ and $\delta<<1$ so we can neglect the dissipative term. Similarly, we can add to equation (5) a coupling between its elements via a Laplacian matrix $L=[L_{ij}:i,j=1,...,n]$ for the graph $\mathcal{G}$. In that sense, $L_{ij}=a_{ij}$ if $i\neq j$ and $L_{ii}=\sum_{j=1}^{N}a_{ij}$. In this way we obtain the following equation

which can be rewritten as

where $A_{i}=\sum_{j=1}^{N}a_{ij}$ for $i=1,...,N$. In the case that we have $a_{ij}=1$, corresponding to a completely connected network of oscillators, the system can be written as

Supposing a synchronised regime has been established, we can have a synchronization frequency determined by

Apart from this, this reduction suggests that the synchronised oscillatory behaviour happens for sufficiently small values of $\mu$ and also when

If this wasn’t the case we would have an exponentially growing behaviour in our system. Otherwise, we can observe that unless all $a_{ij}$ are equal we cannot extract these conclusions. So, to analyse the general case where not all $a_{ij}$ are equal, we could use for instance certain methods that introduce perturbations to the Fokker-Planck equation [21]. Nonetheless, such analysis is outside the scope of this paper.

In the next section we develop a numerical method based on a finite differences algorithm techniques in order to numerically check the analytic solution that we have found. Finally we interpret those results in the context of the synchronised quorum of genetic clocks experiments in synthetic biology [5].

In order to find a solution to our system we develop a numerical simulation based on a finite differences method [22]. According tho this we start with the generic PDE

which is then transformed to lead to our equation (11). We begin by discretizing our temporal domains $[0,T]$ and the spatial ones $[0,L]$ so

and

In this manner we can introduce for a distributed set of points the spacing constants $\Delta t$ y $\Delta x$ within a mesh to get

Then the solution $u(x,t)$ is obtained for the points in the mesh, were we introduce the notation $u_{i}^{n}$, which approximates the exact solution in the points $\left(x_{i},t_{n}\right)$ for $i=0,\ldots,N_{x}$ y $n=0,\ldots,N_{t}$.

The next step is to change the derivatives by finite differences according to

and we introduce the operator notation for the finite differences as

and for the spatial case

We do the same with first order time derivatives obtaining the algebraic version of our PDE from equation (11) as

That equation can be interpreted as a moving stencil which involves the values in the four pointsfrom its neighbourhood $u_{i}^{n+1},u_{i\pm 1}^{n},u_{i}^{n},$ y $u_{i}^{n-1}$, as we can see in the figure 3.

In the same manner we have discretised the derivatives we need to discretize the initial conditions,

that can be written in algebraic notation as

leading to

and for the other initial condition

Finally we propose a recursive algorithm where $u_{i}^{n+1}$ equals to the sum of the following $A,...,J$ terms as:

where $A,...,J$ parameters are determined in our code (see section Supplementary Info.: Numerical Methods Code).

For the equation (40) we have a problem when $n=0$ since the formula for $u_{i}^{1}$ involves $u_{i}^{-1}$, which is not defined within the mesh. However, we can use the initial condition in (38) for $n=0$ to eliminate $u_{i}^{-1}$.

Finally our algorithm can be summarised as:

1. 1.

   Calculate $u_{i}^{0}=I\left(x_{i}\right)$ for $i=0,\ldots,N_{x}$.

2. 2.

   Calculate $u_{i}^{1}$ from (40) and take $u_{i}^{1}=0$ in the boundary points $i=0$ e $i=N_{x}$ for $n=1,2,\dots,N-1$

3. 3.

   For each time step $n=1,2,\dots,N_{t}-1$ apply (40) to find $u_{i}^{n+1}$ for $i=1,\ldots,N_{x}-1$ and take $u_{i}^{n+1}=0$ for the boundary points $i=0,i=N_{x}$.

This algorithm consists essentially in a moving stencil of finite difference elements through all the points in the mesh. Next we show the plotted results in Figure (4).

The synchronised quorum of genetic clocks synthetic system, which exhibits a range of oscillatory behaviours after a critical density of bacteria is reached, can be mainly framed by a network architecture that consists of an activator that activates its own repressor. Such motifs have been found to exist as regulatory modules in many genetic oscillators and circadian clock networks, studied both in vivo and by mathematical models [23].

Here we will focus on the compared analysis between our model, based on a system of spatially extended coupled oscillators, and some experimental realisations of the synchronised quorum system, such as [5][6]. These kinds of experimental designs make use of the previous repressilator-like architecture, additionally coupled to a system that synthesises a molecule known as an autoinducer (AI). This AI has the property that can diffuse through the cell membrane and allows cell-to-cell communication across the bacterial colony, in analogy to the Hooke-like coupling mechanism mediated by the parameter $\mu$ in our model.

The general solution that we found in (15) corresponds to a travelling wave, which can be described as a periodic function in one-dimensional space that moves with constant phase-velocity determined by the dispersion relation $v_{\mathrm{p}}=\omega_{0}/\sqrt{2}$, where $\omega_{0}$ corresponds to the angular frequency parameter selected within equation (11). This fact has the implication that the wave travels undistorted since $v_{\mathrm{p}}=v_{\mathrm{g}}$, being the group velocity defined by $v_{\mathrm{g}}\equiv\partial\omega/\partial k$.

The solution which was obtained by imposing a coupling mechanism in between elements within the discretised version of the hybrid Rayleigh-van der Pol equation (5) and then taking the continuous limit, does not depend on the coupling constant $\mu$. This might seem a counter-intuitive property of the system, nonetheless it seems to be a consequence of the robustness of the limit cycle under perturbations. Moreover, it makes sense from an evolutionary perspective, where one expects that the speed of propagation of a certain signal does not depend on the detailed features of the coupling, thus providing a structural advantage. Furthermore, we can also notice that the solution in equation (11) is also independent of the parameter $\varepsilon$, meaning that the limit cycle is unconstrained by the damping mechanism.

In addition, we can observe that switching the unitary term in equation (11) by an extra parameter $A$, such that

also allows for a more general travelling wave solution in the form of $u=Asin(\xi)$, $\xi\equiv\omega_{0}t+\sqrt{2}x$, besides the trivial solutions $\mu=0$ and $A=0$. As a consequence, by selecting the parameters $\omega_{0}$ and $A$ within equation (41), one can recover a family of harmonic waves with varying angular frequencies and maximum amplitudes constrained to a fixed wavenumber.

Similarly, a wide range of spatio-temporal dynamics of the synchronised oscillators can be accomplished by varying some experimental parameters, such as the diffusion of the AI molecule, which turns out to modulate the amplitude and the period of the wave. Thus, for a high dissipation rate of AI both the amplitude and period increase. However, this result which is consistent with ‘degrade-and-fire’ oscillations observed previously in intracellular oscillators [24], is not recovered in our model since its amplitude is ultimately decoupled from the frequency.

Additionally, we can note that the experimental results show a transition from an unsynchronised to a synchronised regime as the strength of coupling increases, caused by an increase in cell density. This characteristic property is recovered in our model since the conditions for the emergence of the limit cycle are fulfilled once we take the continuum limit, allowing the number of oscillators to become infinite while the coupling remains local.

Finally, since the speed of wave front propagation in our system does not depend on the coupling parameter $\mu$, in contrast to the experimental case where $v\sim\sqrt{D}$ (being $D$ related to the diffusive coefficient of AI), a subsequent work will focus on studying the synchronization of other diffusive-like coupling mechanisms, where synchronization is a consequence of the dissipation in the coupling along with cases where there is an interaction between the inherent damping in the system and the coupling.

In relation to the topics discussed in this paper, there is an issue regarding the emergence of the harmonic travelling wavefront solution that may be further clarified. The shape of the dispersion relation that we found, which turns out to be independent of the coupling proportionality constant between oscillators and also of the damping coefficient, may allude to an underlying robustness mechanism that needs to be interpreted. A possible way to do so would be in the context of biological soliton-like structures, defined as self-reinforcing solitary waves that travel at constant speed and without changing its shape [25]. One procedure to asses the existence of this kind of structures would be to check out if our solution does not obey the superposition principle, which implies that the travelling wave propagates unperturbed in the presence of other wave structures.

In general terms, whilst analytically solvable examples of coupled-oscillator models may be limited in each biological topic under study, the continuum-limit picture itself, i.e the one-oscillator extended to a collective ensemble preserving the local coupling in a self-consistent manner, may remain useful for qualitative understanding the spatio-temporal dynamics of signal or energy propagation within large groups of coupled oscillators.