# Self-sustained oscillations of active viscoelastic matter

Emmanuel L. C. VI M. Plan  Institute of Theoretical and Applied Research, Duy
Tan University, Ha Noi 100 000, Viet Nam  Faculty of Natural Science, Duy Tan
University, Da Nang 550 000, Viet Nam  Huong Le Thi  Department of Mathematics
and Informatics, Thang Long University, Ha Noi 100 000, Viet Nam  Julia M.
Yeomans  The Rudolf Peierls Centre for Theoretical Physics, Department of
Physics, University of Oxford,
Clarendon Laboratory, Oxford OX1 3PU, United Kingdom  Amin Doostmohammadi  The
Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, 2100
Copenhagen, Denmark

###### Abstract

Models of active nematics in biological systems normally require complexity
arising from the hydrodynamics involved at the microscopic level as well as
the viscoelastic nature of the system. Here we show that a minimal, space-
independent, model based on the temporal alignment of active and polymeric
particles provides an avenue to predict and study their coupled dynamics
within the framework of dynamical systems. In particular, we examine, using
analytical and numerical methods, how such a simple model can display self-
sustained oscillations in an activity-driven viscoelastic shear flow.

##  I Introduction

Active nematics are a class of living or synthetic materials that possess some
form of orientational order and that are continuously driven by their
constituent living or artificially active elements [1, 2]. The orientational
order could arise from the elongated shape of the particles, such as in
filamentous cytoskeletal elements [3, 4, 5], rod-like bacteria [6, 7, 8] or
spindle-shaped fibroblasts [9] and neural progenitor stem cells [10], or could
be an emergent feature of deformable particles such as in monolayers of
epithelial cells [11, 12]. Whether it is moving bacteria or eukaryotic cells,
or subcellular filaments that are put in motion by motor proteins, the common
feature of all these systems is the input of energy at the particle level,
that governs system-wide patterns of motion and self-organisation in the
active material. While active nematic models have been indispensable in
characterising various systems of active matter, their inherent complexity,
that requires coupling between density, velocity, and the orientational field,
still limits theoretical advances and often requires extensive numerical
computations in order to obtain generic and fundamental rules of dynamical
self-organisation in active matter.

To add to this complexity, most realisations of active matter do not simply
exist in homogeneous and isotropic environments. Rather, the microenvironment
surrounding living materials is often characterised by its own complex
structure and responses that can feedback to the active matter dynamics.
Prominent among these is viscoelasticity that allows living materials to
constantly deform their surroundings and then to exploit the deformation and
relaxation dynamics of their viscoelastic environment as an external spatio-
temporal input for their own self-organisation. Examples include bacterial
biofilm formation in extracellular matrices [13, 14, 15], epithelial and
fibroblast cells self-propelling within interconnected networks of collagen
[16, 17, 18], and subcellular filaments self-organising within the crowded
environment of the cytoplasm [19, 20]. Recent theoretical and experimental
studies have revealed the emergence of various forms of collective behavior
for active particles suspended within - or in contact with - a viscoelastic
environment. In particular, numerical models of active viscoelastic nematics
[21] and active nematics in contact with viscoelastic surroundings [22] have
demonstrated the impact of polymeric stresses and viscoelastic relaxation
times on active nematic patterns in confined geometries and in models of cell
self-propulsion and cell division. Oscillations have been shown using a full
space-dependent hydrodynamic models of active matter which includes
viscoelasticity [21]. Moreover, combined experiments and theory have recently
shown the emergence of directional switching in bacterial collectives
suspended within a viscoelastic medium and confined in circular geometries
[23]. Consistent with the experiments, a simplified model that maps into a
Fitzhugh-Nagumo oscillator was deduced [23].

Here we show that an alternative model of active viscoelastic matter relying
only on the orientation of both active and viscoelastic components along with
their respective relaxation timescales can also produce oscillatory dynamics.
This model was first introduced in a recent paper [24] for active matter in
contact with viscoelastic surroundings to reproduce numerical observations of
the reversals of the flow field in a nematic driven by pulsatile activity. We
explore the generic dynamics of this minimal, space-independent model within
the framework of dynamical systems, and prove analytically that sustained
oscillatory behaviour can arise spontaneously as an effect of the viscoelastic
feedback to the active forcing. We provide a phase diagram for the emergence
of the oscillations in terms of the dominant time scales of the system, and
differentiate their form from sinusoidal oscillations. The equations we
consider are not algebraic and cannot be naively truncated nor mapped to well-
known oscillators. Also, we remark that the physics described in this work is
entirely different from that seen in the oscillations in purely active models
[25, 26], because of the presence here of a polymeric component.

##  II Model

Consider an active nematogen and a polymer, oriented at angles $\theta_{n}$
and $\theta_{p}$ respectively with respect to the $x$-axis. Then, modelling
both the nematogen and the polymer as rigid rods rotating in 2D in the shear
$\dot{\gamma}$ created by the active stresses and the polymer relaxation, we
consider the simplified, space-independent, equations of motion (first
proposed in [24])

$\displaystyle\dot{\theta}_{n}$ $\displaystyle=$
$\displaystyle\dot{\gamma}\lambda_{n}\cos 2\theta_{n}-(\theta_{n}/t_{n}),$ (1)
$\displaystyle\dot{\theta}_{p}$ $\displaystyle=$
$\displaystyle\dot{\gamma}\lambda_{p}\cos 2\theta_{p}-(\theta_{p}/t_{p}),$ (2)
$\displaystyle\dot{\gamma}$ $\displaystyle=$
$\displaystyle\frac{\zeta}{4\nu}\sin 2\theta_{n}-\frac{\nu_{p}}{4\nu
t_{p}}\sin 2\theta_{p}.$ (3)

The first terms on the right-hand side of Eqs. (1) and (2) describe the
nematic and polymeric response to the shear, $\dot{\gamma}$, in the limit that
vorticity can be neglected. $\lambda_{n},\lambda_{p}>0$ are phenomenological
parameters related to the aligning/tumbling properties, see for example [27].
The second terms model exponential relaxation to equilibrium, with
$t_{n},t_{p}$ the relaxation timescales of the nematic and polymer. The shear
rate, Eqn. (3), is the off-diagonal term in the stress balance $\Pi_{\rm
visc}=\Pi_{\rm active}+\Pi_{\rm polymer}$, where $\nu$ is the fluid viscosity
and $\nu_{p}$ is the polymeric contribution to viscosity [24]. The dynamics of
the nematogen and polymer are coupled through the shear which describes the
balance between the extensile ($\zeta>0$) stress resulting from the activity
of the active nematic and the contractile stress created by the polymer
relaxation.

These equations represent a simplification motivated by the full hydrodynamic
model investigated numerically in [24] which showed that activity pulses in a
system interacting with a viscoelastic environment can spontaneously generate
flow reversals driven by the stress built up in the surrounding, passive,
medium. In particular, replacing the space-dependent flow by a uniform shear
results in only two dynamical variables, which allows us to fully investigate
the phase space and to prove the existence of oscillations.

##  III Results and Discussion

We discuss the results in two steps: first we will fix the value of the
polymer viscosity $\nu_{p}$ and then we keep the elastic modulus, i.e., the
ratio of the polymer viscosity to the polymer relaxation time
$A_{p}=\nu_{p}/t_{p}$ constant.

###  III.1 Fixed polymer viscosity $\nu_{p}$

We first consider the case when polymer contribution to viscosity $\nu_{p}$ is
fixed and the elastic modulus $A_{p}$ varies inversely with $t_{p}$. Note that
in the absence of activity ($\zeta=0$) the system relaxes to the equilibrium
and trivial fixed point $\theta_{n}=\theta_{p}=0$. Moreover, for finite
activity, in the limit of infinite relaxation times, only the active force
dictates the shear and thus the steady-state solution is
$\theta_{n}=\pm\pi/4,\theta_{p}=\pm\pi/4$. Henceforth, without losing
generality we restrict our analysis of the system to the box
$\Omega=[-\pi/4,\pi/4]^{2}$. This two-dimensional system
$\bm{\theta}^{T}=(\theta_{n},\theta_{p})$ can then be written as
$\dot{\bm{\theta}}=\bm{f}(\bm{\theta})$, with
$\bm{f}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$. By introducing the time
scale of active stress generation $t_{a}=\nu/|\zeta|$, Eqs. (1),(2) can be
expressed in terms of dimensionless numbers $T_{p}=t_{p}/t_{a}$ and
$T_{n}=t_{n}/t_{a}$:

$\displaystyle\dot{\theta}_{n}$ $\displaystyle=f_{n}(\theta_{n},\theta_{p})$
$\displaystyle=\frac{1}{t_{a}}\left(\frac{\lambda_{n}}{4}\cos 2\theta_{n}\sin
2\theta_{n}-\frac{\lambda_{n}\nu_{p}}{4\nu T_{p}}\cos 2\theta_{n}\sin
2\theta_{p}-\frac{\theta_{n}}{T_{n}}\right),\hfill$ (4a)
$\displaystyle\dot{\theta}_{p}$ $\displaystyle=f_{p}(\theta_{n},\theta_{p})$
$\displaystyle=\frac{1}{t_{a}}\left(\frac{\lambda_{p}}{4}\cos 2\theta_{p}\sin
2\theta_{n}-\frac{\lambda_{p}\nu_{p}}{4\nu T_{p}}\cos 2\theta_{p}\sin
2\theta_{p}-\frac{\theta_{p}}{T_{p}}\right).$ (4b)

From Eqs. (4), it is easy to see that the active timescale modulates the speed
of the evolution of the system towards any asymptotic limit, with larger
$t_{a}$ meaning slower dynamics.

To examine long-term dynamics, a linear stability analysis is performed by
examining at any fixed point $(\theta_{n}^{*},\theta_{p}^{*})$ the eigenvalues
$e_{i}=(\operatorname{Tr}J\pm\sqrt{(\operatorname{Tr}J)^{2}-4\operatorname{Det}J})/2,$
where $J_{ij}=\partial f_{i}/\partial\theta_{j}$ is the Jacobian. The origin
$(0,0)$ is a trivial fixed point for all values of the parameters, and a short
calculation shows that a necessary condition for stability is
$\operatorname{Tr}J=\frac{\lambda_{n}}{2}-\frac{1}{T_{n}}-\frac{\lambda_{p}\nu_{p}}{2\nu
T_{p}}-\frac{1}{T_{p}}<0$. For simplicity, we assume
$\lambda_{n}=\lambda_{p}=\nu=\nu_{p}=1$ for the rest of this work. The
condition for stability then simplifies to $T_{p}<3T_{n}/(T_{n}-2)$.

Thus any trajectory that begins in the neighborhood of the origin either
asymptotically approaches this point if it is stable, or wanders away from it
if $\operatorname{Tr}J>0$. This is indeed the case, as is illustrated in a
typical phase diagram in Fig. 1a, where the red solid line gives the line
$T_{p}=3T_{n}/(T_{n}-2)$ and the blue circles above it denote the region where
a trajectory is not attracted to the trivial fixed point. If the trivial fixed
point is attractive (points below the red line), the angles
$\theta_{n},\theta_{p}$ asymptotically align towards 0, possibly in an
oscillatory motion. Physically speaking, small values of either $T_{n}$ or
$T_{p}$ imply rapid orientational relaxation dominating over the active
forces.

Figure 1: (a) Phase diagram in the $T_{p}$-$T_{n}$ space for
$\lambda_{n}=\lambda_{p}=\nu=1$. Crosses ($\times$) denote the region where
the only fixed point is the trivial one, whether stable or unstable. Circles
($\circ$) denote the region where the trivial fixed point is unstable. The red
shaded region is where oscillations occur. The solid red line is
$T_{p}=3T_{n}/(T_{n}-2)$ and the dashed line is a numerical estimate of the
second border of the oscillatory region. (b) Bifurcation diagram showing
$\theta_{n}(\theta_{p})$ as a function of $T_{n}$, with $T_{p}=60$. Stable
points are denoted by circles ($+$ marks) in blue and unstable points are
denoted by cross (square) marks in red. The gray-shaded region ($T_{n}\lesssim
2.11$) is where the trivial fixed point is stable, green is where oscillatory
behavior occurs, and the white region ($T_{n}\gtrsim 2.23$) is where more than
one fixed point is present. The fixed points with negative values of
$\theta_{n},\theta_{p}$ will be a reflection of the current figure with
respect to the axis.

Aside from the trivial fixed point, the mix of algebraic and transcendental
functions in $\bm{f}$ prevents us from analytically obtaining other fixed
points, unlike other relaxation oscillators in literature (e.g. the Van der
Pol equation)[28]. We thus approach the problem numerically and obtain both
the location and the stability of any other fixed points of this system, as
demonstrated in the bifurcation diagram in Fig. 1b. To obtain this diagram
trajectories were initiated within the entire space and allowed to evolve
until convergence to unique fixed points. The locations of these fixed points
can be anticipated by drawing the nullclines
$\dot{\theta}_{n}=\dot{\theta}_{p}=0$ for different values of the parameters
and dimensionless values. A similar diagram can be constructed for when the
bifurcation parameter is $T_{p}$.

There are three distinct regions in this figure. The first region, when the
dimensionless nematic relaxation time, $T_{n}$, is sufficiently small (gray
region), corresponds to when the trivial fixed point is stable. Here there are
no other fixed points, and the origin provides the global attracting point of
the system. The trivial fixed point loses its stability with increasing
$T_{n}$. By contrast, if $T_{n}$ is sufficiently large (white region), a
saddle node bifurcation occurs at the intersection of the nullclines
$\dot{\theta}_{n}=\dot{\theta}_{p}=0$ and spontaneously gives birth to a pair
of stable and unstable nodes [28]. The stable node approaches
$(\pm\pi/4,\pm\pi/4)$ as $T_{n}\rightarrow\infty,$ as expected from our
earlier arguments.

For a narrow intermediate range of the dimensionless number $T_{n}$ (the green
region) the only fixed point is that at the origin, which is now unstable. We
show that in this region the dynamical system undergoes an oscillatory motion
in the nematic and polymer alignment angles $\theta_{n},\theta_{p}$. To do
this, we will utilize the Poincaré-Bendixson Theorem [28, 29], which can be
invoked given a two-dimensional, continuously differentiable differential
equation defined on an open set encompassing a compact subset $\Omega^{*}$
that contains a trajectory confined within $\Omega^{*}$. The Poincaré-
Bendixson Theorem asserts that if $\Omega^{*}$ does not contain any fixed
points, then $\Omega^{*}$ contains a closed orbit. In our case, this closed
orbit corresponds to an oscillatory motion in the time evolution of the
angles.

Note that the $f_{i}$ as defined in Eqs. (4) are continuous in
$\mathbb{R}^{2}\supset\Omega$ and it can be easily shown (see Appendix) that
$\Omega$ is a trapping region, i.e. all trajectories at the boundary point
towards the interior of $\Omega$, confining all trajectories within $\Omega$.
By erecting a small neighborhood $\mathcal{B}$ of radius $\epsilon$ around the
repeller $(0,0)$, then
$\Omega^{*}=\Omega\setminus\mathcal{B}(\bm{0};\epsilon)$ is a compact trapping
region. As a consequence of the Poincaré-Bendixson Theorem, there is a closed
orbit in this region, which implies oscillatory motion in
$\theta_{n},\theta_{p}$, as exemplified in the time series in Fig. 2a and the
red orbit in Fig. 2b. The vector field in this figure further supports the
attractive nature of this limit cycle.

Figure 2: (a) Time series displaying oscillatory motion for
$T_{n}=2.2,T_{p}=60,t_{a}=40$. (b) Corresponding vector plot and trajectory
(red dotted line) of Eqs. (4) in the $\theta_{n}$-$\theta_{p}$ space. Note
that the vectors in the figure are rescaled to the maximum magnitude in the
plot using the Matlab visualization package Streak [30].

Note that the bifurcation diagram in Fig. 1b is simply a slice of the red
region in the full phase diagram in Fig. 1a, taken at a dimensionless polymer
relaxation time, $T_{p}=60$. Indeed it is in this red region of the phase
diagram, corresponding to only one fixed point which is unstable, where the
conditions of the Poincaré-Bendixson Theorem are satisfied. The transition
from an attractive fixed point into a limit cycle shows that the curve
$T_{p}=3T_{n}/(T_{n}-2)$ marks a supercritical Hopf bifucation.

###  III.2 Fixed elastic modulus $A_{p}$

Next we study the case with fixed elastic modulus $A_{p}=\nu_{p}/t_{p}$. In
contrast to the Eqs. (4) above, the dimensionless relaxation variable $T_{p}$
ceases to appear in the shear terms, and there is a less simple dependence on
the active timescale $t_{a}$. In this case, Eqs. (1)-(3) become

$\displaystyle\dot{\theta}_{n}$
$\displaystyle=\frac{1}{t_{a}}\left(\frac{\lambda_{n}}{4}\cos 2\theta_{n}\sin
2\theta_{n}-\frac{\lambda_{n}A_{p}t_{a}}{4\nu}\cos 2\theta_{n}\sin
2\theta_{p}-\frac{\theta_{n}}{T_{n}}\right),\hfill$ (5a)
$\displaystyle\dot{\theta}_{p}$
$\displaystyle=\frac{1}{t_{a}}\left(\frac{\lambda_{p}}{4}\cos 2\theta_{p}\sin
2\theta_{n}-\frac{\lambda_{p}A_{p}t_{a}}{4\nu}\cos 2\theta_{p}\sin
2\theta_{p}-\frac{\theta_{p}}{T_{p}}\right).$ (5b)

As in the previous case of constant polymer viscosity, the origin is a fixed
point for all values of $T_{n},T_{p}$. However, in the limit of vanishing
relaxation $T_{n},T_{p}\rightarrow\infty$, all the points in the curve
$4\dot{\gamma}=\sin 2\theta_{n}-A_{p}\sin 2\theta_{p}/\nu=0$ are fixed points
of the system, and the domain to be considered cannot a priori be limited to
$\Omega$. Nevertheless, in the finite $T_{n},T_{p}$ case, it will be possible
to restrict the analysis to $\Omega$ and to prove that oscillatory motion may
still occur.

Performing the same stability analysis at the origin (with finite
$T_{n},T_{p}$), we now require two necessary conditions to establish
stability: $\operatorname{Tr}J<0$ and $\operatorname{Det}J>0$, which are both
dependent on the values of $A_{p}$ and $t_{a}$. As before
$0>\operatorname{Tr}J=\frac{\lambda_{n}}{2}-\frac{1}{T_{n}}-\frac{\lambda_{p}A_{p}t_{a}}{2\nu}-\frac{1}{T_{p}}$
which simplifies to $\frac{1}{T_{n}}+\frac{1}{T_{p}}>\frac{1-A_{p}t_{a}}{2}$
if $\lambda_{n}=\lambda_{p}=\nu=1$. Note that if $A_{p}t_{a}\geq 1$, the
origin remains stable for all values of $T_{n},T_{p}$, which implies that
either the elastic modulus or the active time scale (or both) must be small to
ensure stability. The other necessary condition for stability
($\operatorname{Det}J>0$), for the case $\lambda_{n}=\lambda_{p}=\nu=1$, and
under the assumption that $T_{n}>2$, is $T_{p}>\frac{T_{n}-2}{A_{p}t_{a}}.$
This condition is shown as a pink line in Fig. 3 in a typical $T_{n}$-$T_{p}$
phase space when $A_{p}t_{a}<1$. Together the conditions now clearly dictate
the boundary of the stable region.

We numerically explore different combinations of $A_{p}t_{a}$ and confirm that
attractive orbits still exist as was found in the fixed $\nu_{p}$ case. As
before, oscillations occur when the trivial fixed point is unstable and when
no other singularity is present in the domain $\Omega$, as seen in the typical
phase diagram for $A_{p}t_{a}<1$ in Fig. 3, although oscillations are observed
for a possibly wider range $k_{1}<T_{n}<k_{2}$, where $k_{1},k_{2}$ are
dependent on $A_{p}t_{a}$. The nematic relaxation $t_{n}$ must remain
sufficiently longer than $t_{a}$, else the activity-induced shear dominates
suppressing the oscillations. The values of $T_{p}$ that permit oscillations,
in contrast to the fixed polymer viscosity case, can be much larger, owing to
the fact that polymer stress in Eq. (3) is now independent of $T_{p}$.

Figure 3: (a) Phase diagram for $A_{p}=0.5,t_{a}=1.5$ for
$\lambda_{n}=\lambda_{p}=\nu=1$. The red curve is $\operatorname{Tr}J=0$ and
the pink curve is $\operatorname{Det}J=0$. Crosses ($\times$) denote the
region where the only fixed point is the trivial one, whether stable or
unstable. Circles ($\circ$) denote the region where the trivial fixed point is
unstable. The red shaded region is where oscillations occur. The dashed line
is an approximate estimate of the boundary demarcating the region displaying
oscillations. (b) The oscillation period for the active time scale $t_{a}=1.5$
as a function of $T_{n}$ (blue, with $T_{p}=100$) or $T_{p}$ (red, with
$T_{n}=20$).

Note that it would not be possible to guarantee a trapping region in the
entire domain of $\theta_{n},\theta_{p}$, and the presence of several fixed
points will anyway prevent the use of the Poincaré-Bendixson Theorem. However
the Poincaré-Bendixson Theorem can again be used if we restrict the initial
conditions to a smaller domain, e.g. to $\Omega$. The resulting vector
diagrams and bifurcation diagrams are similar to the previous case. We also
note that oscillations occur right before the bifurcations appear
spontaneously in the bifurcation diagram, rather than the origin splitting
into a pair of fixed points. The calculations to establish the trapping
region, needed to utilise the Poincaré-Bendixson Theorem, are exactly the same
as in the previous case, thanks to the same terms vanishing.

###  III.3 Periodicity

Previous models and experiments have shown that the oscillations of active
matter in viscoelastic environments exhibit non-sinusoidal trajectories,
consistent with the behaviour of dynamical systems as they move farther from a
bifurcation point [23]. Here we characterize the oscillatory motion of our
model and confirm that a similar behaviour can be observed for the periodicity
and the shape of the orbits for different values of $T_{n}$ and $T_{p}$.

We first report results for fixed polymer viscosity. The period of the
oscillatory behavior, calculated numerically, increases nonlinearly as either
$T_{n}$ and $T_{p}$ increases (Fig. 4a). The active timescale $t_{a}$ simply
rescales the period for a fixed $T_{n}$-$T_{p}$ pair, as evidenced by the
perfect collapse of the measured periods.

Figure 4:  (a) The oscillation period for the fixed polymer viscosity case,
with active time scale $t_{a}=120$ and as a function of $T_{n}$ (blue, with
$T_{p}=55$) or $T_{p}$ (red, with $T_{n}=2.2$), obtained by rescaling
simulations from $t_{a}=20(\triangle)$, $t_{a}=30(+)$, $t_{a}=40(\times)$,
$t_{a}=60(\circ)$. (b) Points indicate the normalized discrete Fourier
transform of $\theta_{n}$ as a function of frequency; the dashed lines show
the sinusoidal signal of corresponding periodicity (with an exponential tail
arising from the numerical approximation). The parameters are
$t_{a}=60,T_{n}=2.2$ and $T_{p}=50$ (blue $+$), and $T_{p}=70$ (red
$\triangle$), and the values are normalized for each case by the maximum value
of the signal. The insets indicate the corresponding time series of
$\theta_{n}$ and orbits in the $\theta_{n}$-$\theta_{p}$ space.

As expected, the increase in the period is accompanied by a decreasingly
sinusoidal shape in the time series of the angles. For instance, the insets at
the top and bottom of Fig. 4b illustrate the evolution of the nematic
alignment angle $\theta_{n}$ for identical values of $T_{n}$. The blue curve
($T_{p}=50$) is barely indistinguishable from a sinusoidal curve whereas the
red curve ($T_{p}=70$) is reminiscent of the Fitzhugh-Nagumo or the Van der
Pol oscillators, which are characterised by alternating slow-fast relaxation
[31, 28].

To quantify the deviation of these trajectories from a sinusoid, we
numerically calculate their Fourier spectra and compare them to the spectra of
sinusoids with corresponding periods (Fig. 4b), with the peak in the
frequencies indicating the inverse of the period. Oscillations that result
from Eqs. (4) (in markers) indicate that the frequency distributions differ
from the discretized sinusoidal counterparts (dashed lines) with more
significantly-activated frequencies as $T_{p}$ or $T_{n}$ increases. This
irregularity can be visually corroborated by looking at the orbits in the
$\theta_{n}$-$\theta_{p}$ phase space as provided in the accompanying insets.
A sinusoidal trajectory would correspond to an elliptical shape, whereas the
slow-fast dynamics manifests as a rugby-ball-shaped orbit, with nearly pointed
tips. Moreover, we remark that the behaviour is closer to sinusoidal for pairs
$T_{n}$-$T_{p}$ near to the supercritical Hopf bifurcation (where
$\operatorname{Tr}J=0$, red solid line in Fig. 1a), consistent with the
findings in [23].

For the case of fixed elastic modulus, similar non-sinusoidal orbits are
observed. While the period increases with $T_{n}$ like in the fixed $\nu_{p}$
case, it decreases with $T_{p}$ and approaches a fixed value in this case
(Fig. 3b). This asymptotic decrease in period is attributed to the polymer
orientation being dominated by the fluctuating shear, which is independent of
$T_{p}$ in the fixed elastic modulus case. Slower polymer relaxation results
in the slow-fast transition being attained faster.

##  IV Conclusion

Our results help to elucidate the physics that underlies the existence of
oscillatory dynamics. We first note that our minimal model of an active
viscoelastic fluid displays oscillatory behavior when $T_{n}<T_{p}$. As such,
the nematic director relaxes and seemingly equilibrates faster than the
polymer, but polymer orientation attains larger values because of the stronger
tendency to align to the shear (weaker relaxation term). Since the shear rate
is an intricate balance between the active and polymeric stresses, then as the
nematic angle approaches its (non-interacting) steady state, it gets pushed
away from this equilibrium by feedback arising from the slower but stronger
polymer force. Indeed, if $\theta_{p}$ is sufficiently larger than
$\theta_{n}$, then $\sin 2\theta_{p}>(T_{p}\nu/\nu_{p})\sin 2\theta_{n}$ (for
fixed polymer viscosity, $\nu_{p}$) or $\sin 2\theta_{p}>(\nu/A_{p}t_{a})\sin
2\theta_{n}$ (for fixed elastic modulus, $A_{p}$) thereby creating a polymeric
stress large enough to flip the sign of the shear, $\dot{\gamma}.$ The nematic
angle $\theta_{n}$ then starts to realign to the new direction and overshoots
beyond 0, flipping the sign of the active stress and the intensifying the
shear rate in that opposite direction, resulting in a slow-fast dynamics.

To summarise, in this paper, we show how the viscoelastic components of a
fluid can feedback on the shear flow generated by an active component and
result in self-sustained oscillatory dynamics. The oscillations occur for
specific ranges of dimensionless numbers relating the relaxation timescales of
the active and viscoelastic components with the active timescale. In this
range of values, linear stability analysis shows that the trivial fixed point
of the system loses its stability, while numerical simulations verify that no
other fixed point exists in the domain spanned by the orientation angles of
the nematic and polymer, $\theta_{n}$ and $\theta_{p}$, respectively. Whether
the polymer viscosity or the elastic modulus is fixed, the Poincaré-Bendixson
Theorem can then be used to establish the existence of a limit cycle in (a
subset of) the two-dimensional phase plane and, consequently, the oscillatory
behavior in the shear rate and flow direction. For the case where the elastic
modulus $A_{p}$ is kept constant, further constraints on the value of the
product $A_{p}t_{a}$ are necessary to observe the oscillatory motion. The
period for a full cycle increases and the oscillations become less sinusoidal
as either the nematic or the polymer relaxation becomes slower. A slow-fast
dynamics is observed, with the faster dynamics occurring once $\theta_{n}$ has
flipped its orientation with respect to its equilibrium position. The physical
mechanism behind this phenomenon is explained by considering the relative
contributions of the active and polymeric stresses to the shear rate. We
emphasize that the nonlinearity in the model does not result in any known
oscillator models when truncated at any specific order.

Oscillatory dynamics is ubiquitous in biological systems, from simple
predator-prey models to the more complex descriptions of how neurons transmit
electrical signal spikes, such as the Hodgkin-Huxley model or its two-
dimensional counterpart, the Fitzhugh-Nagumo model [32, 31]. Indeed, recently,
a sustained oscillatory behavior of bacteria in a viscoelastic fluid was
reported and studied using both a full hydrodynamic active viscoelastic model
[23], and the Fitzhugh-Nagumo oscillator. Our two-dimensional model, motivated
by concepts from sheared active nematics and in spite of only accounting for
rotational relaxation and activity, captures the essential ingredients that
will result in oscillatory behaviour: activity to generate a spontaneous flow
and viscoelasticity that will generate sufficient amount of feedback. Our
study thereby may provide a minimal framework to predict and analyze generic
dynamical properties of biological systems, particularly those involving
nematic active matter in a complex environment.

## Acknowledgements

A. D. acknowledges funding from the Novo Nordisk Foundation (Grant no.
NNF18SA0035142), Villum Fonden (Grant no. 29476).

## Competing Interests

The authors declare no competing interests.

## Authors contribution

ELCMP, AD, and JMY conceptualized the study. ELCMP and HLT performed the
analyses and simulations. All authors participated the preparation and
submission of the manuscript.

## Verification of the trapping region

Consider the domain $\Omega$ and the vector field that arises from Eqs. (4).
The following observations taken together mean that any trajectory beginning
from the boundary remains confined within $\Omega$.

  * •

If $\theta_{n}=-\pi/4$, then, from Eq. 4a, $\cos 2\theta_{n}=0$ so
$\dot{\theta}_{n}=|\theta_{n}|/T_{n}>0$. Any trajectory starting from the left
side of the box $\Omega$ will point to the right.

  * •

If $\theta_{n}=\pi/4$, then, from Eq. 4a,
$\dot{\theta}_{n}=-\theta_{n}/T_{n}<0$. Any trajectory starting from the right
side of the box $\Omega$ will point to the left.

  * •

If $\theta_{p}=-\pi/4$, then, from Eq. 4b, $\cos 2\theta_{p}=0$ and so
$\dot{\theta}_{p}=|\theta_{p}|/T_{p}>0$. Any trajectory starting from the
bottom side of the box $\Omega$ will point upwards.

  * •

If $\theta_{p}=\pi/4$, then, from Eq. 4b,
$\dot{\theta}_{p}=-\theta_{p}/T_{p}<0$. Any trajectory starting from the top
side of the box $\Omega$ will point downwards.

We further remark that the same analysis but using (5)a,b in the domain
$\Omega$ results in the same observation that $\Omega$ is an external trapping
region.

## References

  * [1] Marchetti M C, Joanny J F, Ramaswamy S, Liverpool T B, Prost J, Rao M, et al. 2013  Hydrodynamics of soft active matter  Reviews of Modern Physics 85 1143–1189
  * [2] Doostmohammadi A, Ignés-Mullol J, Yeomans J M, Sagués F 2018  Active nematics  Nature Communications 9 1–13
  * [3] Sanchez T, Chen D T N, DeCamp S J, Heymann M, Dogic Z 2012  Spontaneous motion in hierarchically assembled active matter  Nature 491 431
  * [4] Kumar N, Zhang R, de Pablo J J, Gardel M L 2018  Tunable structure and dynamics of active liquid crystals  Science Advances 4 eaat7779
  * [5] Maroudas-Sacks Y, Garion L, Shani-Zerbib L, Livshits A, Braun E, Keren K 2020  Topological defects in the nematic order of actin fibres as organization centres of Hydra morphogenesis  Nature Physics 17 251–259
  * [6] Volfson D, Cookson S, Hasty J, Tsimring L S 2008  Biomechanical ordering of dense cell populations  Proceedings of the National Academy of Sciences 105 15346–15351
  * [7] Copenhagen K, Alert R, Wingreen N S, Shaevitz J W 2020  Topological defects promote layer formation in Myxococcus xanthus colonies.  Nature Physics 17 211–215
  * [8] Meacock O J, Doostmohammadi A, Foster K R, Yeomans J M, Durham W M 2020  Bacteria solve the problem of crowding by moving slowly  Nature Physics 17 205–210
  * [9] Duclos G, Blanch-Mercader C, Yashunsky V, Salbreux G, Joanny J F, Prost J, et al 2018  Spontaneous shear flow in confined cellular nematics  Nature Physics 14 728–732
  * [10] Kawaguchi K, Kageyama R, Sano M 2017  Topological defects control collective dynamics in neural progenitor cell cultures  Nature 545 327–331
  * [11] Saw T B, Doostmohammadi A, Nier V, Kocgozlu L, Thampi S, Toyama Y, et al. 2017  Topological defects in epithelia govern cell death and extrusion  Nature 544 212–216
  * [12] Blanch-Mercader C, Yashunsky V, Garcia S, Duclos G, Giomi L, Silberzan P 2018  Turbulent dynamics of epithelial cell cultures  Physical Review Letters 120 208101
  * [13] Nadell C D, Drescher K, Wingreen N S, Bassler B L 2015  Extracellular matrix structure governs invasion resistance in bacterial biofilms  The ISME Journal 9 1700–1709
  * [14] Hobley L, Harkins C, MacPhee C E, Stanley-Wall N R 2015  Giving structure to the biofilm matrix: an overview of individual strategies and emerging common themes  FEMS Microbiology Reviews 39 649–669
  * [15] Vidakovic L, Singh P K, Hartmann R, Nadell C D, Drescher K 2018  Dynamic biofilm architecture confers individual and collective mechanisms of viral protection  Nature Microbiology 3 26–31
  * [16] Carey S P, Martin K E, Reinhart-King C A 2017  Three-dimensional collagen matrix induces a mechanosensitive invasive epithelial phenotype  Scientific Reports 7 42088
  * [17] Mereness J A, Bhattacharya S, Wang Q, Ren Y, Pryhuber G S, Mariani T J 2018  Type VI collagen promotes lung epithelial cell spreading and wound-closure  PLOS ONE 13 e0209095
  * [18] Chaudhuri O, Cooper-White J, Janmey P A, Mooney D J, Shenoy V B 2020  Effects of extracellular matrix viscoelasticity on cellular behaviour  Nature 584 535–546
  * [19] Fletcher D A, Mullins R D 2010  Cell mechanics and the cytoskeleton  Nature 463 485–492
  * [20] Pritchard R H, Huang Y Y S, Terentjev E M 2014  Mechanics of biological networks: from the cell cytoskeleton to connective tissue  Soft Matter 10 1864–1884
  * [21] Hemingway E J, Cates, E M, Fielding S M 2016  Viscoelastic and elastomeric active matter: Linear instability and nonlinear dynamics  Physical Reviews E 93 032702
  * [22] Plan E L C V M, Yeomans J M, Doostmohammadi A 2020  Active matter in a viscoelastic environment  Physical Review Fluids 5 023102
  * [23] Liu S, Shankar S, Marchetti M C, Wu Y 2021  Viscoelastic control of spatiotemporal order in bacterial active matter  Nature 590 80–84
  * [24] Plan E L C V M, Yeomans J M, Doostmohammadi A 2021  Activity pulses induce spontaneous flow reversals in viscoelastic environments  Journal of the Royal Society Interface 18 20210100
  * [25] Giomi L, Mahadevan L, Chakraborty B, Hagan M F 2011  Excitable Patterns in Active Nematics  Physical Review Letters 106 218101
  * [26] Woodhouse F G, Goldstein R E 2012  Spontaneous Circulation of Confined Active Suspensions  Physical Review Letters 109 168105
  * [27] Aigouy B, Farhadifar R, Staple D B, Sagner A, Röper JC, Jülicher F, et al. 2010  Cell flow reorients the axis of planar polarity in the wing epithelium of Drosophila  Cell 142 773–786
  * [28] Strogatz S H 1994  Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, And Engineering (Massachusetts: Perseus Books)
  * [29] Perko L 2001  Differential Equations and Dynamical Systems (New York: Springer Verlag)
  * [30] Dano B 2021  StreakArrow.  MATLAB Central File Exchange.  Available from: https://www.mathworks.com/matlabcentral/fileexchange/22269-streakarrow
  * [31] FitzHugh R 1955  Mathematical models of threshold phenomena in the nerve membrane  Bulletin of Mathematical Biophysics 17 257–278.
  * [32] Hodgkin A L, Huxley A F 1952  Currents carried by sodium and potassium ions through the membrane of the giant axon of Loligo  The Journal of Physiology 116 449–472