Activity pulses induce spontaneous flow reversals in viscoelastic environments

Complex interactions between cellular systems and their surrounding extracellular matrices are emerging as important mechanical regulators of cell functions, such as proliferation, motility and cell death, and such cellular systems are often characterized by pulsating actomyosin activities. Here, using an active gel model, we numerically explore spontaneous flow generation by activity pulses in the presence of a viscoelastic medium. The results show that cross-talk between the activity-induced deformations of the viscoelastic surroundings and the time-dependent response of the active medium to these deformations can lead to the reversal of spontaneously generated active flows. We explain the mechanism behind this phenomenon based on the interaction between the active flow and the viscoelastic medium. We show the importance of relaxation time scales of both the polymers and the active particles and provide a phase space over which such spontaneous flow reversals can be observed. Our results suggest new experiments investigating the role of controlled pulses of activity in living systems ensnared in complex mircoenvironments.


Introduction
The study of biological systems as active materials has made tremendous advances in the past decades [1][2][3][4][5]. The 'activity' describes the ability of living systems to extract chemical energy from their surrounding environment and convert it into mechanical work. This happens at the level of individual constituents of the matter in, for example, sperm cells thrusting forwards by the rotation of their flagella, bacterial self-propulsion, eukaryotic cells migrating within extracellular networks and the cytoskeletal machinery inside cells that is powered by motor proteins. As such, the overarching theme in various kinds of living systems is the local activity generation that drives the entire system far from thermodynamic equilibrium, resulting in the collective patterns of motion observed in cellular tissues, bacterial colonies and sub-cellular flows [3,6,7].
The cross-talk between the mechanical micro-environment of living matter and this intrinsic ability of living systems to actively generate self-sustained motion governs pattern formation and self-organization in important biological processes including collective transport of sperm cells in confined tubes [8], shaping bacterial biofilms [9,10], tissue regeneration [11] and sculpting organ development [12]. Not only does the mechanical micro-environment provide geometrical constraints for active materials [13], but it is also often endowed with viscoelastic properties that allow for time-dependent responses to activity-induced stresses and deformations [2,14]. Significant examples are the extracellular matrices, surrounding cells and tissues that play a key role in cell death and proliferation, stem cell differentiation, cancer migration and cell response to drugs [15]. It is thus essential to understand the dynamic interconnection between the activity-induced stresses and the mechanical response of the viscoelastic medium.
Indeed, several recent studies have made the first attempts in this direction, showing that accounting for viscoelastic effects of the medium around living matter results in significant modification in the patterns of motion generated by continuous activity-induced stresses [16][17][18][19][20][21][22]. Importantly, however, in various biological contexts the activity generation is not constant and continuous, but rather is characterized by changes in the activity level and even activity pulses. Striking examples are the well-documented actomyosin contractility pulses that power the activity of epithelial cells and that have been shown to be essential in tissue elongation during development [23][24][25]. Therefore, here we examine the impact of activity pulses on the behaviour of active matter surrounded by a viscoelastic medium.
In order to investigate the fundamental impact of activity pulses, we employ a continuum model of active matter based on the theory of active gels, which has proven very successful in describing several aspects of the physics of active systems including actomyosin dynamics at the cell cortex [5,26], actomyosin-induced cell motility [27,28], actin retrograde flows [29,30] and the topological characteristics of actin filaments [31,32]. One important prediction of active gel models is the emergence of spontaneous flow generation in a confined active gel [33], which has been further experimentally validated in different biological systems [34,35] and is a generic feature of confined active materials. Here, we consider a simplified set-up of an active gel confined within viscoelastic surroundings and study the emergence of a spontaneous flow by introducing activity pulses and varying the relaxation time of the viscoelastic medium. We show that introducing activity pulses can result in the reversal of the spontaneous flow direction accompanied by the rearrangement of the orientation of active constituents. The mechanistic basis for this reversal is explained based on the feedback between the active flows and the viscoelastic deformation, particularly in between activity pulses. We further provide a simple model that reproduces the essential dynamics of the flow reversal and shows its dependence on the relevant time scales through a stability diagram.
The paper is organized as follows. In §2, we describe the details of the simulation set-up and introduce the governing equations of motion for the active gel, the surrounding viscoelastic medium and the coupling between the two. The results of the simulation together with the physical mechanism of flow reversal and its associated phase space are presented in §3. Finally, concluding remarks and a discussion of the broader impacts of the results are provided in §4.

Methods
The active gel is simulated in two dimensions as a horizontal stripe within a passive viscoelastic region, and is differentiated via an indicator function ϕ whose value is ϕ = 1 in the active region and ϕ = 0 in the passive viscoelastic region. The indicator function ϕ is only defined to distinguish between the active and passive regions and as such it is fixed, without any dynamical evolution. Activity and viscoelasticity are incorporated by introducing a generic two-phase model of active matter in viscoelastic domains [22,36], which is summarized below.

Active region
Following its success in describing the dynamics of the cell cytoskeleton, bacterial colonies and confluent cell layers, we use liquid crystal nematohydrodynamics to model the active region [4,37,38]. Within this framework, each active particle generates a dipolar flow field with the axis along its direction of alignment. The alignment direction is nematic, i.e. it has a head-tail symmetry. This can be captured at a coarse-grained level by the order parameter tensor Q through its principal eigenvector, which describes the nematic director orientation, and the associated eigenvalue, which describes the degree of alignment.
The free energy f Q for two-phase nematic models follows the Landau-De Gennes description where A Q describes the stability of the nematic or isotropic configurations, with the former being favoured when η > 2.7. The elastic coefficient K Q penalizes gradients in Q, and a positive (negative) L 0 enforces nematic orientation parallel ( perpendicular) to the active-viscoelastic interface.
In the presence of a velocity field u, the nematic tensor is evolved according to the equation where the left-hand side is the usual material advective derivative. The co-rotational term describes nematic reorientation in response to both vorticity Ω and flow strain D, with the tumbling parameter ξ determining whether the directors align or tumble in the flow. Γ Q controls the speed of relaxation towards the free energy minimum determined by the molecular field H Q = −δf Q /δQ. The typical nematic relaxation time scale t n when confined in a channel of width W is given by W 2 /Γ Q K Q and the dynamical equation for Q can thus be rewritten as Since individual components of the active region generate dipolar forces with the axis along their direction of alignment, the corresponding active stress is proportional to the orientation tensor [1,30,37] s active ¼ ÀzfQ, such that gradients in the orientation field generate forces on the fluid and drive active flows. The activity parameter, ζ, measures the strength of the active driving.

Viscoelastic region
The passive region is endowed with viscoelasticity that is described by the conformation tensor C, which characterizes the polymer orientation by its principal eigenvector, and the (square of the) polymer length, by its trace. Here we use the Oldroyd-B model to reproduce simple viscoelastic effects, i.e. polymer relaxation is linear with respect to the elongation and is governed by a single relaxation time τ [39]. The free energy associated with an Oldroyd fluid, governs the polymer relaxation to its unstretched equilibrium C = I (where I is the identity tensor). Here, the modulus of elasticity A C = ν/τ is the ratio of the polymer contribution to viscosity ν to the polymer relaxation time τ. The corresponding molecular field H C = −δf C /δC appears in the dynamical equation governing the evolution of C: , whered enotes the matrix transpose. Similarly to equation (2.2) for the orientation tensor Q, the evolution of C accounts for the advection of the polymer conformation C, its response to velocity gradients through S C ¼ CV À VC þ CD þ D`C`and its relaxation to royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210100 equilibrium at a rate Γ C = ν −1 . The viscoelastic time scale is clearly seen when this equation of motion is simplified as The polymer contribution to the stress, within the assumptions of the Oldroyd-B model, is where the factor 1 − ϕ ensures that the polymer stress only acts within the passive region (where ϕ = 0).

Coupling and simulation details
The active and viscoelastic regions interact with each other through the incompressible velocity field u, which obeys the Cauchy momentum equation where ρ is the fluid density, p is the pressure and σ is the sum of viscous, capillary, elastic, active (equation (2.3)) and polymer (equation (2.6)) stresses. The other stresses take the forms [22]: Here, f is the total free energy, μ = −δf/δϕ is the molecular potential and K f is a parameter controlling the active-polymeric interface width. Both the active and viscoelastic regions exert stresses on the fluid, and the resulting velocity field couples the two regions through advective and co-rotational terms in equations (2.2) and (2.5).
Equations (2.2), (2.5) and (2.7) are evolved using a hybrid lattice Boltzmann method [40,41]. The simulation domain has dimensions L × H, is periodic in the x-direction and has noslip boundary conditions in the y-direction (see schematic in figure 1a), as in typical experimental and numerical studies of active matter (e.g. [14]). The dynamics is not affected by the length of the channel L because of periodicity; our results here use L = 10, H = 100. The width of the active region is fixed at W = 20.
In order to obtain a unidirectional flow, the parameters for the active region are chosen to be A Q = 1.0 and K Q = 0.2, and we use ξ = 0.7 corresponding to the flow-aligning regime as previous studies of confined active nematics have shown that a unidirectional shear-like flow of active nematics in a confined channel can only be achieved for flow alignment [42]. η(ϕ) is chosen such that η = 2.95 within the active region ϕ = 1 [42], while the passive viscoelastic region, ϕ = 0, is in the isotropic royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210100 phase with η < 2.7. The polymer contribution to viscosity is fixed at ν = 1 and A C = 1/τ. The values of the flow parameters are ρ = 1, p = 0.25, and the Newtonian contribution to the viscosity is ν flow-= 2/3. We also enforce weak nematic anchoring at the boundary, L 0 = 0.05, for stability. This was verified to have no qualitative effect on the flow reversal dynamics. Each simulation begins with equilibrated polymers C = I and nematic directors with small random perturbations about the x-axis. Active stresses are applied until the system establishes a steady-state, unidirectional flow, after which the activity is temporarily turned off at t off = 10 5 for a duration of d = t on − t off time steps. The equations are solved until the flow is re-established.

Results
When activity is turned off at a time t off , and turned back on at a time t on , the flow is re-established in the same or, surprisingly, in the opposite direction. This is not a random choice but depends sensitively on the parameters setting the relevant nematic and viscoelastic time scales. For example, the movie provided as electronic supplementary material shows several successive reversals in the direction of the velocity. We first explain how this dependence comes about, and then present a simple model which illustrates the underlying physics. The mechanism of the flow reversal is summarized in figure 1b-d. Figure 2 shows the variation of the mean velocity in the channel, 〈u x 〉, and the angle θ n that the mean director field at the interface makes with the channel axis as a function of time for selected simulation parameters. Consider first figure 2a where there are no polymers in the passive region. Activity is switched on at time t = 0. This drives the active nematic instability and active stresses resulting from gradients in the director field set up a linear flow along the stripe. The spontaneous flow is stabilized by the channel interfaces, and the resulting steady-state flow profile corresponds to directors aligning at the Leslie angle to the local shear [43,44]. When the activity is switched off the flow royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210100 decays faster than the relaxation time of the nematic director. In the absence of polymers, if activity is switched back on before the end of the decay the residual director rotation ensures that the flow re-starts in the same direction as before.
The presence of polymers changes the response of flow to activity pulses. In figure 2b-d, the passive medium is viscoelastic and the shear flow induced by the activity stretches the polymers, thus storing energy. The initial build-up of the flow is similar to the no polymer case, but slower, as work is done to stretch the polymers, and the steady-state value of the flow in the channel is lower because the polymers impose stresses that oppose the flow. This also means that after t off the decays of the velocity and director fields to equilibrium are faster and, in particular, stress imposed by the relaxation of the surrounding polymers may be strong enough to reverse the flow in the channel.
For example in figure 2b both the residual flow and the director have reversed at t on , making a velocity reversal inevitable. By comparison, for the same polymer relaxation time but faster switching of activity in figure 2c, at t on the director distortion is still just positive but the velocity has been reversed by the polymer stresses. The activity tries to rebuild the flow in the original direction whereas the remaining elastic energy stored in the surrounding polymers pushes the fluid in the opposite direction. The flow slowly reverses as the instability is ( just) overcome by the residual polymer stresses. The director distortion reverses and the velocity slowly increases until it eventually attains its steady-state value, but in the opposite direction. By contrast, figure 2d shows an example of faster polymer relaxation, compared with figure 2c, where the polymer stress is not sufficiently strong to cause flow re-orientation and both nematic directors and mean velocity regain their original direction.
Such a flow reversal mechanism that depends on stored polymer stresses and residual director orientation implies that the relevant time scales here are the polymer relaxation time τ, the nematic relaxation time t n and the period of inactivity d = t on − t off . To further explain the phenomenon of flow reversal and highlight the competition between these varying time scales, we construct simplified, space-independent, dynamic equations for the evolution of the active nematic and polymeric particles.
To this end, we consider the dynamics of the alignment for the angles θ n , θ p formed by the active nematic directors and the polymer directors, respectively, with respect to the direction of a simple shear flow, u ¼ ( _ gy, 0) [45], These equations are approximations of the orientation of rodlike particles with tumbling parameters ξ n , ξ p in response to a shear flow with rate _ g. In the absence of shear, _ g ¼ 0, the angles will relax exponentially to zero.
A time-dependent shear rate _ g couples the two equations and can be determined from balancing the viscous stress, σ viscous = 2ν flow D, with the active (equation (2.3)) and polymer (equation (2.6)) stresses. Solving for the off-diagonal term of the rate of strain tensor D and constructing the tensors Q ij = n i n j − δ ij /2 and C ij = p i p j from the unit directorsn ¼ ( cos u n , sin u n ) (and analogously for p), we obtain the simple time-dependent shear in terms of the nematic and polymer directors. For simplicity, we fix ν flow = 2/3, ν = 1 and set ξ n = 1.1 and ξ p = 0.275 and take initial conditions θ n = θ p = 0.01. It is clear from (3.3) that, when the activity ζ = 0 is turned off, the flow reverses in the opposite direction owing to the presence of the polymers. Moreover, a simulation of this simpler set of equations features the same essential reason for the reversal as in the full simulations: there exists sufficient polymer stress to reorient the nematics during the period between activity pulses. Figure 2e-g shows that the nematic angle θ n exhibits trajectories similar to the result from the full equations of motion, figure 2b-d. The reversed flow can drive the nematic angle negative during the period of inactivity, which is a sufficient condition for flow reversal. Moreover, even if θ n remains positive when the activity is switched on, the polymer stress can still overcome the activity (i.e. ζ sin 2θ n < ν sin 2θ p /τ) and continue to drive the nematic director to reverse direction.
Using the simple model, we examined the phase space of flow reversal for varying time scales of nematic and polymer relaxations, t n and τ, respectively, as well as the activity switching time d = t on − t off (figure 3a). In particular, we observe that the reversal can be ensured if three dimensionless ratios are sufficiently large: (i) d/t n : nematic directors have enough time to relax. (ii) τ/t n : polymers retain enough energy to reverse nematic orientation when θ n ≈ 0. (iii) dt=t 2 n : residual polymer stress at the moment of reactivation can help reverse nematic orientation if condition (i) alone is insufficient.
These three constraints are illustrated as dashed black lines on the flow reversal phase diagram in the d/t n −τ/t n phase space and clearly distinguish the parameter space with (blue markers) and without (red markers) flow reversal.
We also plot, in figure 3b, the similar phase diagram obtained from full simulations of the active gel surrounded by a viscoelastic medium showing close qualitative agreement with the simple model, and indicating that the parameter region for flow reversal is suitably captured by these three constraints. The quantitative difference between the simplified model and the full simulations is expected because in the simplified model there is no space dependence and we only account for polymer orientation in the polymer stress, whereas the full simulations are space dependent and have a polymer stress (equation (2.6)) that depends on polymer elongation, which is higher when τ is larger. Notwithstanding these limitations, the close qualitative agreement of the angle profiles and the phase diagram obtained from the simple model with those from the full hydrodynamic simulations of spatio-temporal evolution of the active nematics and polymeric fluids confirms the importance of the time-scale constraints for the flow reversal and the underlying physics of stress balance during the period of inactivity.

Conclusion
In this article, we present how, in the absence of any external forcing, activity pulses in living matter interacting with a viscoelastic environment can spontaneously generate flow reversals. Based on a well-documented active gel theory, a spontaneous steady flow of active matter is achieved even while in contact with polymer-laden surroundings. This flow stretches the polymers near the interface, which, in between periods of activity, relax and produce a weak backflow that may determine the flow direction upon resumption of the active driving. While this work is done in two dimensions, we expect a similar phenomenology to occur in a three-dimensional set-up.
The well-established spontaneous flow generation of confined active matter relies on the level of activity. Here the spontaneous flow reversals hinge on several time scales: the polymer relaxation time, the interval between activity pulses and the relaxation dynamics of nematic active matter, as shown in the phase diagram. Indeed the need for sufficient polymer stretching and feedback as well as quick nematic reordering highlight the time-dependent viscoelastic response in between activity pulses. Our work emphasizes not only the importance of accounting for a viscoelastic environment but also the involvement of several time scales arising from both active matter and its surroundings.
Flow reversals have recently been observed in Escherichia coli bacteria swimming in DNA polymer solutions [14] in the absence of pulsatile activity. Our theoretical work invites several experiments to be performed in order to confirm our mechanistic predictions, provide a possible mapping to the parameters and time scales, and deepen our understanding of the role of a viscoelastic environment in the dynamics of living matter. For instance, cells confined in a monolayer have been shown to display nematic behaviour and exhibit unidirectional flows in channels [46]. They could be repeatedly subjected to cell inhibitors or uncouplers [47][48][49] to simulate gaps between activity pulses. Another option would be to regulate the movement of a colony of elongated bacteria by utilizing phototactic methods [50]. Our result thus holds potential in understanding mechanotaxis and motivating the use of viscoelastic media to control living matter at microscopic scales.
Authors' contributions. All authors contributed to the conceptualization, data analysis and preparation and submission of the manuscript. royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210100