# A shape-driven reentrant jamming transition in confluent monolayers of
synthetic cell-mimics

Pragya Arora  Chemistry and Physics of Materials Unit, Jawaharlal Nehru Centre
for Advanced Scientific Research, Jakkur, Bangalore - 560064, INDIA  Souvik
Sadhukhan  Tata Institute of Fundamental Research, Hyderabad - 500046, INDIA
Saroj Kumar Nandi  Tata Institute of Fundamental Research, Hyderabad - 500046,
INDIA  Dapeng Bi  Department of Physics, Northeastern University, Boston, MA
02115, U.S.A.  A K Sood  Department of Physics, Indian Institute of Science,
Bangalore- 560012, INDIA  International Centre for Materials Science,
Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur, Bangalore -
560064, INDIA  Rajesh Ganapathy  International Centre for Materials Science,
Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur, Bangalore -
560064, INDIA  School of Advanced Materials (SAMat), Jawaharlal Nehru Centre
for Advanced Scientific Research, Jakkur, Bangalore - 560064, INDIA

###### Abstract

Many critical biological processes, like wound healing, require confluent cell
monolayers/bulk tissues to transition from a jammed solid-like to a fluid-like
state. Although numerical studies anticipate changes in the cell shape alone
can lead to unjamming, experimental support for this prediction is not
definitive because, in living systems, fluidization due to density changes
cannot be ruled out. Additionally, a cell's ability to modulate its motility
only compounds difficulties since even in assemblies of rigid active
particles, changing the nature of self-propulsion has non-trivial effects on
the dynamics. Here, we design and assemble a monolayer of synthetic cell-
mimics and examine their collective behaviour. By systematically increasing
the persistence time of self-propulsion, we discovered a cell shape-driven,
density-independent, re-entrant jamming transition. Notably, we observed cell
shape and shape variability were mutually constrained in the confluent limit
and followed the same universal scaling as that observed in confluent
epithelia. Dynamical heterogeneities, however, did not conform to this
scaling, with the fast cells showing suppressed shape variability, which our
simulations revealed is due to a transient confinement effect of these cells
by their slower neighbors. Our experiments unequivocally establish a
morphodynamic link, demonstrating that geometric constraints alone can dictate
epithelial jamming/unjamming.

Confluent cell monolayers and tissues are a maximally crowded environment; the
cell packing fraction is nearly unity. And yet, remodeling and repair occur,
and cancer cells migrate from a tumor and invade distal sites. These processes
require tissue/cell collectives to flow locally [1, 2, 3, 4, 5, 6, 7, 8], and
the jamming-unjamming transition, akin to the one seen in inert particle
assemblies [9], provides a pathway. In inert particle assemblies, the
unjamming transition is density-driven [10, 11], and there is evidence of such
a transition during embryonic morphogenesis [12] and cancer invasion [13, 14,
15]. But cells, unlike inert particles, can deform to surmount the
constraining effects of crowding [6, 7, 8, 16]. In fact, in the vertex model
of confluent epithelia [17, 18], the competing effects of cell contractility
and cell-cell adhesion lead to a density-independent but cell shape change
driven unjamming transition [19]. In the jammed state, the cells have a more
regular shape than in the fluid state, where they are more elongated. A
subsequent model that included cell motility showed the jamming transition
could be driven by additional factors, but the qualitative shape-based nature
remained unchanged [20]. Cell shape is proving to be a structural marker of
unjamming [6, 7, 8, 16, 21], and this is the case even when alignment
interactions between neighboring cells are present [22, 23, 24, 25]. Cell
shape-mediated unjamming is now implicated in the pathophysiology of asthma
[26] and tumor progression [15, 27].

Besides shape, cells in a collective also show substantial shape variability,
which, until recently, was dismissed as noise. Atia et al., [28], observed
that across vastly different epithelial systems, the cell shape distribution
became progressively less skewed as the system jammed. When scaled
appropriately, these distributions collapsed to a universal distribution known
as the $k-$gamma distribution [29]. Remarkably, across these different
systems, the shape variability and the cell shape follow a simple linear
relationship; the larger the cell aspect ratio, the greater the shape
variability. These results suggest that, like in inert particle packings [29,
30], geometric constraints take center stage in tissue jamming/unjamming.

Tissues, however, are complex, and besides cell shape changes, jamming-
unjamming can be driven by cell division, apoptosis, extrusion, cell size
changes, and motile topological defects [31]. These processes may operate
independently or work in cohorts [28], and being impossible to suppress, there
is no consensus if cell shape changes alone drive unjamming [6, 8].
Additionally, cells can regulate their motility, which is another critical
parameter governing glass/jamming physics [32, 33]. Even for the simpler case
of rigid self-propelled particles, tuning the activity at constant density
results in non-intuitive behavior, such as re-entrant dynamics [34, 35, 36],
and it is unknown whether similar physics is at play in cell collectives,
since, here, the means to tune activity systematically are not available.
Synthetic models that embody the key features of confluent epithelia can help
bridge the divide between theory/numerical predictions and experiments on
living systems. They also have the potential to single out the role of
geometric constraints in tissue jamming/unjamming. In fact, a recent theory
inspired by [28] posits that the relationship between shape variability and
cell shape is purely a mathematical property of a confluent monolayer of
closed-loop objects [37], but this idea still needs to be tested in a
synthetic model system.

### Making deformable cell-mimics with tunable activity

Here, we create synthetic active cell-mimics, assemble them into monolayers,
and examine their collective behavior. At the single-cell level, our mimics
possess just two features: deformability and tunable activity. Existing
deformable active matter systems - centimeter-sized bots confined within
paper/metal rings [38] and active colloids confined within vesicles [39] do
not serve our purpose. In the former system, each cell is itself large
($\approx 30$ cm in diameter), making it a challenge to study collective
behavior; with the latter system, methods to create assemblies of these and
regulate their activity are unavailable.

Our cell-mimics instead are flexible paper rings ($3$ cm in diameter) that
enclose a monolayer of 3D-printed granular ellipsoids rendered chiral active
by vertical vibration. These ellipsoids experience both an active torque and
force under such a driving and perform circle active motion with the
handedness of the trajectory, i.e., clockwise $(+)$ or counterclockwise $(-)$,
manually determined during particle placement on our shaker apparatus [42]
(Fig. 1a top panel & see Materials and Methods and Fig. S1-S4). Our
preliminary observations guided us to leverage chiral activity rather than
achiral activity to make our cells self-propelled. On confining achiral polar
active ellipsoids within the paper rings, the interplay of membrane curvature
and particle orientation often resulted in their accumulation at diametrically
opposite ends with their polarity pointing outwards [38, 40, 41]. Besides
resulting in low cell motility, this preferential accumulation led to a
nonuniform cell stiffness (Supplementary Fig. S4 and Supplementary Movie 1).
In contrast, when persistent active torques are also present, particles hug a
wall with their propulsion axis parallel to it and with the direction of
motion along the wall set by the handedness of the activity (Fig. 1a bottom
panel) [43]. Supplementary Movie 2 shows the dynamics of a granular cell with
$N=20$ chiral ellipsoids of the same handedness, $(-)$, within (Figure 1c(i)).
Unlike with the achiral ellipsoids, the particles now uniformly decorated the
membrane interior with the polarized particle current acting like a ``dynamic
internal skeleton". However, this current only caused the paper ring to spin
with a well-defined handedness, but we observed little translation.

To tune the motility of the cell-mimics, we changed the handedness of a few
ellipsoids inside the membrane from $(-)$ to $(+)$ (red ellipsoids in Fig.
1c(ii)-(iv)). We observed counter-propagating particle currents that
suppressed cell spin and made it motile (see Supplementary Movie 2 and
Supplementary Table I). Thus, even when $N=N_{+}+N_{-}$ is held fixed at
twenty particles, changing the magnitude of chirality of the cell interior,
$|\chi_{\text{cell}}|$, helps tune cell activity. Here,
$\chi_{\text{cell}}=\frac{N_{+}-N_{-}}{N_{+}+N_{-}}$, and $N_{+}$ and $N_{-}$
are the number of $(+)$ and $(-)$ ellipsoids, respectively. We note that the
handedness of cell spin depends on the sign of $\chi_{\text{cell}}$, but the
internally generated active force depends only on its magnitude. For example,
the green and yellow cells depicted in Fig. 1d spin in opposite directions due
to an excess of $(-)$ ellipsoids in one and $(+)$ ellipsoids in the other;
however, the active force is determined solely by $|N_{+}-N_{-}|$, which is
the same for both.

We quantified the activity of isolated cells for different
$|\chi_{\text{cell}}|$ values by working in the low area fraction limit of
cells, $\phi<1\%$. We measured the persistence time, $\tau_{p}$ \- the time at
which the cells' mean-squared displacement crossed over from ballistic to
diffusive dynamics (Supplementary Fig. S5a), and also the average cell speed,
$v$ (Supplementary Fig. S5b) [36]. A larger value of $\tau_{p}$ signifies a
stronger departure from equilibrium, i.e., greater activity. When the cell
interior was gradually changed from enantiopure ($|\chi_{\text{cell}}|=1)$ to
racemic ($|\chi_{\text{cell}}|=0)$, both $\tau_{p}$ and the persistence
length, $l_{p}=v\tau_{p}$, increased in a systematic manner (Fig. 1e,
Supplementary Table II). Recent studies have identified $\tau_{p}$ as a
crucial parameter governing active glass physics [34, 33, 44, 36], and our
approach allows for systematically tuning it. Interestingly, the active
torques and forces exerted by the particles on the membrane interior also
influenced the shape variability of isolated cells. We used the standard
deviation, $SD$, of an isolated cells' aspect ratio, $AR_{s}$, as a measure of
shape variability, and this is shown in Fig. 1f for different
$|\chi_{\text{cell}}|$ values. In the limit of an enantiopure cell interior,
the persistent torques due to the unidirectional particle current at the
boundary suppressed cell shape fluctuations, resulting in small $SD$, while on
approaching a racemic cell interior, the more frequent reorganization of the
counter-propagating particle currents made it floppy (large $SD$, see
Supplementary Movie 2).

### Cell-mimic collective

Before we can ascertain if our synthetic cells in the dense limit embody the
crucial features of confluent epithelia, another parameter, namely, the net
chirality, $\chi_{\text{Sys}}$, of the synthetic cell assemblies must be
adjusted. Here,
$\chi_{\text{Sys}}=\frac{|N_{\text{CW}}-N_{\text{CCW}}|}{N_{\text{CW}}+N_{\text{CCW}}}$,
where $N_{\text{CW}}$ and $N_{\text{CCW}}$ are the number of clockwise and
counterclockwise spinning cells. We observed that $\chi_{\text{Sys}}$ had a
qualitative effect on the dynamics; for example, in a nearly confluent
assembly ($\phi\approx 0.94$) of clockwise spinning cells, i.e.,
$\chi_{\text{Sys}}=1$ (green cell in Fig. 1d), there was an emergent edge
current (Supplementary Movie 3) like those seen in other confined chiral
active matter systems [45, 46]. However, except in very rare cases and
specific cell lines [47, 48], most experiments examining the jamming-unjamming
of confluent epithelia do not report such edge currents [1, 2, 26, 28],
suggesting that these assemblies do not possess an overall chirality. To make
contact with these studies, we, therefore, set $\chi_{\text{Sys}}=0$ by having
an equal number of clockwise (green) and counterclockwise (yellow) spinning
cells on the plate for all values of $\phi$ and $\tau_{p}$ studied (see Fig.
1f for $\phi=0.92$). The edge flows were absent in these assemblies.
(Supplementary Movie 4).

### A re-entrant jamming transition mediated by cell shape

Next, by analyzing the dynamical trajectories of the cell centers (see
Materials and Methods), we charted the relaxation dynamics of the granular
cell assemblies for different values of $\phi$ and $\tau_{p}$ (Fig. 2a-c). We
determined the structural relaxation time, $\tau_{\alpha}$, as the time at
which the self-intermediate scattering function $F_{s}(q,t)$ decayed to
$\frac{1}{e}$ (Supplementary Fig. S6). Here, $q$ is the wavevector and was
chosen to correspond to the inverse of the cell diameter. While for a fixed
value of $\tau_{p}$, increasing $\phi$ results in dynamical slowing down, as
expected, the effects of increasing $\tau_{p}$ at fixed $\phi$ is more subtle
(Fig. 2b). For $0.4\leq\phi\leq 0.88$, an increase in $\tau_{p}$ sped up the
dynamics, and over a narrower window $0.4\leq\phi\leq 0.7$, we observed a weak
tendency of the cells to form transient clusters (Fig. 2a and Supplementary
Movie 5) This clustering without attractive interactions is a cardinal feature
of active matter [49, 50] and was evident as a shoulder in the distribution of
Voronoi areas (Supplementary Fig. S7). Notably, over a range of densities
close to confluence ($0.88\leq\phi\leq 0.92$), increasing $\tau_{p}$ resulted
in a re-entrant behavior: structural relaxation was fastest at an intermediate
value of $\tau_{p}$ (horizontal bar in Fig. 2b). This re-entrant behavior is
evident in the cell trajectories, with the cells being strongly caged (glass-
like) for both small and large values of $\tau_{p}$ and ergodic (fluid-like)
for an intermediate value (vertical bar corresponding to $\phi\approx 0.91$ in
Fig. 2c, see Supplementary Movie 6).

On increasing the strength of particle attraction in dense assemblies of hard
passive particles, the competing effects of particle crowding and attraction
often lead to re-entrant dynamics [51, 52]. Similar behavior is seen in dense
assemblies of hard active particles on increasing $\tau_{p}$ [32, 34, 35, 36],
with the difference being that the attraction itself is activity-mediated, and
its strength is proportional to $\tau_{p}$. However, unlike in hard
passive/active particle systems, where the re-entrant behavior begins to
manifest even at moderate densities [52, 34, 36], we found it only on nearing
confluence, which suggests that cell shape changes may have a vital role here.
We calculated the average cell aspect ratio $\overline{AR}$ of the assemblies
for a range of $\phi$ and all values of $\tau_{p}$ studied (Fig. 2d).
Remarkably, over the window of $\phi$ values where we observed a re-entrant
behavior (yellow-filled symbols in Fig. 2d), it is at the intermediate value
of $\tau_{p}$ where the system was fluid-like that $\overline{AR}$ is largest
(horizontal bar in Fig. 2d). Even in our synthetic system, cell shape governs
structural relaxation: assemblies of elongated cells relax faster (Fig. 2f)
[15, 19, 20, 26]. Indeed, for moderate densities ($\phi=0.75\text{, }0.60$),
not only does the lack of re-entrant dynamics coincide with the absence of a
non-monotonicity in $\overline{AR}$ with $\tau_{p}$ (Fig. 2d), the
correspondence between $\tau_{\alpha}$ and $\overline{AR}$ is also weaker
(Fig. 2f).

The above findings suggest that shape variability may also bear the imprint of
the dynamics of these assemblies [28]. We quantified the shape variability of
our cell assemblies via the standard deviation of the aspect ratio $SD(AR)$.
Strikingly, $SD(AR)$ mirrors the behavior of $\overline{AR}$ for different
values of $\phi$ and $\tau_{p}$ (Fig. 2e). Additionally, similar to
$\overline{AR}$, there is a correlation between $SD(AR)$ and $\tau_{\alpha}$,
with faster relaxation observed in assemblies with greater shape variability
(Fig. 2g).

We can now explain the observed re-entrant behavior. For small $\tau_{p}$
values (large $|\chi_{\text{cell}}|$), besides the low cell motility, each
cell also resists significant shape fluctuations (Fig. 1f) and adopts a
circular shape (small $\overline{AR}$). Since disk assemblies freeze into
glasses at densities smaller than elongated particle ones as they require
fewer constraints [8], the system is glassy near confluence (Supplementary
Movie 6). At the largest value of $\tau_{p}$ (i.e., $|\chi_{\text{cell}}|=0$),
individual cell motility is large, and the cells are also floppier (Fig. 1f),
but so is the activity-mediated adhesion between cells. Near confluence, this
cell-cell adhesion brings down $\overline{AR}$ and $SD(AR)$, and the system is
again glassy. However, these effects compete for an intermediate $\tau_{p}$
value: individual cell motility and shape variability are reasonably large,
but adhesion is weak, making the assembly fluid-like. As far as we know, this
is the first observation of a shape-mediated re-entrant jamming.

### Shape and shape variability in confluent cell-mimic monolayers

The strikingly similar behavior of $\overline{AR}$ and $SD(AR)$ for different
$\phi$ and $\tau_{p}$ values (Fig. 2d & e) suggests these two quantities are
interdependent. Indeed, for our granular cells, $SD(AR)$ scales linearly with
$\overline{AR}$ like that observed in vastly different epithelial systems
(Fig. 3a) [28, 37]. We now checked whether this linear scaling was simply a
result of the probability distribution function ($PDF$) of the aspect ratio
being universal and a $k-$gamma distribution [28, 29]. Figure 3b shows the
$PDF$ for different $\tau_{p}$ values for $\phi=0.92$. The $PDF$ is broad with
a large positive skew for the unjammed (fluid-like) state, which, here, is
observed at an intermediate value of $\tau_{p}=20$ s. For the more jammed
states, corresponding to other $\tau_{p}$ values, the $PDF$s are narrower and
with a smaller skew. This behavior is also seen for other values of $\phi$
near confluence (Supplementary Figs. S8-S11). Moreover, the $PDF$s at
different $\tau_{p}$ values collapse onto a single universal distribution on
rescaling $AR$ to $x=\frac{(AR-1)}{(\overline{AR}-1)}$, which ensures that the
scaled $PDF$s start at the origin and have a mean of unity. Within
experimental uncertainty, the scaled $PDF$s are equally well-fit by the
$k-$gamma distribution (1) with the value of $k\approx 2.51$ [28], as well as
the recently posited distribution for a confluent assembly of closed-loop
objects (2), which is only nearly universal [37] (Supplementary Figs. S12).

$P(x,k)=\left[k^{k}/\Gamma(k)\right]x^{k-1}\exp[-kx]$ (1)
$P(AR)=\frac{1}{\mathcal{N}}\left(AR+\frac{1}{AR}\right)^{3/2}\left(1-\frac{1}{AR^{2}}\right)e^{-\alpha\left(AR+\frac{1}{AR}\right)},$
(2)

In the first expression, $\Gamma(k)$, is the Legendre gamma function; in the
second, $N$, is the normalization constant, and $\alpha$ is a system-specific
parameter predicted to correlate with dynamics. The PDF of $x$ obtained from
(2) becomes nearly independent of $\alpha$. For all values of $\phi$ close to
confluence and all $\tau_{p}$ values, $k$ hovers between 2.5$-$2.8 like that
observed in confluent epithelia [28].

Equation (2) follows from a mean-field theory [37], which also anticipates
that for a confluent monolayer of closed-loop objects, which applies to both
living and our synthetic cells, $SD(AR)=0.71\overline{AR}-0.75$ (black line in
Fig. 3a). The agreement between our experiments and theory is excellent even
for densities far from confluence ($\phi=0.74$) (inset to (Fig. 3a)), although
the agreement becomes weaker for smaller $\phi$ values (Supplementary Fig.
S8). Notably, we also find that $\alpha$ scales linearly with
$\text{log}(\tau_{\alpha})$, as predicted [37], indicating a strong link
between structure and dynamics (Supplementary Fig. S13).

###  .1 Dynamical heterogeneities, cell shape and shape variability

Although it is widely recognized that the dynamics of dense passive/active
liquids and glasses is spatiotemporally heterogeneous and consist of domains
with different relaxation times [53, 54, 34, 36], the correlation between cell
aspect ratio and dynamics, as we report here so far, and like in previous
studies [19, 26, 28], has been arrived at by system-wide averaging. In fact,
recent experiments on confluent cell monolayers and 3D tissues show that such
a correlation is also present locally, i.e., cells in the faster relaxing
regions are more elongated than those in the slower relaxing ones [8, 15].
Additionally, since at the system level, $SD(AR)$ and $\overline{AR}$ are
mutually constrained (Fig. 3a), it is reasonable to expect that cell shape
variability in these faster regions is also correspondingly larger, although
this remains untested.

To this end, we followed standard procedures to quantify dynamical
heterogeneities [53, 54] and first identified the top $10$% fastest and
slowest cells over the cage-rearrangement time, $t^{*}$ (Supplementary Fig.
S14). In Fig. 4a, we show the fastest cells (red circles) overlaid on the cell
displacement maps also generated over $t^{*}$ for different $\tau_{p}$ values
for $\phi=0.92$. As is typical of deeply supercooled liquids, these cells are
spatially clustered, which is the case with slow cells as well, and further,
these fast cells within a cluster move together largely as a flock. Also, the
heterogeneity is most pronounced at an intermediate value of $\tau_{p}$
(Supplementary Movie S7). Next, we determined the mean and standard deviation
of the aspect ratio of the cells within these fast and slow clusters. Notably,
fast cells have a larger $\overline{AR}$ than slow cells (Fig. 4b). However,
while for the slow cells, the growth of $SD(AR)$ with $\overline{AR}$ follows
the universal scaling, for the fast cells, surprisingly, it lies below the
predicted line, implying that these cells show reduced shape variability
(solid symbols in Fig. 4c).

To strengthen these findings and examine their generality, we performed
simulations of the Vertex model (VM) [18] and the Cellular Potts model (CPM)
[55], which are highly effective representations of confluent epithelia. The
larger system size and superior temporal statistics of the simulations also
allowed quantifying the $PDF$ for the fast and slow cells. We observed that
while the $PDF$ of the slow cells overlapped with that of the entire system,
the $PDF$ for the fast cells had a larger skew, as is expected for a locally
unjammed region of the system (Fig. 4d and Supplementary Fig. S15). It is
noteworthy that even in the simulations, the growth of $SD(AR)$ with
$\overline{AR}$ for slow cells followed the predicted scaling, while for fast
cells, $SD(AR)$ was lower than anticipated (half-filled symbols in Fig. 4c and
Supplementary Fig. S16).

The fact that most-mobile cells do not conform to the universal scaling is, in
all likelihood, because the scaling comes from a mean-field theory, and hence,
fluctuations are ignored. Mean-field approaches fail to capture the physics of
dynamical heterogeneities even in conventional dense liquids and glasses [53].
Nonetheless, to gain insights into our observations, we made the simplifying
assumption that on a time scale comparable to $t^{*}$, fast cells are
essentially confined within a frozen exterior of slower cells due to their
different relaxation rates. To mimic this effect numerically, we first allowed
the entire system to equilibrate over many times $\tau_{\alpha}$, and then
froze all but a cluster of $n$ cells in the system interior (shown in cyan in
the inset to Fig. 4e). We performed simulations at different temperatures to
let the system access different values of $\overline{AR}$. Trivially, when all
but one cell is frozen ($n=1$), there can be no shape variability, and
$SD(AR)=0$. Indeed, on increasing $n$, $SD(AR)$ increased systematically and
followed the predicted scaling for $n\geq 12$. This observation nicely
demonstrates that the subdued shape variability of the fast cells is an
outcome of the temporary confinement imposed by their slower neighbors.

###  .2 Discussion

Our confluent synthetic cell monolayers are an oversimplification of confluent
epithelia. Unlike in the latter, no explicit cell-cell adhesion exists,
although an effective activity-mediated attraction emerges here. Also, our
cells have a fixed perimeter, which is not the case with live cells. The area
constraint is weaker in our system than confluent epithelia, but it still
exists because of a dynamic internal skeleton. Nonetheless, in the context of
jamming, our synthetic cell collective behaves like a living one, thus proving
that geometric constraints alone can govern the dynamics of deformable active
matter systems, be they inert or living. The exquisite control offered by our
synthetic model system allowed us to discover a re-entrant jamming transition
mediated solely by cell shape changes. Additionally, by dissecting dynamical
heterogeneities, we have shown that transient confinement of the fast cells by
their slower neighbors suppresses shape variability and leads to the violation
of mean-field predictions. Cell shape, on the other hand, correlates with both
global and local dynamics.

On the experimental front, the model system developed here allows for precise
and easy tuning of many critical physical parameters like the chirality of the
assembly, individual cell membrane stiffness, the nature of activity within
the membranes, and even the inter-cell friction. Although each of these
parameters has been identified to substantially impact the form and function
of the cell collective [56, 57, 58], these are impossible to control
systematically in real systems, and some of these, like inter-cell friction,
are yet to be incorporated even within the models of confluent epithelia.
Through our preliminary experiments with magnetic membranes, we can already
modify individual cell properties in real time. The experimental advance made
here now makes it possible to dissect the role of heterogeneous cell
properties on collective behavior [59] and in fundamental processes like cell
sorting [60] from a purely geometric/mechanical perspective [61, 62, 63].

##  I Author Contributions

P.A.: Conceptualization, Methodology (experiments), Software, Validation,
Formal Analysis, Investigation, Data Curation, Visualization, Writing -
Original Draft, Writing - Review & Editing. S.S.: Methodology (simulations),
Software, Validation, Formal Analysis, Data Curation, Writing - Review &
Editing. S.K.N.: Methodology (simulations), Software, Validation, Data
Curation, Formal Analysis, Writing - Review & Editing, Supervision
(simulations), Funding Acquisition. D.B.: Formal Analysis, Validation, Writing
- Review & Editing. A.K.S.: Validation, Writing - Review & Editing. R.G.:
Conceptualization, Methodology (experiments), Validation, Investigation,
Formal Analysis, Visualization, Writing - Original Draft, Writing - Review &
Editing, Supervision (experiments), Project Administration, Funding
Acquisition.

##  II Acknowledgements

P.A. thanks the Jawaharlal Nehru Centre for Advanced Scientific Research,
Bangalore, INDIA for a research fellowship. S.S and S.K.N. acknowledge the
support of the Department of Atomic Energy, Government of India, under Project
Identification No. RTI 4007. SKN thanks SERB for grant via SRG/2021/002014.
A.K.S. thanks the Science and Engineering Research Board, Government of India
for the National Science Chair . R.G. thanks the Department of Science and
Technology, INDIA, for financial support through the SwarnaJayanti Fellowship
Grant (DST/SJF/PSA-03/2017-22).

##  III Data Availability

All relevant data are available from the corresponding author upon reasonable
request.

## References

  * [1] Trepat, X., Wasserman, M. R., Angelini, T. E., Millet, E., Weitz, D. A., Butler, J. P., & Fredberg, J. J. (2009). Physical forces during collective cell migration. Nature Physics, 5(6), 426-430.
  * [2] Angelini, T. E., Hannezo, E., Trepat, X., Marquez, M., Fredberg, J. J., $\&$ Weitz, D. A. (2011). Glass-like dynamics of collective cell migration. Proceedings of the National Academy of Sciences, 108(12), 4714–4719.
  * [3] Sadati, M., Qazvini, N. T., Krishnan, R., Park, C. Y., $\&$ Fredberg, J. J. (2013). Collective migration and cell jamming. Differentiation, 86(3), 121–125.
  * [4] Schoetz, E.-M., Lanio, M., Talbot, J. A., $\&$ Manning, M. L. (2013). Glassy dynamics in three-dimensional embryonic tissues. Journal of The Royal Society Interface, 10(89), 20130726.
  * [5] Garcia, S., Hannezo, E., Elgeti, J., Joanny, J.-F., Silberzan, P., $\&$ Gov, N. S. (2015). Physics of active jamming during collective cellular motion in a monolayer. Proceedings of the National Academy of Sciences, 112(50), 15314–15319.
  * [6] Oswald, L., Grosser, S., Smith, D. M., & Käs, J. A. (2017). Jamming transitions in cancer. Journal of Physics D: Applied Physics, 50(48), 483001.
  * [7] Atia, L., Fredberg, J. J., Gov, N. S., & Pegoraro, A. F. (2021). Are cell jamming and unjamming essential in tissue development? Cells & Development, 168, 203727.
  * [8] Blauth, E., Kubitschke, H., Gottheil, P., Grosser, S., & Käs, J. A. (2021). Jamming in embryogenesis and cancer progression. Frontiers in Physics, 9, 666709.
  * [9] Liu, A. J., & Nagel, S. R. (1998). Jamming is not just cool any more. Nature, 396(6706), 21-22.
  * [10] Liu, A. J., & Nagel, S. R. (2010). The jamming transition and the marginally jammed solid. Annu. Rev. Condens. Matter Phys., 1(1), 347-369.
  * [11] Behringer, R. P., & Chakraborty, B. (2018). The physics of jamming for granular materials: a review. Reports on Progress in Physics, 82(1), 012601.
  * [12] Mongera, A., et al., (2018). A fluid-to-solid jamming transition underlies vertebrate body axis elongation. Nature, 561(7723), 401-405.
  * [13] Han, Y. L., et al., (2020). Cell swelling, softening and invasion in a three-dimensional breast cancer model. Nature Physics, 16(1), 101-108.
  * [14] Ilina, O., et al., (2020). Cell–cell adhesion and 3D matrix confinement determine jamming transitions in breast cancer invasion. Nature Cell Biology, 22(9), 1103-1115.
  * [15] Grosser, S., et al., (2021). Cell and nucleus shape as an indicator of tissue fluidity in carcinoma. Physical Review X, 11(1), 011033.
  * [16] Stephens, M. (2021). Elongated Cells May Unjam Cancers. Physics, 14, s19.
  * [17] Nagai, T., & Honda, H. (2001). A dynamic cell model for the formation of epithelial tissues. Philosophical Magazine B, 81(7), 699-719.
  * [18] Farhadifar, R., Röper, J. C., Aigouy, B., Eaton, S., & Jülicher, F. (2007). The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing. Current Biology, 17(24), 2095-2104.
  * [19] Bi, D., Lopez, J. H., Schwarz, J. M., $\&$ Manning, M. L. (2015). A density-independent rigidity transition in biological tissues. Nature Physics, 11(12), 1074–1079.
  * [20] Bi, D., Yang, X., Marchetti, M. C., $\&$ Manning, M. L. (2016). Motility-driven glass and jamming transitions in biological tissues. Physical Review X, 10(4), 041010.
  * [21] Trepat, X., $\&$ Sahai, E. (2018). Mesoscale physical principles of collective cell organization. Nature Physics, 14(7), 671-682.
  * [22] Malinverno, C., et al., (2017). Endocytic reawakening of motility in jammed epithelia. Nature Materials, 16(5), 587-596.
  * [23] Giavazzi, F., et al., (2018). Flocking transitions in confluent tissues. Soft matter, 14(18), 3471-3477.
  * [24] Wang, X., Merkel, M., Sutter, L. B., Erdemci-Tandogan, G., Manning, M. L., & Kasza, K. E. (2020). Anisotropy links cell shapes to tissue flow during convergent extension. Proceedings of the National Academy of Sciences, 117(24), 13541-13551.
  * [25] Paoluzzi, M., Angelani, L., Gosti, G., Marchetti, M. C., Pagonabarraga, I., & Ruocco, G. (2021). Alignment interactions drive structural transitions in biological tissues. Physical Review E, 104(4), 044606.
  * [26] Park, J.-A., Kim, J. H., Bi, D., Mitchel, J. A., Qazvini, N. T., Tantisira, K., … $\&$ Tschumperlin, D. J. (2015). Unjamming and cell shape in the asthmatic airway epithelium. Nature Materials, 14(10), 1040-1048.
  * [27] Gottheil, P., et al., (2023). State of cell unjamming correlates with distant metastasis in cancer patients. Physical Review X, 13(3), 031003.
  * [28] Atia, L., Bi, D., Sharma, Y., Mitchel, J. A., Gweon, B., A. Koehler, S., DeCamp, S. J., Lan, B., Kim, J. H., Hirsch, R., $\&$ others. (2018). Geometric constraints during epithelial jamming. Nature Physics, 14(6), 613-620.
  * [29] Aste, T., & Di Matteo, T. (2008). Emergence of Gamma distributions in granular materials and packing models. Physical Review E, 77(2), 021309.
  * [30] Boromand, A., Signoriello, A., Ye, F., O’Hern, C. S., & Shattuck, M. D. (2018). Jamming of deformable polygons. Physical Review Letters, 121(24), 248003.
  * [31] Saw, T. B. et al., (2017). Topological defects in epithelia govern cell death and extrusion. Nature, 544(7649), 212-216.
  * [32] Berthier, L., Flenner, E., & Szamel, G. (2017). How active forces influence nonequilibrium glass transitions. New Journal of Physics, 19(12), 125006.
  * [33] Nandi, S. K., Mandal, R., Bhuyan, P. J., Dasgupta, C., Rao, M., & Gov, N. S. (2018). A random first-order transition theory for an active glass. Proceedings of the National Academy of Sciences, 115(30), 7688-7693.
  * [34] Berthier, L. (2014). Nonequilibrium glassy dynamics of self-propelled hard disks. Physical Review Letters, 112(22), 220602.
  * [35] Klongvessa, N., Ginot, F., Ybert, C., Cottin-Bizonne, C., & Leocmach, M. (2019). Active glass: Ergodicity breaking dramatically affects response to self-propulsion. Physical Review Letters, 123(24), 248004.
  * [36] Arora, P., Sood, A. K., & Ganapathy, R. (2022). Motile topological defects hinder dynamical arrest in dense liquids of active ellipsoids. Physical Review Letters, 128(17), 178002.
  * [37] Sadhukhan, S., & Nandi, S. (2022). On the origin of universal cell shape variability in a confluent epithelial monolayer. Elife, 11, e76406.
  * [38] Boudet, J. F.,et al., (2021). From collections of independent, mindless robots to flexible, mobile, and directional superstructures. Science Robotics, 6(56), eabd0272.
  * [39] Vutukuri, H. R., et al., (2020). Active particles induce large shape deformations in giant lipid vesicles. Nature, 586(7827), 52-56.
  * [40] Paoluzzi, M., Di Leonardo, R., Marchetti, M. C., & Angelani, L. (2016). Shape and displacement fluctuations in soft vesicles filled by active particles. Scientific reports, 6(1), 34146.
  * [41] Peterson, M. S., Baskaran, A., & Hagan, M. F. (2021). Vesicle shape transformations driven by confined active filaments. Nature communications, 12(1), 7247.
  * [42] Arora, P., Sood, A. K., & Ganapathy, R. (2021). Emergent stereoselective interactions and self-recognition in polar chiral active ellipsoids. Science Advances, 7(9), eabd0331.
  * [43] Barois, T., Boudet, J.-F., Lintuvuori, J. S., $\&$ Kellay, H. (2020). Sorting and extraction of self-propelled chiral particles by polarized wall currents. Physical Review Letters, 125(23), 238003.
  * [44] Henkes, S., Kostanjevec, K., Collinson, J. M., Sknepnek, R., & Bertin, E. (2020). Dense active matter model of motion patterns in confluent cell monolayers. Nature Communications, 11(1), 1405.
  * [45] Van Zuiden, B. C., Paulose, J., Irvine, W. T. M., Bartolo, D., $\&$ Vitelli, V. (2016). Spatiotemporal order and emergent edge currents in active spinner materials. Proceedings of the National Academy of Sciences, 113(46), 12919-12924.
  * [46] Beppu, K., et al., (2021). Edge current and pairing order transition in chiral bacterial vortices. Proceedings of the National Academy of Sciences, 118(39), e2107461118.
  * [47] Duclos, G., Blanch-Mercader, C., Yashunsky, V., Salbreux, G., Joanny, J. F., Prost, J., & Silberzan, P. (2018). Spontaneous shear flow in confined cellular nematics. Nature Physics, 14(7), 728-732.
  * [48] Yashunsky, V., Pearce, D. J. G., Blanch-Mercader, C., Ascione, F., Silberzan, P., & Giomi, L. (2022). Chiral edge current in nematic cell monolayers. Physical Review X, 12(4), 041017.
  * [49] Fily, Y., & Marchetti, M. C. (2012). Athermal phase separation of self-propelled particles with no alignment. Physical Review Letters, 108(23), 235702.
  * [50] Cates, M. E., & Tailleur, J. (2015). Motility-induced phase separation. Annu. Rev. Condens. Matter Phys., 6(1), 219-244.
  * [51] Pham, K. N., et al., (2002). Multiple glassy states in a simple model system. Science, 296(5565), 104-106.
  * [52] Mishra, C. K., Rangarajan, A., & Ganapathy, R. (2013). Two-step glass transition induced by attractive interactions in quasi-two-dimensional suspensions of ellipsoidal particles. Physical Review Letters, 110(18), 188301.
  * [53] Berthier, L., $\&$ Biroli, G. (2011). Theoretical perspective on the glass transition and amorphous materials. Reviews of Modern Physics, 83(2), 587.
  * [54] Gokhale, S., Sood, A. K., $\&$ Ganapathy, R. (2016). Deconstructing the glass transition through critical experiments on colloids. Advances in Physics, 65(4), 363-452.
  * [55] Noppe, A. R., Roberts, A. P., Yap, A. S., Gomez, G. A., & Neufeld, Z. (2015). Modelling wound closure in an epithelial cell sheet using the cellular Potts model. Integrative Biology, 7(10), 1253-1264.
  * [56] Sato, K., Hiraiwa, T., & Shibata, T. (2015). Cell chirality induces collective cell migration in epithelial sheets. Physical Review Letters, 115(18), 188102.
  * [57] Fuhs, T., Wetzel, et al., (2022). Rigid tumours contain soft cancer cells. Nature Physics, 18(12), 1510-1519.
  * [58] Angelini, T. E., Dunn, A. C., Urueña, J. M., Dickrell, D. J., Burris, D. L., & Sawyer, W. G. (2012). Cell friction. Faraday Discussions, 156(1), 31-39.
  * [59] Li, X., Das, A., & Bi, D. (2019). Mechanical heterogeneity in tissues promotes rigidity and controls cellular invasion. Physical Review Letters, 123(5), 058101.
  * [60] Lawson-Keister, E., & Manning, M. L. (2021). Jamming and arrest of cell motion in biological tissues. Current Opinion in Cell Biology, 72, 146-155.
  * [61] Manning, M. L. (2023). Essay: Collections of deformable particles present exciting challenges for soft matter and biological physics. Physical Review Letters, 130(13), 130002.
  * [62] Graner, F., & Riveline, D. (2017). ‘The Forms of Tissues, or Cell-aggregates’: D'Arcy Thompson's influence and its limits. Development, 144(23), 4226-4237.
  * [63] Thompson, D'Arcy W. (1917). On growth and form.

Figure 1: Making a confluent layer of synthetic cell mimics. (a) Superimposed
snapshots showing a nearly circular path traced by polar chiral active
ellipsoids under vertical agitation. Clockwise $(+)$ and counterclockwise
$(-)$ moving ellipsoids are represented in blue and red colours, respectively.
(b) Polar chiral active ellipsoids align with their propulsion axis parallel
to the wall with the direction of motion along the wall determined by the
ellipsoid's handedness. (c) Granular cell enclosing $N=20$ polar chiral active
ellipsoids. Panels from left to right show cells with the chirality of the
cell interior, $\chi_{\text{Cell}}$, varying from enantiopure to racemic.
While a unidirectional polarized wall current results in cell spin, counter-
propagating currents result in cell motility, the extent of which can be tuned
by systematically decreasing $|\chi_{\text{Cell}}|$. (d) The handedness of
cell spin is determined by the sign of $\chi_{\text{cell}}$. The green and
yellow membranes spin in clockwise and anticlockwise directions due to an
excess of $(+)$ ellipsoids in one and $(-)$ ellipsoids in the other. (e)
Persistence time $\tau_{p}$ (hollow circles) and the persistence length
$l_{p}$ (solid circles) versus $|\chi_{\text{cell}}|$ of isolated cells. (f)
Shape variability of isolated cells, $SD(AR_{s})$, versus
$|\chi_{\text{cell}}|$. The greater the cell activity, the greater the shape
variability. (g) Snapshot of a nearly confluent assembly ($\phi\approx 0.94$).
The net chirality of the system, $\chi_{\text{Sys}}=0$. The plate has an equal
number of clockwise (green) and counterclockwise (yellow) spinning cells.
Figure 2: Cell shape and shape variability are correlated with dynamics for a
confluent cell monolayer. (a) Snapshots of the granular cell assemblies at
$\phi=0.44$ for three representative activities. Note a weak tendency of the
cells to form clusters with increasing $\tau_{p}$. (b) Relaxation dynamics
phase diagram in the ($\phi$, $\tau_{p}$) plane. The squares represent the
$\chi_{\text{Cell}}$ and $\phi$ at which experiments were performed. The black
dashed lines are the isochrones. The colour bar indicates the value of the
relaxation time $\tau_{\alpha}$. The $\tau_{\alpha}$ values between
experimental data points were obtained by linear interpolation. (d) & (e)
$\overline{AR}$ and $SD(AR)$ for different $\tau_{p}$ values are shown for a
range of $\phi$. The window of $\phi$ values where the re-entrant behavior was
seen is represented by yellow-filled symbols. Both $\overline{AR}$ and
$SD(AR)$ are largest at the intermediate value of $\tau_{p}$ (horizontal bar),
where the system exhibits fluid-like behavior. (f) & (g) $\tau_{\alpha}$
versus $\overline{AR}$, and $\tau_{\alpha}$ versus $SD(AR)$ for different
values of $\phi$.  Figure 3: Cell shape and shape variability are mutually
constrained for granular cell-mimics. (a) $SD(AR)$ versus $\overline{AR}$ for
granular cells adhere to the same universal scaling as seen in living systems.
The black line is the theoretical prediction of ref. [37] and has the
functional form ($SD(AR)\simeq 0.71\overline{AR}-0.75$). (b) Probability
density functions (PDFs) of $AR$ for different $\tau_{p}$ values for
$\phi=0.92$. Note that $PDF$ exhibits a broad distribution with a large
positive skew for the unjammed (fluid-like) state at an intermediate value of
$\tau_{p}=20$ s. For the more jammed states, corresponding to other $\tau_{p}$
values, the $PDF$s are narrower, and with a smaller skew. (c) When scaled by
$(AR-1)/(\overline{AR}-1)$, the PDF at different $\tau_{p}$ values collapse.
The black line represents a $k-$Gamma distribution, and the data can be
equally well-fit to the distribution proposed in ref. [37].  Figure 4: Fast
cells show suppressed shape variability. (a) Cell velocity maps over the cage
breaking time, $t^{\star}$, for different $\tau_{p}$ values for $\phi=0.92$.
The color bar denotes the displacement magnitude. The red circles represent
the top 10 $\%$ most-mobile cells over $t^{\star}$. (b) $\overline{AR}$ for
different values of $\tau_{p}$ for the top 10 $\%$ most-mobile (hollow
circles) and least-mobile cells (squares) cells over $t^{\star}$. Fast cells
typically have a larger $\overline{AR}$ than the slow cells. (c) $SD(AR)$
versus $\overline{AR}$ for the fast (circles) and slow cells (squares). The
fast cells show a reduced shape variability in both experiments (filled
symbols) and Vertex Model simulations (half-filled symbols). (d) PDF of $AR$
for all cells (hollow squares), 10 $\%$ fast (red circles) and 10 $\%$ slow
cells (solid squares) from Vertex model simulations at $T=0.009$. (e) Inset:
Mimicking dynamical heterogeneities. After equilibrating the system, a cluster
of $n$ cells (shown in cyan) is allowed to evolve while the rest of the system
is frozen (gray). Main figure: $SD(AR)$ versus $\overline{AR}$ for different
$n$ and at different $T$. The black line is the universal scaling predicted in
ref. [37]. With increasing $n$, $SD(AR)$ approaches the universal scaling from
below, demonstrating that the subdued shape variability of fast cells is due
to their confinement by their slower surrounding. Also, at low $T$, where the
mobility of the system as a whole is smaller, the deviation from the line is
also small. At high $T$, on the other hand, the confinement effect is more
pronounced.


11institutetext:  Departamento de Ciencias Exactas y Tecnología, Centro
Universitario de los Lagos, Universidad de Guadalajara, Enrique Díaz de León
1144, Colonia Paseos de la Montaña, Lagos de Moreno, Jalisco, Mexico
Lobachevsky State University of Nizhny Novgorod, Nizhny Novgorod, Russia
Centro de Tecnología Biomédica, Universidad Politécnica de Madrid, Campus de
Montegancedo, Pozuelo de Alarcón, 28223 Madrid, Spains
Innopolis University, Universitetskaya Str. 1, 420500 Innopolis, Republic of
Tatarstan, Russia

Synchronization: coupled oscillators Nonlinear dynamics and chaos

# Dynamics of a ring of three unidirectionally coupled Duffing oscillators
with time-dependent damping

J. J. Barba-Franco 11 A. Gallegos 11 R. Jaimes-Reátegui 22 S. A. Gerasimova 22
A. N. Pisarchik 22334411223344

###### Abstract

We study dynamics of a ring of three unidirectionally coupled double-well
Duffing oscillators for three different values of the damping coefficient:
fixed dumping, proportional to time, and inversely proportional to time. The
dynamics in all cases is analyzed through time series, Fourier and Hilbert
transforms, Poincaré sections, as well as bifurcation diagrams and Lyapunov
exponents with respect to the coupling strength. In the first case, we observe
a well-known route from a stable steady state to hyperchaos through Hopf
bifurcation and a series of torus bifurcations, as the coupling strength is
increased. In the second case, the system is highly dissipative and converges
into one of stable equilibria. Finally, in the third case, transient toroidal
hyperchaos takes place.

###### pacs:

05.45.Xt

###### pacs:

05.45.-a

##  1 Introduction

Duffing equations describe dynamics of a systems with cubic nonlinearity,
which can have either single-well or double-well potential. As distinct from a
harmonic oscillator described by a linear second-order differential equation,
the Duffing oscillator in its original form has essentially only one extra
nonlinear stiffness term. Despite its enigmatic simplicity, the Duffing
oscillator was successfully used to model various physical processes such as
stiffening strings, beam buckling, superconducting Josephson parametric
amplifiers, ionization waves in plasmas, as well as biomedical processes (see
[1] and references therein). The Duffing equations are easily implemented in
electronic circuits. Several physical realizations of the Duffing equation in
laboratory-based experimental models are described in the Virgin’s book [2]
and other review papers (see, e.g., [3, 4] and references therein). In
particular, a double-well Duffing oscillator was widely used to evaluate a
large variety of nonlinear systems including slender aerostructures that may
buckle under loads [3], microelectromechanical switches [5], vibration-based
energy harvesters [6, 7], electrical circuits [8], and optical systems [9].

Coupled Duffing oscillators have attracted special attention due to their
intriguing synchronization behavior [10], including two-state intermittency
[11], transition to hyperchaos [12], intermittent lag synchronization [13],
attractor annihilation in stochastic resonance [14], etc. The coupled Duffing
oscillators were extensively studied with respect to their parameters, such as
the coupling strength, nonlinearity stiffness term, external force or
modulation of parameters accessible to the system have been developed. In
particular, a ring of unidirectionally coupled Duffing oscillators exhibit the
most interesting dynamics, such as a transition from periodic to chaotic and
hyperchaotic behavior and so-called rotating wave [15, 16, 17]. The stability
of this system was estimated using Lyapunov exponents [18, 19]. Nevertheless,
little attention was paid to a study of the effect of the damping term,
although ring-coupled overdamped Duffing oscillators were investigated in the
presence a delay in coupling [20], multistability [21, 22], and even proposed
for spectrum-sensing technology application [23].

It is well known when a motion takes place in the environment, the system has
dissipation which is modeled by a damping term related to the velocity [24].
In this regard, several papers were devoted to a study of Duffing oscillators
with linear and nonlinear damping terms (see [4] and references therein). The
knowledge of the potential can be useful not only for conservative quantum
systems, but also for understanding dynamics of dissipative systems. Although
dissipative systems cannot be described by a proper potential [25], however,
in some cases the potential can still be found. For example, in the case of a
linear time-dependent damping term, the system can be viewed as an undamped
oscillator but with a variable mass and therefore the corresponding potential
can be obtained [26, 27]. In addition, analytical studies of transitions which
occur between three possible dynamical states (cluster synchronization,
complete synchronization, and instability) were performed in a ring of $N$
diffusely coupled Duffing oscillator [28, 29].

In this Letter, we study the dynamics of three double-well Duffing oscillators
coupled in a cyclic ring. As distinct for previous studies of the same
coupling configuration, in our work we focus on the effects of the damping
term. We consider three different cases of the damping coefficient: fixed
damping, damping proportional to time (overdamping case), and damping
inversely proportional to time (quasi-conservative case). To characterize the
system dynamics we use time series, power spectra, Poincaré sections,
bifurcation diagrams, and Lyapunov exponents. The third case is the most
interesting, because we observe, for the first time to the best of our
knowledge, _transient toroidal hyperchaos_.

The rest of the Letter is organized as follows. First, we describe the model
of three unidirectionally ring-coupled damped Duffing oscillators. Then, we
consider the case of the fixed damping term and demonstrate the existence of
the rotating wave in the ring. After that, we study the case when the damping
coefficient if proportional to time and demonstrate the overdamped dynamics
leading to a stable fixed point. Finally, we analyze the system dynamics when
the damping term is inversely proportional to time and demonstrate transient
toroidal hiperchaos. At the end of the Latter, we summarize the main results.

##  2 Model

Let us consider the simplest form of the undamped Duffing oscillator without a
driving force, which can be given as

$\ddot{x}+\omega_{0}^{2}x+\delta x^{3}=0,$ (1)

where the point means time derivative and $\omega_{0}^{2}$ and $\delta$ are
non null real constants. It is easy to find the corresponding potential
associated with the motion eq. (1), which has the following form

$V(x)=\frac{1}{2}\omega_{0}^{2}x^{2}+\frac{1}{4}\delta x^{4}.$ (2)

The shape of the potential function eq. (2) depends on the values of the
parameters $\omega_{0}^{2}$ and $\delta$. Basically, there are four different
cases:

  * •

If $\omega_{0}^{2}>0$ and $\delta<0$ then the potential has a double-hump well
with a local minimum at $x=0$ and two maxima at
$\pm\sqrt{\omega_{0}^{2}/\left|\delta\right|}$.

  * •

If $\omega_{0}^{2}>0$ and $\delta>0$ then the potential has a single well with
a local minimum at $x=0$.

  * •

If $\omega_{0}^{2}<0$ and $\delta<0$ then the potential has a single hump with
a local maximum at $x=0$.

  * •

If $\omega_{0}^{2}<0$ and $\delta>0$ then the potential has a double well with
two minima at $\pm\sqrt{\left|\omega_{0}^{2}\right|/\delta}$ and a local
maximum at $x=0$.

Here, we are interested in the last case, i.e. bistability. Since a damping
term is included in eq. (1), the expected motion is ruled by the potential eq.
(2) but with dissipation until the trajectory is attracted to a stable fixed
point. The question is if it is possible to define a proper potential for this
kind of damped systems.

It is well known that in general the damped systems do not always have a
defined potential [25]. However, it is possible to find potentials for
dissipative systems if the damping coefficient can be written as logarithmic
derivative of certain function (see, e.g., [27, 26]). More explicitly, if we
have a motion equation of the type

$\ddot{x}+\alpha(t)\dot{x}+\beta(x,t)=0,$ (3)

it can be described by a Langrangian of the form

$L=\frac{1}{2}m(t)\dot{x}^{2}-V(x,t),$ (4)

where

$\alpha(t)=\dot{m}/m,\qquad V(x,t)=m(t)\int^{x}\beta(z,t)dz,$ (5)

and $\alpha(t)$ is a time-dependent damping coefficient.

In other words, eq. (3) can be viewed as an undamped motion with variable mass
$m(t)$. Therefore, including $\alpha(t)$ to eq. (1) we get

$\ddot{x}+\alpha(t)\dot{x}+\omega_{0}^{2}x+\delta x^{3}=0,$ (6)

and the corresponding potential is given as

$V(x,t)=m(t)\left[\frac{1}{2}\omega_{0}^{2}x^{2}+\frac{1}{4}\delta
x^{4}\right],$ (7)

where $\alpha(t)=\dot{m}/m$. In other words, we are talking about the same
potential eq. (2) with similar fixed points but with a scaling factor $m(t)$.
However, we will show in the next sections that different types of time
dependence of the damping coefficient can change the dynamics of the ring-
coupled Duffing oscillatory system.

The ring of three unidirectionally coupled Duffing oscillators is described as
follows [27, 22]

$\displaystyle\ddot{x}_{1}+\alpha(t)\dot{x}_{1}+\omega_{0}^{2}x_{1}+\delta
x_{1}^{3}+\sigma(x_{1}-x_{3})$ $\displaystyle=0,$ (8)
$\displaystyle\ddot{x}_{2}+\alpha(t)\dot{x}_{2}+\omega_{0}^{2}x_{2}+\delta
x_{2}^{3}+\sigma(x_{2}-x_{1})$ $\displaystyle=0,$
$\displaystyle\ddot{x}_{3}+\alpha(t)\dot{x}_{3}+\omega_{0}^{2}x_{3}+\delta
x_{3}^{3}+\sigma(x_{3}-x_{2})$ $\displaystyle=0,$

where $\sigma$ is the coupling strength.

In the rest of the paper, we set $\omega_{0}^{2}=-0.25$, $\delta=0.5$ which
are parameters related to a bistable Duffing oscillator that has three fixed
points associated with the potentials eq. (2) or (7) given as

$\displaystyle x_{u,d}$
$\displaystyle=\pm\sqrt{\left|\omega_{0}^{2}\right|/\delta}=\pm\sqrt{1/2},$
(9) $\displaystyle x_{0}$ $\displaystyle=0,$

where fixed points $x_{u,d}$ are stable (subindices $u$ and $d$ correspond to
positive and values) and $x_{0}$ is unstable.

By converting the second-order eq. (8) into first order equations using the
change of the variable $\dot{x}=y$, the dynamics of the $j$th oscillator in
the ring can be described by the following pair of first order dimensionless
ordinary differential equations:

$\displaystyle\dot{x}_{j}$ $\displaystyle=y_{j},$ (10)
$\displaystyle\dot{y}_{j}$
$\displaystyle=-\alpha(t)y_{j}+\omega_{0}^{2}x_{j}-\delta
x_{j}^{3}+\sigma(x_{j-1}-x_{j}),$

where $\sigma$ is the coupling coefficient for each oscillator $j=1,2,3$.

In the next sections, we will consider three cases of the time-dependent
damping coefficient: $\alpha=0.4$ (constant coefficient), $\alpha=t/4$
(linearly increasing in time), and $\alpha=1/t$ (linearly decreasing in time).

Figure 1: (a) Bifurcation diagram of the local max of $x_{1}$ and (b) four
largest Lyapunov exponents $\lambda$ versus coupling strength $\sigma$ for
$\omega_{0}^{2}=-0.25$, $\delta=0.5$, and $\alpha=0.4$.

##  3 Dynamics of the system with fixed damping

First, we consider the case of the fixed damping coefficient $\alpha(t)=0.4$.
Due to the symmetrical coupling of the three Duffing oscillators eq. (8), the
analysis is made with the bifurcation diagram of the local maxima of the
amplitude of one of the oscillators (e.g., $x_{1}$) and four largest Lyapunov
exponents $\lambda$ with respect to the coupling strength $\sigma$. These
diagrams are shown in fig. 1. The observed bifurcation scenario from the
equilibrium point to chaos and hyperchaos via subsequent Hopf bifurcations is
in a good agreement with the Landau-Hopf transition to turbulence [30, 31] and
the Newhouse, Ruelle, and Takens theorem [32] stated that just after the
successive Hopf bifurcation a torus decays into a strange chaotic attractor.

(a) (b) (c) (d) (e)

Figure 2: (i) Time series, (ii) Poincaré sections, and (iii) power spectra for
different values of the coupling strength (a) $\sigma=0.354$, (b)
$\sigma=0.5$, (c) $\sigma=1.078$, (d) $\sigma=1.98$, and (d) $\sigma=3.25$.
$\omega_{0}^{2}=-0.25$, $\delta=0.5$, $\alpha=0.4$.

In fig. 1 we observe the well-known scenario from a stable steady state to
hyperchaos when the coupling strength $\sigma$ is increased. The details of
dynamical regimes on this route are presented in fig. 2 with the time series,
Poincaré sections, and power spectra. Starting from $\sigma=0$, first, a
stable fixed point transforms into a limit cycle in the Hopf bifurcation when
the coupling reaches $\sigma_{1}\approx 0.35$, where the largest Lyapunov
exponent reaches zero (black line in Fig. 2(b)). This regime is illustrated in
fig. 2(a) and maintained within a relatively small region of
$0.35<\sigma<0.5$. Then, at $\sigma_{2}=0.5$ the limit cycle transforms into a
quasiperiodic regime (2D torus) when the second frequency appears in the power
spectrum (Fig. 2(b)). This quasiperiodic regime appears when the second
largest Lyapunov exponent becomes zero and the third largest Lyapunov exponent
approaches zero. When $\sigma$ is further increased, a 3D torus arises at
$\sigma_{3}=1.0$ when the third largest Lyapunov exponent approaches zero.
This regime is observed for $1.0<\sigma<1.75$ and characterized by a large
number of frequencies in the power spectrum (fig. 2(c)). At $\sigma_{4}=1.75$
the system becomes chaotic because the largest Lyapunov exponent becomes
positive (fig. 2(b)) and the power spectrum is wide (fig. 2(d)). A further
increase in the coupling strength leads to hiperchaos at $\sigma_{5}$ when two
largest Lyapunov exponents become positive; this regime is illustrated in fig.
2(e)).

Another interesting dynamical feature of the ring-coupled oscillators is the
existence of the rotation wave in the quasiperiodic and chaotic regimes which
consists in the constant phase difference between the envelope (second
frequency) of quasiperiodic and chaotic oscillations of each oscillator. This
fascinating synchronization state was first discovered in ring-coupled Chua
oscillators [33, 34] and then found in ring-coupled Lorenz [35, 36, 37] and
Duffing oscillators [15, 16, 17]. In fig. 3(a) we plot the maximum spectral
component $S_{0}$ at the dominant frequency $\Omega_{0}$ (main oscillation
frequency) and the rotating wave spectral power $S_{W}$ at the wave frequency
$\Omega_{W}$ (envelope freqeuncy) as functions of the coupling strength. The
latter is obtained from the power spectrum of the Hilbert transform and shown
in fig. 3(b). One can see from fig. 4(a) that both spectral components
increases as $\sigma$ is increased up to $\sigma=3$ and then the powers become
decrease. This decrease happens because in this hyperchaotic regime the
spectral energy is distributed over a wide frequency range, as seen in fig.
2(e)(iii)). The wave frequency $S_{W}$ is almost independent of $\sigma$ in
the chaotic and hyperchaotic regions, i.e. for $\sigma>1.75$, although for
very large couplings ($\sigma>3$) it is not well determined.

Figure 3: (a) Maximum ($S_{0}$) and rotating wave ($S_{W}$) spectral
components and (b) their frequencies ($\Omega_{0}$ and $\Omega_{W}$) of
$x_{1}$ versus coupling strength $\sigma$. $\omega_{0}^{2}=-0.25$,
$\delta=0.5$, $\alpha=0.4$.

The underlying mechanism of the rotating wave stability is explained as the
effect of an additional rotational degree of freedom and the symmetric
structure of the coupling. Previous results associate the rotational wave with
the space-time symmetry property of an invariant ring of identical oscillators
under cyclic group [38, 39]. Our results obtained with double-well Duffing
oscillators confirm the previous findings.

##  4 Damping coefficient proportional to time ($\alpha(t)=t/4$)

Let us now consider the second case, when the damping coefficient is
proportional to time, specifically, $\alpha(t)=t/4$. In this case, the system
is extremely dissipative and therefore the rotating wave does not appear.
Instead, starting from different initial conditions the system tends to one of
two stable equilibrium points ($x_{1}=\pm 0.71,y_{1}=0$) of the potential
$V(x_{1},t)$ (see eqs. (7) and (9)), as illustrated with the time series and
phase portraits in figs. 4(a) and (b), respectively.

Figure 4: (a) Time series and (b) phase portraits on the ($x_{1}$, $y_{1}$)
plane for $\omega_{0}^{2}=-0.25$, $\delta=0.5$, and $\alpha(t)=t/4$.

Figures 5(a) and (b) show, respectively, the bifurcation diagram of oscillator
$x_{1}$ and three largest Lyapunov exponent as functions of the coupling
strength $\sigma$. One can see in fig. 5(a) that the system switches to
another coexisting fixed point although the initial conditions are fixed. This
happens because for these values of the coupling parameters the initial
conditions hit to the basin of attraction of another equilibrium. As seen from
fig. 5(b), all Lyapunov exponents are negative for any value of the coupling
strength, and the system becomes more stable as $\sigma$ is increased. So, the
high damping coefficient does not allow the rotating wave.

Figure 5: (a) Bifurcation diagram of $x_{1}$ and (b) three largest Lyapunov
exponents versus coupling strength for $\omega_{0}^{2}=-0.25$, $\delta=0.5$,
and $\alpha(t)=t/4$.

##  5 Damping coefficient inversely proportional to time ($\alpha(t)=1/t$)

Finally, we consider the case when the damping coefficient is inversely
proportional to time, specifically, $\alpha(t)=1/t$, i.e., the damping tends
to zero over time. The time series and corresponding instantaneous frequency
of $x_{1}$ calculated using the Hilbert transform are shown in figs. 6(a) and
(b), respectively. One can see that such damping function results in the
transient behavior followed by infinity. Moreover, in fig. 6(b) one can
observe continuous changes of the instantaneous frequency during all time
evolution. This evidences the existence of a permanent hyperchaotic
transitory.

Figure 6: (a) Time series of $x_{1}$ and (b) corresponding instantaneous
frequency, representing transient hyperchaos. $\omega_{0}^{2}=-0.25$,
$\delta=0.5$, $\sigma=3.25$, and $\alpha(t)=1/t$.

The transient hyperchaos manifests itself as follows. The dissipative system
eq. (8) turns into conservative over time, causing hyperchaotic transient
behavior, because a non-attracting chaotic set coexists with a chaotic
attractor, i.e., there are two distinct forms of chaotic behavior. The
trajectory proceeding from randomly chosen initial conditions looks chaotic
for a sufficiently long period of time, during which it visits various chaotic
and quasiperiodic states rather abruptly, as seen from the time series in fig.
7 corresponding to different time periods.

Figure 7: Time series in different periods of time representing (a) chaotic,
(b) quasiperiodic of $x_{1}$ and (b) corresponding instantaneous frequency,
representing transient hyperchaos. $\omega_{0}^{2}=-0.25$, $\delta=0.5$,
$\sigma=3.25$, and $\alpha(t)=1/t$.

Transient chaotic behavior appears in a conservative system where the phase-
space volume is constant under time evolution [40]. In the Poincaré section on
($x_{1}$, $y_{1}$) plane (fig. 8) one can distinguish two regions. The region
of many scattered points represents non-attracting chaotic set called chaotic
saddle [41] at the beginning of the transients (large amplitude oscillations
at $t<2500$), whereas the central bagel of high density points corresponds to
the resting part of the transients and represents toroidal hyperchaos.

Figure 8: Poincaré section on ($x_{1}$, $y_{1}$) plane representing transient
toroidal hyperchaos for $\omega_{0}^{2}=-0.25$, $\delta=0.5$, $\sigma=3.25$,
and $\alpha(t)=1/t$.

Figure 9(a) and (b) show the bifurcation diagram of local maxima of $x_{1}$
and four largest Lyapunov exponents versus $\sigma$, respectively. One can see
that when the coupling strength is increased, the system exhibits a fast
transition to hyperchaos from a fixed point via a crisis bifurcation at a very
small coupled strength $\sigma=0.056$. For larger $\sigma$, two largest
Luyapunov exponents are always positive and the third exponent is zero. This
means that the system is hyperchaotic.

Figure 9: (a) Bifurcation diagram of $x_{1}$ and (b) largest Lyapunov
exponents versus coupling strength for $\omega_{0}^{2}=-0.25$, $\delta=0.5$,
$\sigma=3.25$, and $\alpha(t)=1/t$.

##  6 Conclusion

In this Letter, we have studied the dynamics of a system of three double-well
Duffing oscillators unidirectionally coupled in a ring. We have considered
three cases depending on the damping coefficient: constant damping, damping
proportional to time ($\alpha(t)=t/4$), and damping inversely proportional to
time ($\alpha(t)=1/t$). The system dynamics have been analyzed using time
series, Fourier and Hilbert transforms, Poincaré sections, and Lyapunov
exponents. In the first case, we observed the route from a steady state to
hyperchaos through a series of torus bifurcations, as the coupling strength is
increased, as well as the existence of a rotating wave with a fixed phase
between neighboring oscillators, which persists for quasiperiodic, chaotic and
hyperchaotic regimes. These results are in a good agreement with previous
studies of single-well Duffing oscillators in the same coupling configuration.
Such similarity is explained by the fact that in the quasiperiodic and chaotic
regimes, the system becomes monostable.

In the second case, the system becomes highly dissipative and therefore does
not generate a rotating wave. As a result, the system, is evolved to one of
the stable steady states which depends on the initial conditions. During the
transient behavior, the system can occasionally switch between two coexisting
states.

Finally, in the third case, the system becomes conservative when the time $t$
involves to large values, and then goes to infinity. The transient toroidal
hyperchaotic behavior was observed. During transients, the system can visit
different unstable periodic orbits.

Although this work was devoted to the study of the simple network motif of
only three coupled oscillator, we suppose that similar dynamics may occur in
larger oscillatory networks (see, e.g., [42]), especially in a larger ring of
unidirectionally coupled Duffing oscillators. This is a promising topic for
future research.

###### Acknowledgements.

J. J. B. F. thanks the National Council for Science and Technology of Mexico
(CONACYT) for the financial support granted through the scholarship number
924190. S. A. G. and A. N. P. acknowledge support from the Lobachevsky
University Competitiveness Program in the frame of the 5-100 Russian Academic
Excellence Project.

## References

  * [1] M. Lakshmanan K. Murali Chaos in Nonlinear Oscillators: Controlling and Synchronization World Scientific 1996\.
  * [2] L. N. Virgin Introduction to Experimental Nonlinear Dynamics: a Case Study in Mechanical Vibration Cambridge Univ. Press 2000\.
  * [3] L. N. Virgin Vibration of Axially-Loaded Structures Cambridge Univ. Press 2007\.
  * [4] I. Kovacic M. J. Brennan The Duffig Equation: Nonlinear Oscillators and Their Behaviour John Wiley & Sons 2011\.
  * [5] J. Qiu J. H. Lang A. H. Slocum J. Microelectromech. Syst.132004137.
  * [6] T. J. Kazmierski S. Beeby Energy Harvesting Systems Springer 2014\.
  * [7] R. L. Harne K. Wang Smart Mater. Struct.222013023001.
  * [8] G. Debnath T. Zhou F. Moss Phys. Rev. A3919894323.
  * [9] M. Dykman A. Velikovich G. Golubev D. Luchinskii S. Tsuprikovs J. Exper. Theor. Phys. Lett.531991193.
  * [10] S. Boccaletti A. N. Pisarchik C. I. Del Genio A. Amann Synchronization: From Coupled Systems to Complex Networks Cambridge Univ. Press, 2014\.
  * [11] R. Jaimes Reategui A. N. Pisarchik Phys. Rev. E692004067203.
  * [12] T. Kapitaniak Phys. Rev. E471993R2975.
  * [13] A. N. Pisarchik R. Jaimes-Reategui Phys. Lett. A3382005141.
  * [14] A. N. Pisarchik U. Feudel Phys. Rep.5402014167.
  * [15] P. Perlikowski S. Yanchuk M. Wolfrum A. Stefanski P. Mosiolek T. Kapitaniak Chaos202010013111.
  * [16] L. Borkowski P. Perlikowski T. Kapitaniak A. Stefanski Phys. Rev. E202015062906.
  * [17] L. Borkowski A. Stefanski Eur. Phys. J. Spec. Top.22920202249.
  * [18] A. Dabrowski Nonlin. Dyn.672012283.
  * [19] M. Balcerzak D. Pikunov A. Dabrowski Nonlin. Dyn.9420183053.
  * [20] M. V. Tchakui P. Woafo Chaos262016113108.
  * [21] C. Meena P. D. Rungta S. Sinha Eur. Phys. J. B932020210.
  * [22] R. Jaimes-Reátegui J. M. Castillo-Cruz J. H. García-López G. Huerta-Cuellar L. A. Gallegos-Infante A. N. Pisarchik Cybern. Phys.9202031.
  * [23] Z.-L. Tang S.-M. Li L.-J. Yu Sensors162016844.
  * [24] L. D. Landau, E. M. Lifshiftz Mechanics Pergamon Press 1969\.
  * [25] R. Graham T. Tél Phys. Rev. Lett.5219849.
  * [26] J. L. Cieśliński T. Nikiciuk J. Phys. A432010175205.
  * [27] J. J. Barba-Franco P. B. Espinoza A. Gallegos R. Jaimes-Reategui J. E. Macias-Diaz Rev. Mex. Fis.662020440.
  * [28] Y. C. Kouomou P. Woafoz Phys. Rev. E662003046205.
  * [29] V. T. Yolong P. Woafo Phys. Scr.742006591.
  * [30] L. D. Landau Dokl. Akad. Nauk USSR441944311.
  * [31] E. Hopf Commun. Pure Appl. Math.11948303.
  * [32] S. Newhouse D. Ruelle F. Takens Commun. Math. Phys.64197835.
  * [33] M. A. Matías J. Göémez V. Pérez-Muñuzuri I. P. Mariño M. N. Lorenzo V. Pérez-Villar Europhys. Lett.371997379.
  * [34] I. P. Mariño, V. Pérez-Muñuzuri, V. Pérez-Villar, E. Samchez, M. A. Matías Physica D1281999224.
  * [35] M. A. Matías J. Göémez Phys. Rev. Lett.8119984124.
  * [36] X. L. Deng H. B. Huang Phys. Rev. E652002055202.
  * [37] E. Sánchez D. Pazo M. A. Matías Chaos162006033122.
  * [38] J. J. Collins I. Stewart Biol. Cybern.71199495.
  * [39] D. Pazó E. Sánchez M. A. Matías Intern. J. Bifurc. Chaos1120012683.
  * [40] Y.-C. Lai T. Tél Transient Chaos: Complex Dynamics on Infinite Time Scales Springer Science & Business Media 2011\.
  * [41] S. Sabarathinam K. Thamilmaran Chaos Solitons Fractals732015129.
  * [42] R. Milo S. Shen-Orr S. Itzkovitz N. Kashtan D. Chklovskii U. Alon Science2982002824.