Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

Saptorshi
Ghosh
^{a},
Chaitanya
Joshi
^{b},
Aparna
Baskaran
*^{a} and
Michael F.
Hagan
*^{a}
^{a}Martin Fisher School of Physics, Brandeis University, Waltham, Massachusetts 02453, USA. E-mail: aparna@brandeis.edu; hagan@brandeis.edu
^{b}Department of Physics and Astronomy, Tufts University, Medford, Massachusetts 02155, USA

Received
8th May 2024
, Accepted 12th August 2024

First published on 14th August 2024

We apply optimal control theory to a model of a polar active fluid (the Toner–Tu model), with the objective of driving the system into particular emergent dynamical behaviors or programming switching between states on demand. We use the effective self-propulsion speed as the control parameter (i.e. the means of external actuation). We identify control protocols that achieve outcomes such as relocating asters to targeted positions, forcing propagating solitary waves to reorient to a particular direction, and switching between stationary asters and propagating fronts. We analyze the solutions to identify generic principles for controlling polar active fluids. Our findings have implications for achieving spatiotemporal control of active polar systems in experiments, particularly in vitro cytoskeletal systems. Additionally, this research paves the way for leveraging optimal control methods to engineer the structure and dynamics of active fluids more broadly.

Recent experimental advances have put the objective of control within our reach. Experiments have demonstrated that light can be used as a control field to assemble self-limited functional structures in active colloids^{39,40} and to exert spatiotemporal control of motility induced phase-separation (MIPS).^{41,42} By shining sequences of light that vary in space and time on active materials constructed with light-activated motor proteins, researchers can control the average speed of active flows^{43,44} and steer defects in active nematics,^{45} and drive the formation and movement of asters in isotropic suspensions.^{46}

Theoretical progress toward functionalizing active matter has taken two paths. The first has been to impose spatiotemporal activity patterns and observe their effects on the system dynamics,^{45,47–49} thereby identifying easily accessible target states. The second is to use the framework of optimal control^{50,51} and optimal transport^{52} to identify spatiotemporal activity patterns that will drive the system to a pre-chosen target dynamics.^{47,53,54}

The work described in this article belongs in this second class. We study the optimal control theory of the classic active matter theory, a dry 2D active polar fluid, first considered by Toner and Tu.^{55,56} We treat the convective speed of the active particles as the control parameter and identify control solutions that drive the system to targeted steady states, including forcing an aster to move to a given spatial location, causing propagating stripes to reorient along a particular direction, and driving a propagating stripe to convert into a stationary aster. In the previous work most closely related to this one, Norton et al.^{53} obtained a control solution to switch a confined active nematic between two symmetric attractors (changing handedness of a circular flow state). Here, we show that an optimal control framework can solve much more diverse problems, including switching among attractors with very different broken symmetry patterns and driving the system into states which are not stable attractors at a given set of parameters. Further, we analyze the identified optimal solutions to uncover generic control principles that are broadly applicable to active polar fluids.

This paper is laid out as follows. In Section II, we review the hydrodynamic theory and describe the key features of the steady states that arise in the absence of control. Section III describes the method for implementing optimal control theory. In Section IV, we report the results of the control solutions for driving the system to particular steady states or switching between them. Then, we analyze the control solutions in terms of the dynamical equations of motion to identify the essential mechanisms that drive the system into the desired behavior, with the goal of identifying generic principles. We also investigate how robust the optimal control solutions are to errors arising from experimental noise and inaccuracies in model parameters. Finally Section V concludes with a discussion on testing these results in experiments.

∂_{t}ρ = −∇·(ωτ − D∇ρ) | (1) |

∂_{t}τ + λ_{1}τ·∇τ = −ν(a_{2}(ρ) + a_{4}(ρ)|τ|^{2})τ − ∇(ωρ) + K∇^{2}τ + λ_{2}τ_{α}∇τ_{α} + λ_{3}τ∇·τ. | (2) |

The density dynamics eqn (1) is an advection-diffusion equation where the advective velocity is proportional to the orientational order parameter P(r,t). The dynamics of the orientational order has three features: (i) It has a self-convection term with coefficent λ_{1}, encoding the absence of Galilean invariance in the ‘dry’ theory.^{55,56} (ii) It has terms consistent with model A dynamics that drive the system downhill on a free energy landscape where , encoding the fact that the flocking occurs due to spontaneous symmetry breaking. (iii) It contains a hydrostatic pressure term of the form , where λ_{2} encodes the tendency for elongated self-propelled particles to splay due to recollision events.^{62,66,67}

In this study of spatiotemporal control, we choose a_{2}(ρ) = (1 − ρ/ρ_{c}) and a_{4}(ρ) = (1 + ρ/ρ_{c})/ρ^{2}, yielding a continuous mean-field transition from an isotropic τ to a homogeneous, polar or a swarming state (|τ| > 0) at the critical density ρ = ρ_{c}. For simplicity, we choose D = K. Without loss of generality we set the critical density ρ_{c} = 1. We further reduce the number of independent parameters and set λ_{1} = λ_{2} = λ_{3}. We set the unit of time as τ_{0} = ν^{−1}, the relaxation time scale of the orientation field, and the unit of length as . The simplified dimensionless equations are:

∂_{t}ρ = −∇·(ωτ − ∇ρ) | (3) |

∂_{t}τ = −(a_{2}(ρ) + a_{4}(ρ)|τ|^{2})τ − ∇(ωρ) + ∇^{2}τ + λ(τ_{α}∇τ_{α} + τ∇·τ − τ·∇τ). | (4) |

The phenomenology of this model is described in ref. 62 and summarized in Fig. 1a. For the purposes of this work, we note that the dynamics of this system admits two inhomogeneous steady states: (i) propagating stripes composed of ordered swarms moving through a disordered background, which are referred to as polar drops elsewhere in the literature,^{63,71} when ω/λ ≫ 1, and (ii) a stationary high density aster, again in an isotropic background, when λ/ω ≫ 1. Both of these states arise close to the threshold density for orientational ordering (which we set to ρ_{c} = 1), and correspond to the system phase separating into a dense ordered phase and a dilute disordered phase. In this study, we fix the homogeneous density at ρ = 1.07 and consider the problem of controlling the inhomogeneous steady states using spatiotemporal patterning of the convective strength ω, referred to as the activity or control parameter in the rest of this paper.

This model has the virtue of mathematical simplicity and a minimal number of parameters, thus enabling physical insight from the optimal control solutions. At the same time, the states we seek to control are directly realizable in experiments (see Section V). Further, in the ESI,†^{} we show that this theory is linearly controllable on short length scales. Thus it is ideal for investigating applications of optimal control in active systems.

(5) |

The terms (ρ − ρ*)^{2} and (τ − τ*)^{2} measure the deviations from the target state, ρ* and τ*. We constrain our search of optimal state trajectories to those that obey the system dynamics by introducing Lagrange multipliers, η and ν, which are adjoint variables for ρ and τ. These dynamical constraints are enforced in the optimization by adding them to the cost function as

(6) |

The necessary condition for optimality is ,^{72,73} so , , , , , , . The first two conditions return eqn (3) and (4) governing ρ and τ. The following two conditions yield the dynamical equations for the adjoint variables η and ν,

∂_{t}η = C(ρ − ρ*) − ∇^{2}η − ^{2}∇·ν + (δ_{ρ}a_{2}(ρ) + δ_{ρ}a_{4}(ρ)|τ|^{2})(ν·τ). | (7) |

∂_{t}ν = D(τ − τ*) + (a_{2}(ρ) + a_{4}(ρ)|τ|^{2})ν + 2a_{4}(ρ)τ(ν·τ) − λ[−τ∇·ν + 2ν∇·τ − ∇(ν·τ) + τ·∇ν − ν_{α}∇τ_{α}] − ∇^{2}ν − ^{2}∇η. | (8) |

2A(^{2} − _{0}^{2}) − 2B∇^{2}^{2} − 2K(d^{2}^{2}/dt^{2}) − 2τ·∇η − 2ν·∇ρ = 0. | (9) |

We use the direct-adjoint-looping (DAL) method^{74} to minimize the cost function under the constraint that the dynamics satisfies eqn (3) and (4), to yield the optimal schedule of activity in space and time (see ESI† Sections SIV and SV for more details). Specifically, we construct an initial condition by performing a simulation with unperturbed dynamics (eqn (3) and (4)) until reaching steady-state, at a parameter set that leads to a desired initial behavior. We construct a target configuration in the same manner, using a different parameter set that leads to the desired target behavior. We also specify a time duration t_{F} over which the control protocol will be employed, and an initial trial control protocol. We then perform a series of DAL iterations; in each iteration the system and the adjoint dynamics are integrated from the initial condition for time t_{F} under the current control protocol, and the cost function (eqn (5)) is computed from the resulting trajectory. The adjoint equations are integrated backwards in time to propagate the residuals. After each backward run, the control protocol is updated via gradient descent, , to minimize the cost function. We employ Armijo backtracking^{75} to adaptively choose the step-size for gradient descent and to ensure convergence of the DAL algorithm. Once the optimal spatiotemporal activity pattern has been computed, performing a forward integration using eqn (3) and (4), the same initial condition, and ω = ^{2} will yield the solution trajectory. We have implemented this entire calculation in the open-source Python finite element method solver FEniCS.^{76} We have provided the associated source code and data here.

Fig. 2 summarizes the results of this computation. Fig. 2a and b respectively show the time evolution of the system configuration and the applied control field that drives the transformation. At early times, the applied control is strongest at the aster core while it is lowest in front of the aster along the direction we seek to move it. The aster then elongates to assume a comet-like shape, with a denser, polar-ordered head (see snapshots at t = 50, 100), as it advects toward the target location. Note that the activity is largest to the rear of the aster in this region. Thus, the control solution pushes (rather than pulls) the aster.

Fig. 2 Advecting an aster. The control solution for moving an aster from (x_{I} = 60, y_{I} = 40) at t = 0 to (x_{F} = 60, y_{F} = 70) at time t_{F} = 2000. All length and timescales are presented in dimensionless units, which are defined in Section II. (a) and (b) Snapshots of (a) the density ρ (color map) and polarization τ (arrows) profiles and (b) the control solution activity field ω (color map). The ‘x’ symbols in the snapshots in (a) and (b) show the position of the aster core. At times t = 50, 100 the upper quadrant of the aster unwinds and the aster becomes prolate while the aster core maintains the +1 defect. At t = 2000 the aster is reformed at the target location. (c) Analysis of the aster shape: the left y-axis shows the asphericity of the aster (defined as the ratio of eigenvalues of the shape tensor, see Section IV A) and the right y-axis shows the y-coordinate of the center of mass of the aster as a function of time. (d) The aster profile, as characterized by the polarization direction θ as a function of the azimuthal angle ϕ around the aster center (defined as the position at which the density is maximum, which coincides with the defect core, where τ = 0). The measurement is taken at radius r = 20 from the core. (e) The active torque, ∇ω × τ, integrated over the left, Ω_{L}, and right, Ω_{R}, subdomains of the aster (see Section IV A) as a function of time, showing the driving forces for unwinding and closing of the aster in each subdomain due to activity gradients. The objective function parameters are {A, B, C, D, K, ω_{0}} = {0.1, 1.0, 2.0, 2.0, 0, 0.05}, λ = 0.8, and the simulation box size is 128 × 128. A video of this trajectory is in ESI† Movie S1. |

To quantify how the position and profile of the aster change over the course of advection, we measure its center-of-mass position (x_{COM}, y_{COM}) and asphericity. Here, we track the y-coordinate of the center of mass, which is calculated as , and the asphericity is given by the ratio of eigenvalues of shape tensor: , where Latin indices denote grid points and Greek indices correspond to Cartesian coordinates, and r_{ij} is the distance of the ijth grid point from the center of mass. As shown in Fig. 2c, the control window naturally partitions into two stages. During the initial stage (0 < t ≲ 100) the aster rapidly changes shape into the comet-like configuration, as seen by the decrease in its asphericity, while simultaneously undergoing advection, moving towards the target point. Then, over the remaining long time window (100 ≲ t < t_{F}) the aster reforms slowly, the asphericity increases back to 1, and it moves the remaining small distance to the target position.

For further description of the aster profile during these stages, we present the angle of the polarization field θ as a function of the azimuthal angle ϕ around the aster core at four time points in Fig. 2d. At the initial time (t = 0) the system is radially symmetric with polarization vectors pointing toward the aster center, while by t = 50 and t = 100 the symmetry of θ in the top [0,π] and bottom [−π,0] quadrants is broken, with more polarization at the bottom points toward the advection direction, and the front-end remains aster-like with a radial configuration. The orientation returns to the aster configuration by t = t_{F}.

Further, we can understand the control solution physically by noting that the dynamics of τ is such that gradients in the control field ω create a torque on the orientation field, i.e., ∂_{t}τ ∼ −ρ∇ω or equivalently ∂_{t}θ ∼ ∇ω × τ. We then calculate the integral of the torque over two domains, Ω_{L}: [x ∈ {x_{0} − 10, x_{0}}, y ∈ {0, L_{y}}] on the left and and Ω_{R}: [x ∈ {x_{0}, x_{0} + 10}, y ∈ {0, L_{y}}] to the right, with x_{0} = 40 in our case (Fig. 2e). When the aster unwinds and advects, the region to the left has a positive torque (countercockwise rotation) and the region on the right experiences negative torque (clockwise rotation), which correspond to the partial unwinding of the aster. This is also illustrated by the snapshot shown for t = 35, where the dashed line (x = 40) represents the axis along the aster's motion. During the subsequent reformation phase, as the aster regains circular symmetry, the profile winds back such that the left/right subdomains experience clockwise/counterclockwise torque respectively. The snapshot shown at t = 1000 illustrates this behavior. Eventually, the system relaxes sufficiently close to an unperturbed aster configuration that the net torque becomes zero.

1. Specifying the trajectory of aster advection.
A limitation of the approach described thus far is that convergence of the control solution becomes unreliable when trying to advect the aster over distances significantly greater than its size (Δȳ ≈ 30). This is because the gradients in the objective function are extremely shallow for the initial stages of the trajectory when the target state is far from the initial state. While techniques to find global minima can potentially overcome this problem, an alternative approach is to change the objective function to ensure sufficient gradients at all stages. A simple example of the latter strategy is to prescribe the entire trajectory of the aster. This approach has the added benefit of controlling the translocation speed, but has the potential drawback of arriving at a suboptimal solution (either slower translocation or higher control cost), since the problem is more constrained.

We applied the latter strategy to the problem of translocating an aster a distance of Δȳ = 120. We formulated the control problem to translocate the aster at a constant speed for a time t_{adv}, followed by a time t_{reform} for reformation of the aster. We set the target trajectory as a series of configurations in which the initial state with an aster has been translated by 1 additional unit along the translocation axis. The configurations are separated by a time increment Δ = t_{adv}/Δȳ, so that in the final target state configuration at time t = t_{adv}, the aster center-of-mass will have traveled Δȳ, with a mean advection speed of . We then allow the additional time t_{reform} for the aster to acquire its target profile. We find that specifying the path in this way allows specifying a target distance that is arbitrarily far without any difficulties in achieving convergence of the control solution.

Fig. 3a and b show the system configurations and corresponding control solution ω for an example with = 1/20, t_{adv} = 2400, and t_{reform} = 1200. At t = 0, the activity is maximum at the core of the aster, but unlike the previous setup where only initial and final state of the aster are specified, the order of magnitude remains same throughout the advection phase of aster. During the first phase of the solution (constant advection), the aster undergoes partial dissolution and, as intended, a roughly steady rate of translation toward the target (see snapshots at t = 1200, 1900). However, because we specified the trajectory at discrete intervals spaced by Δȳ = 1, the optimal control field oscillates with a period of about Δ ≈ Δȳ/ = 20. These oscillations arise due to the discrete nature of the objective function that we have imposed, and thus they are not eliminated by penalizing time variations in activity with nonzero K.^{77} The oscillatory behavior is evident in Fig. 3c, which shows the positions of the maxima of ρ and ω and the minimum of the y-component of the torque, τ_{y}, as a function of time. The maximum in the control solution ω exhibits strong oscillations of 20 length units between the front and rear of the aster (while the activity remains low at the aster core, see Fig. 3b), whereas the density maximum moves at a nearly continuous speed toward the target. The minimum τ_{y} tracks polarization toward the −ŷ direction and it consistently coincides with the high activity point at the front of the aster. Taken together, these observations show that the control solution pushes the aster from the rear, while exerting torques at the front that maintain aster-like polarization. This is captured in Fig. 3d which shows θ as a function of the azimuthal angle ϕ at two intermediate times during the advection phase, t = 1200 and 1900. Finally, during the second (reformation) phase, the aster re-acquires its radially symmetric steady state density and polarization profile. The dynamical interplay among these forces can be seen in the ESI.†^{}

Fig. 3 Prescribing the path and speed of aster advection. The control problem is formulated in two stages: in advection, the aster moves at a mean speed of for t_{adv} time units; in reformation, the aster reacquires its steady-state profile over a timescale t_{reform}. The figure shows two examples. In both cases, the control task is to move the aster by 120 length units in the ŷ direction. (a)–(d) Example 1: slower advection, with = 1/20, t_{adv} = 2400, t_{reform} = 1200. (e)–(h) Example 2: faster advection, with = 1/5, t_{adv} = 600, t_{reform} = 1200. (a) and (e) Snapshots of the density (color map) and polarization (arrows) profiles for examples 1 (a) and 2 (e). For both examples, snapshots are shown for the initial state t = 0, two intermediate times during the advection phase, and the final point at t = 3600 (slow)/t = 1800 (fast). (b) and (f) Corresponding snapshots for the activity field (color map). The ‘x’ symbols in the snapshots in ((a), (b), (e) and (f)) show the position of the aster core. (c) and (g) Tracking the progress of the aster and the control solution. The plot shows the y-components of the position corresponding to the aster core (density maximum, ρ, green curve); activity maximum (ω, blue curve); and the minimum torque (τ_{y}, red curve). (d) and (h) The polarization direction θ as a function of the azimuthal angle ϕ, measured at a distance r = 20 from the aster core, at indicated times. The objective function parameters for both examples are the objective function parameters are {A, B, C, D, K, ω_{0}} = {0.1, 1.0, 2.0, 2.0, 0, 0.05}, λ = 0.8, and the simulation box size is 120 × 200. Videos corresponding to Examples 1 and 2 are provided in ESI† Movies S2 and S3. |

Since we are specifying the path of the aster, we can investigate how the control solution depends on the chosen advection rate. Fig. 3e–h shows analogous results for a trajectory in which the advection phase is shortened to t_{adv} = 600, forcing a higher translation speed = 1/5. The higher speed leads to a qualitatively different type of trajectory; the aster unwinds into a flock during the advection phase, and then reforms during the second phase. Here the control solution takes a bean-shaped spatial profile, which initially pulsates periodically to unwind the leading edge of the aster and push the aster toward the target position. At early times (by t = 150) the rear of the aster adopts a flock-like state with polarization primarily pointing in the ŷ-direction; by t = 450 most of the aster becomes flock-like. The extent of polarization along ŷ is particularly clear from the plot of θ(ϕ) at t = 450 (Fig. 3h, green triangles).

The trajectories in Fig. 3 demonstrate that imposing different constraints in the control problem, such as the time duration, can lead to very different activity profiles and intermediate states in the optimal solution. We can qualitatively, but not quantitatively, understand the transition from aster-like to flock-like behaviors to result from higher values of activity ω required to meet the imposed advection speed (due to the shorter time duration). In addition to setting the effective self-propulsion speed, the activity ω controls the effective compressibility of the active fluid, with small values of ω corresponding to negative compressibilities that favor asters and large values corresponding to positive compressibilities that favor polar flocks. These behaviors are evident in the phase diagram of the uncontrolled system (Fig. 1a). For the value of the particle interaction coefficient λ = 0.8, the uncontrolled system steady-state transitions from asters to ‘blobs’ to polar flocks at about ω = 0.2 and ω = 0.25. Thus, we might expect the system to at least locally undergo a transition out of the aster state above a threshold advection speed and corresponding local activity. However, due to the significant spatial variations of omega in the control solution, we cannot quantitatively connect the trajectory transition to the phase diagram. In particular, we observe maximum activity values (taken over space and time) of ω_{max} = 0.9 and 1.6 for the slow and fast advection trajectories respectively, meaning that the maximum local activity value is well above the phase transition threshold even for the aster-like trajectory.

Fig. 4 Changing the propagation direction of stripes. (a) and (b) Snapshots of (a) the density (color map) and polarization (arrows) profiles and (b) the activity field ω (color map). The system is initialized in an unperturbed stripe steady-state traveling in the + direction, with parameters {ω, λ} = {0.4, 0.0}. The control solution begins at t = 0 and the system state is shown at indicated times. (c) Contribution of each term of the density dynamics, (3), evaluated by integration across the defined sub-domain: Ω_{δ}: 0 < x ≤ 50 and 0 < y ≤ 100 (see Section IV B). The terms are: (density change), which is driven by −ω∇·τ (self-propulsion), −τ·∇ω (density flow due to activity gradients) and ∇^{2}ρ (diffusion due to density gradients). (d) The active torque ∇ω × τ integrated over the entire domain with time for different target orientations. The objective function parameters are {A, B, C, D, K, ω_{0}} = {0.1, 1.0, 7.0, 7.0, 0, 0.4}; λ = 0.0, and the simulation box size is 256 × 256. Videos S4 and S6 (ESI†) respectively show this trajectory and an independent run in which the stripe is forced to re-orient by 90°. |

While a straightforward route to reorienting a strip would be to melt the stripe to an isotropic domain and then have it reform in the new direction, this is not the optimal solution given by the control theory. Starting with the density equation, eqn (3), we investigate the primary forces influencing density evolution during the stripe reorientation process. For this we choose a subdomain 0 < x ≤ 50, 0 < y ≤ 100 within the simulation box of size 256 × 256. We integrate each of the three terms in eqn (3) over the specified subdomain as a function of time (Fig. 4c): −ω∇·τ, which governs the convection of active particles at convection speed ω; −τ·∇ω, which determines the local density dynamics due to gradients in activity; and ∇^{2}ρ, which determines the diffusion due to density gradients. We find that, at all times, the dominant contribution arises from self-propulsion, −ω∇·τ; contributions from gradients in activity and density have negligible contributions. Thus, we conclude that activity gradients are not the driving force for the density dynamics, but rather lead to the torques that reorient the polarization field, as described in our analysis in Section IV A. To quantify the effect of active torque in this case, we illustrate in Fig. 4c that as the difference in orientation between the initial state and target state increases, the active torque also increases, and as the system settles into the target orientation, the active torque goes to 0.

This result can be qualitatively understood as follows. A trajectory which proceeds through an isotropic intermediate state would require large magnitudes and gradients of activity to melt the stripe and then reform it in a new direction. Instead, the optimal solution takes advantage of the fact that a stripe arises through a spontaneous breaking of rotational symmetry, and thus requires only small torques (and correspondingly small orientational biases in gradients of activity) to reorient its direction.

Fig. 5 Remodelling a propagating stripe to a stationary aster. The system is initialized in an unperturbed stripe steady-state traveling in the + direction, with parameters {ω, λ} = {0.4, 0.6} seeded at (x_{0} = 60, y_{0} = 60). The target state is obtained from a simulation of an aster steady state with parameters {ω, λ} = {0.04, 0.6}. (a) and (b) Snapshots of (a) the density (color map) and polarization (arrows) profiles and (b) the activity field (color map). The final activity is equal to the activity of the steady state aster, ω = 0.04 everywhere except near the defect core, where it approaches zero. The baseline activity ω_{0} = 0.4 was set to the value corresponding to the stripe steady state. (c) The active torque ∇ω × τ integrated over time for each of the two subdomains 0 < x ≤ L_{x}, 0 < y ≤ L_{y}/2 and L_{y}/2 < y ≤ L_{y}. The objective function parameters are {A, B, C, D, K, ω_{0}} = {0.1, 1.0, 5.5, 5.5, 0, 0.4}; λ = 0.6; and the simulation box size is 128 × 128. A video of this trajectory is in ESI† Movie S5. |

We note that this control problem also differs from the previous one (reorienting stripes) in that the initial and final state have different symmetries. However, the optimal solution still avoids isotropic intermediate states, since these would require significantly larger magnitudes and gradients of activity.

(10) |

Fig. 6 Robustness of the control solution to adding noise to the initial condition. (a) Deviation as a function of time between the trajectories integrated computed with and without noise (, eqn (10)). We added Gaussian noise to the initial condition with indicated magnitude , and integrated the dynamics using the control protocol computed in the absence of noise that is shown in Fig. 2. The inset shows the deviation of the final state as a function of noise magnitude . (b) and (c) The initial and final states for (b) small and (c) large noise. We see robust aster structures even when the error is large in the end state for 35% noise. Parameter values are as in Fig. 2. |

To investigate the effect of model parameter errors on the optimal control solution, we performed trials for the aster advection problem shown in Fig. 2 with varying λ. That is, we used the same initial condition and the optimal protocol computed with λ = 0.8 (i.e. the activity sequence shown in Fig. 2b), but we integrated the equations of motion with different values of λ. As shown in Fig. 7, we see small errors for variations in λ smaller than 10%, beyond which the errors rise significantly. However, similar to the noise results, we observe well-formed asters even for cases with large variations, although their final positions deviate from the target state (Fig. 7c).

Fig. 7 Robustness of the control solution to variations in the model parameter λ. (a) The deviation between the trajectory integrated with the optimal control protocol from Fig. 2b but with the indicated value of λ, compared to the one with λ = 0.8 shown in Fig. 2b. (b) and (c) Snapshots showing the initial and final states for (b) relatively small (λ = 0.78) and (c) large (λ = 0.9) parameter variations. We observe well-formed aster structures even though the error is relatively large at the end state for λ = 0.9, although the aster position deviates from the target. Parameter values are as in Fig. 2. |

Further, we show that the optimal control solutions are robust to noise. In particular, perturbing the initial condition by up to 20% leads to minimal quantitative deviations from the target behavior, and the solution remains qualitatively accurate for significantly larger perturbations. Also, we note that additional strategies can be employed for experimental systems where larger noise sources or systematic errors are unavoidable. This includes integrating closed-loop control components. For example, one can observe the current state of the system at regular intervals along a trajectory, and recompute the optimal control solution using the current state as the initial condition. Alternatively, one can add linear feedback terms that analyze deviations from the pre-computed optimal control trajectory.^{78}

In addition to directly applying the computed activity protocols, examination of their forms provides both fundamental and practical insights into controlling active materials. In particular, we show how the spatial gradients in the applied activity field lead to localized torques which rotate polarization directions, leading to the programmed reformulation of the pattern of interest (e.g. aster or stripe). Unsurprisingly, the form of the trajectory is different depending on the task being encoded for – changing the broken symmetry state of the system (e.g. stripe-to-aster, Fig. 5) requires very different spatial arrangements of active torques then advection (Fig. 2 and 3) or reorientation (Fig. 4). Notably however, the applied activity field and corresponding trajectory also depend strongly on the time allowed for the transformation. In the example of advecting the aster over a distance many times its size (Fig. 3), the aster mostly retains its form throughout the course of the trajectory when moving at moderate speed, but when forced to complete the journey 4× faster, the applied activity field reshapes the aster into a localized flock or swarm, which reformulates into an aster upon reaching the target position.

The states we seek to control are realizable in experiments. Propagating concentration waves of aligned self-propelled particles have been observed in dense actin motility assays^{79–82} and in self-chemotactic bacterial systems.^{83–91} Asters are ubiquitous in cell biology in processes such as the formation of the mitotic spindle, oogenesis, and plant cell cytokinesis.^{92–99} They can be reliably obtained in in vitro suspensions of cytoskeletal filaments and motor proteins.^{6,46,95,100–110} In such systems, activity can be controlled in space and time by constructing active materials with optogenetic molecular motors. For example, the closely related microtubule-based active nematic system has been modified to contain light-activated kinesins clusters, so that the activity is proportional to the local light intensity.^{43–46} The researchers then used a digital light projector to shine programmed spatiotemporal sequences of light on the sample.^{43–46} Our optimal control protocol can be used to determine the spatiotemporal light sequence that will drive such a system into a desired state. Due to inevitable accuracies in hydrodynamic descriptions of active matter systems and experimental noise, it is important that we found our optimal control solutions to be robust against noise and parameter inaccuracies (Fig. 6 and 7). Further, optimal control solutions can be corrected by periodically observing the state of the system and applying feedback control. For example, feedback can be implemented by adding linear restoring forces that drive the system toward the computed trajectory, or by recomputing the optimal control solution based on the current state of the system.

The optimal control framework presented here is highly generalizable, and can be readily applied to any system provided there is a means to externally actuate the system and there is a reasonably accurate continuum model. Importantly, the control variable need not be limited to the activity field, since the objective function can be extended to include any property of the material that can be actuated. With the recent success of automated PDE learning tools in discovering continuum models for active systems (e.g. ref. 111–113), applications need not be limited to systems with accurate models already available. Furthermore, since model discovery tools work better when provided with a variety of data, including from non-steady-state observations, we anticipate that combining model discovery tools with optimal control could be a powerful approach to both discover more accurate models and enhance the reliability of the control solutions.

- K.-A. Liu and I. Lin, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 011924 CrossRef PubMed.
- A. Czirók and T. Vicsek, Phys. A, 2000, 281, 17–29 CrossRef.
- C. Pierce, H. Wijesinghe, E. Mumper, B. Lower, S. Lower and R. Sooryakumar, Phys. Rev. Lett., 2018, 121, 188001 CrossRef PubMed.
- D. Needleman and Z. Dogic, Nat. Rev. Mater., 2017, 2, 1–14 Search PubMed.
- G. Sarfati, A. Maitra, R. Voituriez, J.-C. Galas and A. Estevez-Torres, Soft Matter, 2022, 18, 3793–3800 RSC.
- F. Ndlec, T. Surrey, A. C. Maggs and S. Leibler, Nature, 1997, 389, 305–308 CrossRef.
- T. Surrey, F. Nédélec, S. Leibler and E. Karsenti, Science, 2001, 292, 1167–1171 CrossRef CAS PubMed.
- M. L. Gardel, K. E. Kasza, C. P. Brangwynne, J. Liu and D. A. Weitz, Methods Cell Biol., 2008, 89, 487–519 CAS.
- M. L. Gardel, I. C. Schneider, Y. Aratyn-Schaus and C. M. Waterman, Annu. Rev. Cell Dev. Biol., 2010, 26, 315–333 CrossRef CAS PubMed.
- M. Dogterom and G. H. Koenderink, Nat. Rev. Mol. Cell Biol., 2019, 20, 38–54 CrossRef CAS.
- M. Soarese Silva, M. Depken, B. Stuhrmann, M. Korsten, F. C. MacKintosh and G. H. Koenderink, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 9408–9413 CrossRef PubMed.
- C. G. Wagner, M. M. Norton, J. S. Park and P. Grover, Phys. Rev. Lett., 2022, 128, 028003 CrossRef CAS.
- W. Wang, W. Duan, S. Ahmed, A. Sen and T. E. Mallouk, Acc. Chem. Res., 2015, 48, 1938–1946 CrossRef CAS PubMed.
- J. R. Gomez-Solano, S. Samin, C. Lozano, P. RuedasBatuecas, R. van Roij and C. Bechinger, Sci. Rep., 2017, 7, 14891 CrossRef.
- O. Hallatschek, S. S. Datta, K. Drescher, J. Dunkel, J. Elgeti, B. Waclaw and N. S. Wingreen, Nat. Rev. Phys., 2023, 1–13 Search PubMed.
- S. Liu, S. Shankar, M. C. Marchetti and Y. Wu, Nature, 2021, 590, 80–84 CrossRef CAS PubMed.
- A. E. Hamby, D. K. Vig, S. Safonova and C. W. Wolgemuth, Sci. Adv., 2018, 4, eaau0125 CrossRef CAS PubMed.
- B. Ni, R. Colin, H. Link, R. G. Endres and V. Sourjik, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 595–601 CrossRef CAS.
- E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601 CrossRef.
- D. Gonzalez-Rodriguez, K. Guevorkian, S. Douezan and F. Brochard-Wyart, Science, 2012, 338, 910–917 CrossRef CAS PubMed.
- S. Henkes, K. Kostanjevec, J. M. Collinson, R. Sknepnek and E. Bertin, Nat. Commun., 2020, 11, 1405 CrossRef CAS PubMed.
- E. Fodor and M. C. Marchetti, Phys. A, 2018, 504, 106–120 CrossRef.
- J. W. Nichol and A. Khademhosseini, Soft Matter, 2009, 5, 1312–1319 RSC.
- G. Zhang, R. Mueller, A. Doostmohammadi and J. M. Yeomans, J. R. Soc., Interface, 2020, 17, 20200312 CrossRef.
- A. Doostmohammadi, S. P. Thampi and J. M. Yeomans, Phys. Rev. Lett., 2016, 117, 048102 CrossRef PubMed.
- T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, Nature, 2017, 544, 212–216 CrossRef PubMed.
- X. Trepat, M. R. Wasserman, T. E. Angelini, E. Millet, D. A. Weitz, J. P. Butler and J. J. Fredberg, Nat. Phys., 2009, 5, 426–430 Search PubMed.
- X. Trepat and E. Sahai, Nat. Phys., 2018, 14, 671–682 Search PubMed.
- G. Duclos, S. Garcia, H. Yevick and P. Silberzan, Soft Matter, 2014, 10, 2346–2353 RSC.
- C. Blanch-Mercader, V. Yashunsky, S. Garcia, G. Duclos, L. Giomi and P. Silberzan, Phys. Rev. Lett., 2018, 120, 208101 CrossRef PubMed.
- S. Garcia, E. Hannezo, J. Elgeti, J.-F. Joanny, P. Silberzan and N. S. Gov, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 15314–15319 CrossRef CAS PubMed.
- M. Ghibaudo, A. Saez, L. Trichet, A. Xayaphoummine, J. Browaeys, P. Silberzan, A. Buguin and B. Ladoux, Soft Matter, 2008, 4, 1836–1843 RSC.
- T. Vicsek and A. Zafeiris, Phys. Rep., 2012, 517, 71–140 CrossRef.
- H. Chaté, Annu. Rev. Condens. Matter Phys., 2020, 11, 189–212 CrossRef.
- A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans and F. Sagués, Nat. Commun., 2018, 9, 1–13 CrossRef CAS PubMed.
- C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006 CrossRef.
- G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U. B. Kaupp and L. Alvarez, et al. , J. Phys.: Condens. Matter, 2020, 32, 193001 CrossRef PubMed.
- M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143 CrossRef.
- A. Aubret, M. Youssef, S. Sacanna and J. Palacci, Nat. Phys., 2018, 14, 1114–1118 Search PubMed.
- J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine and P. M. Chaikin, Science, 2013, 339, 936–940 CrossRef PubMed.
- J. Arlt, V. A. Martinez, A. Dawson, T. Pilizota and W. C. Poon, Nat. Commun., 2018, 9, 768 CrossRef.
- G. Frangipane, D. Dell’Arciprete, S. Petracchini, C. Maggi, F. Saglimbeni, S. Bianchi, G. Vizsnyiczai, M. L. Bernardini and R. Di Leonardo, eLife, 2018, 7, e36608 CrossRef PubMed.
- L. M. Lemma, M. Varghese, T. D. Ross, M. Thomson, A. Baskaran and Z. Dogic, PNAS Nexus, 2023, 2, pgad130 CrossRef PubMed.
- Z. Zarei, J. Berezney, A. Hensley, L. Lemma, N. Senbil, Z. Dogic and S. Fraden, Soft Matter, 2023, 19, 6691–6699 RSC.
- R. Zhang, S. A. Redford, P. V. Ruijgrok, N. Kumar, A. Mozaffari, S. Zemsky, A. R. Dinner, V. Vitelli, Z. Bryant and M. L. Gardel, et al. , Nat. Mater., 2021, 20, 875–882 CrossRef CAS.
- T. D. Ross, H. J. Lee, Z. Qu, R. A. Banks, R. Phillips and M. Thomson, Nature, 2019, 572, 224–229 CrossRef CAS.
- S. Shankar, L. V. Scharrer, M. J. Bowick and M. C. Marchetti, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2400933121 CrossRef CAS PubMed.
- M. Nasiri, H. Löwen and B. Liebchen, Europhys. Lett., 2023, 142, 17001 CrossRef.
- M. Kneević, T. Welker and H. Stark, Sci. Rep., 2022, 12, 19437 CrossRef PubMed.
- S. L. Brunton and J. N. Kutz, Data-driven science and engineering: Machine learning, dynamical systems, and control, Cambridge University Press, 2022 Search PubMed.
- G. E. Dullerud and F. Paganini, A course in robust control theory: a convex approach, Springer Science & Business Media, 2013, vol. 36 Search PubMed.
- C. Villani, et al., Optimal transport: old and new, Springer, 2009, vol. 338 Search PubMed.
- M. M. Norton, P. Grover, M. F. Hagan and S. Fraden, Phys. Rev. Lett., 2020, 125, 178005 CrossRef CAS.
- C. Sinigaglia, F. Braghin and M. Serra, Phys. Rev. Lett., 2024, 132, 218302 CrossRef CAS PubMed.
- J. Toner and Y. Tu, Phys. Rev. Lett., 1995, 75, 4326 CrossRef CAS PubMed.
- J. Toner and Y. Tu, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 58, 4828 CrossRef CAS.
- H. Y. Lee and M. Kardar, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 64, 056113 CrossRef CAS PubMed.
- K. Husain and M. Rao, Phys. Rev. Lett., 2017, 118, 078104 CrossRef PubMed.
- D. Geyer, A. Morin and D. Bartolo, Nat. Mater., 2018, 17, 789–793 CrossRef CAS.
- V. M. Worlitzer, G. Ariel, A. Be’er, H. Stark, M. Bär and S. Heidenreich, New J. Phys., 2021, 23, 033012 CrossRef CAS.
- S. Mishra, A. Baskaran and M. C. Marchetti, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 061916 CrossRef.
- A. Gopinath, M. F. Hagan, M. C. Marchetti and A. Baskaran, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 061903 CrossRef PubMed.
- J.-B. Caussin, A. Solon, A. Peshkov, H. Chaté, T. Dauxois, J. Tailleur, V. Vitelli and D. Bartolo, Phys. Rev. Lett., 2014, 112, 148102 CrossRef.
- H. Reinken, S. H. L. Klapp, M. Bär and S. Heidenreich, Phys. Rev. E, 2018, 97, 022613 CrossRef CAS.
- W. Ngamsaad and S. Suantai, Phys. Rev. E, 2018, 98, 062618 CrossRef CAS.
- E. Bertin, M. Droz and G. Grégoire, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 022101 CrossRef.
- A. Peshkov, I. S. Aranson, E. Bertin, H. Chaté and F. Ginelli, Phys. Rev. Lett., 2012, 109, 268701 CrossRef.
- A. Baskaran and M. C. Marchetti, Phys. Rev. Lett., 2008, 101, 268101 CrossRef.
- P. Srivastava, R. Shlomovitz, N. S. Gov and M. Rao, Phys. Rev. Lett., 2013, 110, 168104 CrossRef PubMed.
- G. I. Menon, Rheology of complex Fluids, 2010, pp. 193–218 Search PubMed.
- A. P. Solon and J. Tailleur, Phys. Rev. Lett., 2013, 111, 078101 CrossRef PubMed.
- D. E. Kirk, Optimal control theory: an introduction, Courier Corporation, 2004 Search PubMed.
- S. Lenhart and J. T. Workman, Optimal control applied to biological models, Chapman and Hall/CRC, 2007 Search PubMed.
- R. Kerswell, C. C. Pringle and A. Willis, Rep. Prog. Phys., 2014, 77, 085901 CrossRef PubMed.
- A. Borzí and V. Schulz, Computational optimization of systems governed by partial differential equations, SIAM, 2011 Search PubMed.
- T. Dupont, J. Hoffman, C. Johnson, R. C. Kirby, M. G. Larson, A. Logg and L. R. Scott, The fenics project, Chalmers Finite Element Centre, Chalmers University of Technology, 2003 Search PubMed.
- At least up to K = 1, larger values led to convergence problems. We expect that interpolating the target state would lead to a smoother solution profile.
- J. Bechhoefer, Control Theory for Physicists, Cambridge University Press, 2021 Search PubMed.
- V. Schaller, C. Weber, C. Semmrich, E. Frey and A. R. Bausch, Nature, 2010, 467, 73–77 CrossRef CAS.
- L. Huber, R. Suzuki, T. Krüger, E. Frey and A. Bausch, Science, 2018, 361, 255–258 CrossRef CAS.
- A. Sciortino and A. R. Bausch, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2017047118 CrossRef CAS.
- R. Suzuki and A. R. Bausch, Nat. Commun., 2017, 8, 41 CrossRef.
- J. Saragosti, V. Calvez, N. Bournaveas, B. Perthame, A. Buguin and P. Silberzan, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 16235–16240 CrossRef PubMed.
- O. Pohl and H. Stark, Phys. Rev. Lett., 2014, 112, 238303 CrossRef PubMed.
- T. Bhattacharjee, D. B. Amchin, R. Alert, J. A. Ott and S. S. Datta, eLife, 2022, 11, e71226 CrossRef PubMed.
- R. Alert, A. Martínez-Calvo and S. S. Datta, Phys. Rev. Lett., 2022, 128, 148101 CrossRef PubMed.
- A. V. Narla, J. Cremer and T. Hwa, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2105138118 CrossRef CAS.
- J. Cremer, T. Honda, Y. Tang, J. Wong-Ng, M. Vergassola and T. Hwa, Nature, 2019, 575, 658–663 CrossRef CAS.
- T. Bhattacharjee, D. B. Amchin, J. A. Ott, F. Kratz and S. S. Datta, Biophys. J., 2021, 120, 3483–3497 CrossRef CAS PubMed.
- M. P. Brenner, L. S. Levitov and E. O. Budrene, Biophys. J., 1998, 74, 1677–1693 CrossRef CAS PubMed.
- C. Fei, S. Mao, J. Yan, R. Alert, H. A. Stone, B. L. Bassler, N. S. Wingreen and A. Komrlj, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 7622–7632 CrossRef CAS.
- T. U. Mayer, T. M. Kapoor, S. J. Haggarty, R. W. King, S. L. Schreiber and T. J. Mitchison, Science, 1999, 286, 971–974 CrossRef CAS PubMed.
- S. C. Schaffner and J. V. José, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 11166–11171 CrossRef.
- R. Loughlin, R. Heald and F. Nédélec, J. Cell Biol., 2010, 191, 1239–1249 CrossRef PubMed.
- M. Shirasu-Hiza, Z. E. Perlman, T. Wittmann, E. Karsenti and T. J. Mitchison, Curr. Biol., 2004, 14, 1941–1945 CrossRef PubMed.
- S. R. Heidemann and M. W. Kirschner, J. Cell Biol., 1975, 67, 105–117 CrossRef PubMed.
- S. R. Heidemann and M. W. Kirschner, J. Exp. Zool., 1978, 204, 431–443 CrossRef CAS PubMed.
- A.-C. Schmit, M. Vantard, J. De Mey and A.-M. Lambert, Plant Cell Rep., 1983, 2, 285–288 CrossRef CAS PubMed.
- F. Kumagai and S. Hasezawa, Plant Biol., 2001, 3, 4–16 CrossRef.
- B. Najma, W.-S. Wei, A. Baskaran, P. J. Foster and G. Duclos, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2300174121 CrossRef CAS.
- V. Wollrab, J. M. Belmonte, L. Baldauf, M. Leptin, F. Nédeléc and G. H. Koenderink, J. Cell Sci., 2019, 132, jcs219717 CrossRef CAS.
- P. J. Foster, S. Fürthauer, M. J. Shelley and D. J. Needleman, eLife, 2015, 4, e10837 CrossRef.
- A. Colin, P. Singaravelu, M. Théry, L. Blanchoin and Z. Gueroui, Curr. Biol., 2018, 28, 2647–2656 CrossRef CAS PubMed.
- W. Luo, C.-H. Yu, Z. Z. Lieu, J. Allard, A. Mogilner, M. P. Sheetz and A. D. Bershadsky, J. Cell Biol., 2013, 202, 1057–1073 CrossRef CAS PubMed.
- S. Stam, S. L. Freedman, S. Banerjee, K. L. Weirich, A. R. Dinner and M. L. Gardel, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E10037–E10045 CrossRef PubMed.
- J. Berezney, B. L. Goode, S. Fraden and Z. Dogic, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2115895119 CrossRef.
- T. Thoresen, M. Lenz and M. L. Gardel, Biophys. J., 2011, 100, 2698–2705 CrossRef PubMed.
- D. V. Köster, K. Husain, E. Iljazi, A. Bhat, P. Bieling, R. D. Mullins, M. Rao and S. Mayor, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, E1645–E1654 CrossRef PubMed.
- M. Glaser, J. Schnauß, T. Tschirner, B. S. Schmidt, M. Moebius-Winkler, J. A. Käs and D. M. Smith, New J. Phys., 2016, 18, 055001 CrossRef.
- F. Nédélec, T. Surrey and E. Karsenti, Curr. Opin. Cell Biol., 2003, 15, 118–124 CrossRef.
- C. Joshi, S. Ray, L. M. Lemma, M. Varghese, G. Sharp, Z. Dogic, A. Baskaran and M. F. Hagan, Phys. Rev. Lett., 2022, 129, 258001 CrossRef CAS PubMed.
- R. Supekar, B. Song, A. Hastewell, G. P. Choi, A. Mietke and J. Dunkel, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2206994120 CrossRef CAS.
- M. Golden, R. O. Grigoriev, J. Nambisan and A. Fernandez-Nieves, Sci. Adv., 2023, 9, eabq6120 CrossRef CAS PubMed.

## Footnote |

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sm00547c |

This journal is © The Royal Society of Chemistry 2024 |