Optimal Control of Short-Time Attractors in Active Fluid Flows

Long-term material transport is a fundamentally Lagrangian phenomenon, originally studied by keeping track of the longer-term redistribution of individual tracers released in the flow. In that setting, the Finite Time Lyapunov Exponent (FTLE) and Lagrangian coherent structures (LCSs) have been efficient predictors of tracer behavior see e.g. [haller2015lagrangian, shadden2005definition, hadjighasem2017critical, serra2017uncovering, serra2020dynamic, serra2021defect]. An alternative to Lagrangian approaches is to find their instantaneous limits purely from Eulerian observations, thereby avoiding the pitfalls of trajectory integration, while predicting short-term material transport. Additionally, LCSs are generally impractical to control – no literature exists – because they are defined as nonlinear functions of fluid trajectories, which are themselves integrals of the Eulerian velocity $\cv{v}$.

Short-time attractors – originally defined as Attracting OECSs [serra2016objective] – govern material transport in fluid flows over short-times, revealing critical information in challenging problems such as search and rescue operations at sea [serra2020search] and oil-spill containment [duran2021horizontal]. A simpler, more controllable alternative to attracting OECSs for locating short-time attractors is the instantaneous Lyapunov Exponent (iLE)[Nolan2020], defined as the instantaneous limit of the well known FTLE. The iLE locates short-time attractors as trenches – or negative regions – of the smallest eigenvalue $s_{1}$ of the rate-of-strain tensor of the fluid velocity $\cv{v}$. \textcolorblackFor example, Figs.1a-b show short-term attractors marked by trenches of $s_{1}$ (scalar field) in an experimental velocity field (black vector field) describing the motion of thousands of cells during chick gastrulation [Rozbicki2015]. A strong trench of $s_{1}$ marks a short-term attractor along the anterior-posterior (AP) axis corresponding to the forming primitive streak [serra2020dynamic] (panel b), while remaining not identifiable from the inspection of the corresponding velocity field (panel a and its inset). Similar results hold in different flows (see e.g., Fig. 1 of [serra2016objective] and Figs. 4-5 of [serra2020search]). To confirm the effect of short-term attractors, panel b shows that tracers (green dots) released from a uniform grid of initial conditions and advected by $\mathbf{v}(\mathbf{x},t)$ for a short time accumulate on the $s_{1}$ trench.

blackClassical asymptotic attracting structures (Fig. 1c) include sink-type fixed points (red dot) and the unstable manifold (red curve) of saddle-type fixed points of steady velocity fields. These structures, however, cannot move or deform in time and influence material transport only for steady velocity fields. \textcolorblackAs shown in Figs.1a-b, inspection and control of the velocity field $\mathbf{v}$ is sub-optimal to create material traps in general unsteady flows. First, because the velocity field and its streamlines are not objective, i.e., they depend on the choice of reference frame used to describe motion (see also SM Section 6 and Figs. S1-S2). By contrast, the location of material accumulation is frame invariant [haller2015lagrangian, serra2016objective], as any quantification of the material response of a deforming continuum must be according to a fundamental axiom of mechanics [gurt]. Second, it might be an unnecessarily strong requirement, or uncompliant with boundary conditions, to prescribe $\mathbf{v}(\mathbf{x},t)$ directly. \textcolorblackOn the other hand, short-time attractors (Figs. 1b,d) are frame invariant, can move and deform, and apply to general unsteady flows. Here we devise an optimal control framework that uses the active stress intensity $m$ as the control input to achieve a target short-time attractor, defined by a desired distribution of $s_{1}$ or iLE, in active viscous flows.

We adopt a simplified version of the mechanochemical model developed in [serra2021mechanochemical] consisting of an active stokes flow characterized by the passive viscous stress $\mathbf{\sigma}_{v}=-p\mathbf{I}+2\mu\mathbf{S}_{d}$ and active stress $\mathbf{\sigma}_{a}=m(\mathbf{B}-\mathbf{I}/2)$, where $\mathbf{S}_{d}$ is the deviatoric rate-of-strain tensor, $\mathbf{B}=\cv{e}\otimes\cv{e}$ characterizes the orientation of active elements $\cv{e}=(\cos(\phi)\,\sin(\phi))^{\top}$, $m$ denotes the intensity of active stress and $\mathbf{I}$ the identity tensor. To account for flow compressibility, we use a simple continuity equation $\mathbf{\nabla\cdot v}=c(-2p-p_{0}m)$ where positive isotropic viscous stress ($p>0$), and isotropic contractile-type ($m>0$) active stress contribute to negative flow divergence via the bulk viscosity $1/c$ and a nondimensional parameter $p_{0}$. Biologically, $p_{0}$ modulates the cell propensity to ingress into the third dimension given active isotropic apical contraction. The resulting system of PDEs in nondimensional form [serra2021mechanochemical] is

where $p_{1}=\mu c$ is a second nondimensional parameter characterizing the ratio of the shear to bulk viscosity, and $\bb{A}=2p_{1}\mathbf{B}+p_{1}(p_{0}-1)\mathbf{I}$. The first two terms in the force balance describe passive forces due to shear and compressibility, while the active forces $\mathbf{g}$ arise from spatial variations of $m$ and $\phi$. To set up our optimal control problem, we consider a simplified dynamics where the orientation of active elements $\phi(\bb{x})$ is time-independent and prescribed, while the active stress intensity $m(\mathbf{x},t)$ is the control input.

s a consequence, the state equation has the same structure found in linear elasticity theory i.e. Navier-Cauchy equations which, however, describe the displacement dynamics.

To control short-time attractors, the OCP involves steering the minimum eigenvalue of the rate of strain tensor towards a target function while minimizing the overall control effort and its gradient:

where $z(\mathbf{x})$ represents a scalar target for the minimum eigenvalue $s_{1}(\mathbf{x})$ of the rate-of-strain tensor $\mathbf{S}$ of $\cv{v}$, and $\cv{d}(\mathbf{x},t)$ is a disturbance on the state dynamics.

Following a Lagrangian variational approach, we derive a system of first-order necessary optimality conditions for the OCP problem. We note that $J_{s}$ is nonlinear and nonquadratic due to the relation $s_{1}(\mathbf{S}(\cv{v}))$. Consequently, sufficiency is not guaranteed, as typical in nonconvex problems. We define the Lagrangian as

where $\cs{\lambda}$ is the adjoint function, and obtain the strong form of the optimality conditions by imposing the first variations of the Lagrangian to be zero for all allowed variations of the state, adjoint and control functions. We provide the complete derivation in the SM Sections 2-4, and summarize here the necessary conditions as a coupled system of PDEs. The optimal pair for the OCP \eqrefOCP_l $(\cv{v},\,m)$ should satisfy the following system of first-order necessary conditions

where $\cs{\xi}_{1}(\bb{x})$ is the eigenvector function associated to $s_{1}(\bb{x})$. We note that \eqrefnc requires an iterative method due to the nonlinear relationship between $\cv{v}$ and the forcing term of the adjoint equation involving $s_{1}$ and $\cs{\xi}_{1}(\bb{x})$.

We numerically solve \eqrefnc using the Finite Element Method (FEM) and a gradient-based algorithm (SM Sec. 5). We consider a circular domain with zero velocity boundary conditions and note that our algorithm applies to arbitrary-shaped domains. We set the target shape $z$ as a scaled indicator function of a rectangle so that the target value is $-10$ inside the rectangle and zero elsewhere. We set the cable orientation to a constant value $\phi=\frac{\pi}{4}$ from the x-axis and choose the control weighting parameter $\beta=10^{-6}$ and the nondimensional model parameters $p_{0}=10,\ p_{1}=0.5$. $p_{1}$ modulates the overall fluid compressibility while high $p_{0}$ induces high negative divergence in regions with higher $m$ [serra2021mechanochemical]. We select the space-time varying disturbance force as \textcolorblack$\cv{d}(\mathbf{x},t)=d\,e^{-(\frac{r(t)}{\sigma})^{2}}\,[-(y-y_{c}(t)),x-x_{c}(t)]$, where $\cv{x}_{c}(t)=[-0.5+t,0.5]$, $r(t)=||\cv{x}-\cv{x}_{c}(t)||$, and set the intensity $d=50$, and standard deviation $\sigma=0.2$. Figure 2 shows the resulting optimal $s_{1}$ along with a grid of particles advected over short times for two different initialization times. Figura 4:

Figura 3:

Figure 3 shows the optimal state-control pair and its associated disturbance forcing $\cv{d}$. The control $m$ acts through $\mathbf{g}(m)=\nabla\cdot(\mathbf{A}m)$, and therefore both $m$ and $\mathbf{\nabla}m$ contribute to the state dynamics. This is why the control generates sharp gradients to track the target eigenvalue accurately. Furthermore, the control field rejects the disturbance by generating a vortex-shaped counteracting effect. Overall the disturbance strongly influences the optimal velocity. Figure 4 shows the interplay between the control weight $\beta$, the disturbance $\mathbf{d}$ and the accuracy of the tracking objective that is key to generating a short-time attractor. As in the earlier figures, each row (a-b) corresponds to a different time ($t_{i},\ i=1,2$).

To present an additional test case, we select a different disturbance compared to Figs. 2-3, generating two vortex-shaped force-field streamlines (Figure 4, first column) using the same functional form in the previous test case. The target eigenvalue $z$ is the same as in Figs. 2-3. The uncontrolled dynamics ($m=0$) does not generate attraction, as shown by short-term advected particles (yellow dots) and the $s_{1}(\mathbf{x})$ field in Figure 4, second column. A high $\beta=0.1$ (or weak control) steers the eigenvalue towards the target, but its effect is too mild to generate short-time attraction (Figure 4, third column). Indeed, a higher control weight would result in a lower relevance of the tracking term in the cost function and hence poor tracking performances. Figure 4 fourth column shows the results for a low $\beta=10^{-6}$ (or strong control), as in Figure 2, resulting in the desired short-time attractor that serves as a material trap while rejecting disturbances. Different rows (a-b) show the same results for two different initialization times as in Figs. 2-3.

figure Optimal fields for $\phi=\frac{\pi}{4}$, $p_{0}=5$,$p_{1}=0.5$ and $\beta=10^{-3}$ where the target $z$ is a combination of radial basis functions with zero mean. The broken symmetry in the top fields is due to the directionality imposed by $\phi$.

figure Optimal state and control norm as a function of the control weighting $\beta$ and the parameter $p_{0}$. In the $\beta$ sweep $p_{0}=5$ while in the $p_{0}$ sweep $\beta$ is set at $10^{-3}$. \crlCapire il minimo di norma di controllo per $p_{0}=4$

By combining recent advances in coherent structures from nonlinear dynamics, control theory and active matter, we formulated, analyzed and numerically solved an optimal control problem that generates material short-time attractors at desired locations while rejecting disturbances. We used a simplified version of model [serra2021mechanochemical] describing a compressible, viscosity-dominated, active fluid representing a planar multicellular flow, and used the active stress intensity as the control input. We identify short-time attractors as the smallest instantaneous Lyapunov exponent [Nolan2020], the instantaneous limit of the well-known FTLE, which is computable from the smallest eigenvalue of the rate of strain tensor of the flow velocity.

Short-time attractors are frame invariant and predict the correct location of material attraction, which may be undetected from the inspection of the frame-dependent velocity field (Fig. 1, [serra2016objective, serra2020search]). Using the same framework, one can control the position of material repellers, which, together with attractors, shape complex motion in synthetic active matter [serra2021defect] and embryo morphogenesis across species [serra2020dynamic, lange2023zebrahub]. Controlling attractors and repellers offers a new, robust, and frame-invariant perspective to steer motion in synthetic and natural active matter. It will enable the creation of material traps for medical applications as well as enhance our understanding and the possibility of manipulating multicellular flows in morphogenesis. For example, it will shed light on how myosin activity (active stress intensity) generates the required motion that compartmentalizes the embryo, segregating distinct cell types (repellers) and steering specific cells to precise locations (attractors). In future work, we plan to consider the explicit orientational dynamics of the active stress anisotropy, the effect of inertial forces, and the control of Lagrangian Coherent Structures that shape fluid motion over longer time scales.