# Optimal crawling: From mechanical to chemical actuation

P Recho, L Truskinovsky

► To cite this version:

P Recho, L Truskinovsky. Optimal crawling: From mechanical to chemical actuation. Physical Review E, 2024, 110, pp.054413. 10.1103/physreve.110.054413. hal-05399181

HAL Id: hal-05399181

https://hal.science/hal-05399181v1

Submitted on 4 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Optimal crawling: From mechanical to chemical actuation

P. Recho<sup>1,*</sup> and L. Truskinovsky<sup>2,†</sup>

<sup>1</sup>LIPhy, CNRS UMR 5588, Université Grenoble Alpes, F-38000 Grenoble, France

<sup>2</sup>PMMH, CNRS UMR 7636, ESPCI PSL, F-75005 Paris, France

<figure><img src="image_2.png" /></figure>

(Received 16 May 2024; revised 10 September 2024; accepted 16 October 2024; published 21 November 2024)

Taking inspiration from the crawling motion of biological cells on a substrate, we consider a physical model of self-propulsion where the spatiotemporal driving can involve both a *mechanical* actuation by active force couples and a *chemical* actuation through controlled mass turnover. When the material turnover is slow and the mechanical driving dominates, we find that the highest velocity at a given energetic cost is reached when the actuation takes the form of an active force configuration propagating as a traveling wave. As the rate of material turnover increases, and the chemical driving starts to dominate the mechanical one, such a peristalsis-type control progressively loses its efficacy, yielding to a standing-wave-type driving which involves an interplay between the mechanical and chemical actuation. Our analysis suggests a paradigm for the optimal design of crawling biomimetic robots where the conventional purely mechanical driving through distributed force actuators is complemented by a distributed chemical control of the material remodeling inside the force-transmitting machinery.

DOI: 10.1103/PhysRevE.110.054413

I. INTRODUCTION

In living systems the crawling mode of motility is ubiquitous and its thorough understanding on both mechanical and biochemical levels constitutes an important fundamental challenge [1–4]. The parallel problem of the design of soft robots that can efficiently crawl by themselves is an equally important engineering problem [5–10] which has been mostly studied from the perspective of a purely mechanical driving in the form of distributed actuators generating “active” force couples [11–18]. In this paper, taking inspiration from the importance of chemical processes in the crawling motion of biological cells, we consider the situation where, in addition to a mechanical actuation, the crawler material can turn over through a process driven by an out-of-equilibrium chemical reaction.

There exist two fundamentally different paradigms to represent the overdamped crawling motion of an active object on a solid substrate. One of them, which builds on older swimming theories [19–21], assumes that the shape of a deformable object is dynamically actuated to harvest friction forces with its environment. The motion is supported by the breaking of the time-reversal symmetry of each periodic stroke [22]. As initially formulated by Purcell in the context of swimming at a low Reynolds number, this necessary condition to obtain a net motion over a stroke is known as the scallop theorem [23–25]. The optimality of such motion as a function of the actuation can be quantified by the Stokes efficiency, which compares how much energy is dissipated during each stroke with the power to move the crawler center of mass at a certain velocity against the frictional background [26–29]. Limbless

animals such as millipedes, caterpillars, or earthworms are good natural representatives of this category of crawlers.

The development of active gel theories allowing one to model the motility of various types of biological cells [30] has paved the way for an alternative paradigm where the crawler is constituted of a mechanical skeleton (scaffold) which can chemically turn over through its polymerization and depolymerization. The fundamental point here is that while the skeleton is in frictional contact with the external environment, its depolymerized building blocks are not, in the same manner as the lower part of the belt of a rolling tank tread transmits tractions to the ground while the upper part of the belt does not. The mechanical actuation of the skeleton by distributed contractile forces then creates a steady-state flow [31]. The spatial analog of the scallop theorem would then state that to ensure self-propulsion the flow of skeleton must be asymmetric with respect to the crawler center. As a result, the friction forces on the background will also lose their symmetry, leading to directional motion [32,33]. The fact that the skeleton is being permanently advected in a particular direction requires to “close the stroke.” This means that there must be sources and sinks allowing material renewal and depletion where necessary, through an implicit chemical reaction [34]. Again, the Stokes efficiency for such type of motion has been introduced and the corresponding optimal regimes were identified [35]. While the implied modeling approach is most simply formulated in the context of an object crawling on a stiff substrate, it was also shown to be fully capable of explaining various regimes of swimming in a liquid [36].

The goal of the present paper is to bring together these two paradigms within a single simple prototypical framework. Specifically, we consider a one-dimensional (1D) soft elastic body with free boundaries, whose material can be chemically driven to undergo internal mass redistribution. In addition,

<footnote>*Contact author: pierre.recho@univ-grenoble-alpes.fr</footnote>

<footnote>†Contact author: lev.truskinovsky@espci.fr</footnote>
it can be mechanically driven through a distributed field of active stresses. Such a system may then be actuated both *mechanically* and *chemically*. It can be seen as representing in a very simplified way the propulsion machinery of a cell. It is then implied that the active mechanical stress is controlled by a distribution of internally driven molecular motors cross-linking skeleton filaments, while the active mass redistribution is governed by a chemical potential representing an internal out-of-equilibrium chemical reservoir of monomers that can be polymerized into filament form. In this respect, it is worth mentioning that the chemical pathways regulating both the cytoskeleton contractility and its turnover can now be externally controlled with light in both space and time using two different optogenetics constructions [37,38].

One of our main results is the realization that the interplay between the chemical and mechanical modalities of active control crucially depends on the relative rate of the corresponding kinetic processes. We characterize the relative importance of chemical versus mechanical activity by a single dimensionless parameter $\lambda$, and study its role in the choice of the optimal crawling gait. More precisely, we ask how the optimal actuation strategy adjusts as $\lambda$ increases from zero to infinity and whether one can identify transitions between different crawling gaits triggered by the variation of $\lambda$.

In particular, we show that in order to ensure optimal efficacy, our two modalities of external driving, mechanical and chemical, must cooperate. In other words, the two controlling agents must conspire if the goal is to achieve the best performance at a fixed total metabolic cost. How the corresponding mechanical and biochemical pathways are organized to reach the necessary level of coherency is beyond the scope of this work, where we neglect the chemomechanical feedbacks between the active agents and the dynamical variables representing, for instance, the flow or the density of material.

To formulate the optimal control problem we assume that our driving mechanisms arrive with an energetic cost. We show how the latter can be specified based on some basic thermodynamic arguments. We then use a simple close-to-equilibrium Onsager formalism to introduce the corresponding kinetic processes. This allows us to specify the two timescales characterizing the mechanical and chemical drivings and formally define the parameter $\lambda$ as a ratio of these timescales.

As the most basic assumption, we associate the performance of the propulsion with the average velocity of the crawler and solve in two limiting cases the mathematical optimization problem delivering the best actuation strategy at a fixed energetic cost. In agreement with previous investigations [14,39,40], we find that, in the case of slow material turnover, the most effective time-periodic actuation strategy is a *traveling wave* propagating along the body. In the opposite limit, when turnover is fast and both mechanical and chemical actuators are relevant, the optimal driving is represented by a *standing wave*. In the latter case, the combined chemomechanical actuation and the associated cooperation between chemistry and mechanics in driving the internal dynamics allow the system to reach a higher performance than in the absence of turnover at the same energetic cost.

Using a simple analytical expression linking the two limiting cases, where either mechanics or chemistry dominates, we study the crossover between the *traveling-wave-type* and the *standing-wave-type* modalities of crawling as the parameter $\lambda$ varies continuously. We show that such a crossover is accompanied by a switch in the allocation of the energetic resources from the purely mechanical driving towards an optimal balance between mechanical and chemical energy inputs. More specifically, we show that, as the rate of material turnover increases, and the chemical driving starts to dominate the mechanical one, the commonly accepted peristalsis-type control progressively loses its efficacy, yielding to a *standing-wave-type* driving which now involves a constructive interplay between the mechanical and chemical actuation.

In addition to offering insights regarding the different fundamental modalities of the functioning of living matter, our analysis can be viewed as providing an alternative paradigm for the design of soft robots that crawl on rigid surfaces. The main difference is in the replacement of the conventional purely mechanical driving through distributed force actuators by a chemomechanical driving, involving the possibility to chemically activate or deactivate the dynamic renewal of the internal force-transmitting scaffold.

The paper is organized as follows. In Sec. II we derive a thermodynamically consistent chemomechanical model of a driven one-dimensional crawling segment with free boundaries. In Sec. III, we introduce the main nondimensional parameter. The time-periodic chemomechanical actuation fields that control the dynamics are introduced in Sec. IV. In Sec. V, based on the energy balance discussed in Sec. II, we formulate the optimization problem whose solution allows one to specify the optimal chemomechanical control. Such a general problem is formulated in the weak actuation regime in Sec. VI, where the controls are considered to trigger only small deformations of the material. We then study two limiting cases where the optimal actuation protocol can be found explicitly. Specifically, in Sec. VII we assume that material turnover is slow compared to mechanical relaxation while in Sec. VIII we consider the opposite limit. Finally, we introduce in Sec. IX a specific form for the actuation fields that interpolates the two above limits. It enables us to study the interplay between the chemical and mechanical actuation when the relative rate of turnover varies continuously and to identify a transition between the two major gaits of chemomechanical crawling. Our conclusions are summarized in Sec. X.

## II. THE MODEL

We consider a prototypical model of a biomimetic object crawling along a one-dimensional track on a rigid substrate. The skeleton of this one-dimensional crawler is effectively represented as a continuum segment with material points indexed by the actual spatial coordinate $x \in [l_{-}(t), l_{+}(t)]$. Here the two moving boundaries representing the front and rear edges of the crawler are $l_{-}(t)$ and $l_{+}(t)$ with $t \ge 0$ denoting the time. The deformed state of the system is described by the time-dependent mapping $x(X, t)$, where $X$ we denote the reference positions of the actual points.
In view of the anticipated role of the chemically driven mass turnover of the skeleton, we describe the mechanical response of the system using the framework of morphoelasticity [41,42]. We start by introducing the deformation gradient $F(x, t) = \partial x/\partial X$ and decompose $F$ multiplicatively into

$$F = AG.$$

Here $A(x, t)$ is the 1D analog of the elastic distortion which can be considered in this setting as purely volumetric. We suppose that such a distortion is counted from an unstressed configuration with a fixed density $\rho_m$. Then the actual (current) density of the material is $\rho(x, t) = \rho_m/A(x, t)$. The scalar function $G(x, t)$ describes the stress-free swelling due to the arrival or departure of the material “building blocks” needed to assemble the configuration with the density $\rho_m$. Such a configuration can be interpreted as intermediate because the primordial configuration of the “building blocks” would be characterized by another density $\rho_0(x, t) = G(x, t)\rho_m$. The function $G(x, t)$ is a 1D analog of the growth tensor in morphoelasticity describing the swelling of the material from density $\rho_0(x, t)$ to the density $\rho_m$ while the scalar $A$ is the deformation gradient necessary to accommodate such a swelling into the actual stressed configuration. The thermodynamical configuration of the resulting material capable of both elastic deformation and inelastic renewal is described by two mechanical (dynamic) variables which may be either $A(x, t)$ and $G(x, t)$ or $\rho(x, t)$ and $\rho_0(x, t)$.

**Mass balance.** In the Eulerian coordinate system, the equation of mass balance can be rewritten in the form

$$\partial_t \rho + \partial_x (\rho v) = r, \quad (1)$$

where $v = \partial_t x(X, t)$ is the velocity of material points, and the mass supply is

$$r = \rho \left( \frac{1}{G} \frac{dG}{dt} \right),$$

where $d/dt = \partial_t + v\partial_x$ denotes the total derivative. The derivation of Eq. (1) is based on the mass conservation between the unstressed configuration and the current one: $\rho = \rho_m/A = \rho_m G/F$. Taking the total differential of this relation, we obtain $d\rho/dt = \rho_m(dG/dt/F - G/F^2dF/dt)$, so that

$$\frac{d\rho}{dt} + \rho_m \frac{G}{F^2} \frac{dF}{dt} = \frac{\rho}{G} \frac{dG}{dt}.$$

Since by a time differentiation of the deformation gradient $\partial_x v = (dF/dt)/F$, we obtain Eq. (1). See, e.g., Ref. [43].

It is clear from Eq. (1) that it is the temporal variation of $G$ which brings the local sources and sinks of mass. For instance, in the context of living cells, this term represents the polymerization and depolymerization of the cytoskeleton filaments, which is controlled by a monomer reservoir [44]. As we do not consider any flux of mass through the boundaries, we will also assume that

$$\dot{l}_{\pm} = v(l_{\pm}(t), t).$$

**Momentum balance.** Since inertia is negligible in our setting, we can write force balance in the form

$$\partial_x \sigma = f_f, \quad (2)$$

where $\sigma(x, t)$ is the axial stress and $f_f$ is the bulk force describing the interaction of the skeleton with the rigid substrate in a thin-film limit. We further assume that

$$\sigma(x, t) = \sigma_e(x, t) + \sigma_a(x, t), \quad (3)$$

where $\sigma_e(x, t)$ is the elastic stress and $\sigma_a(x, t)$ is an active stress, representing the mechanical driving and effectively describing the momentum exchange with an out-of-equilibrium reservoir. One can also think directly in terms of the driving bulk force and the elastic restoring force:

$$f_a(x, t) = \partial_x \sigma_a \quad \text{and} \quad f_e = -\partial_x \sigma_e.$$

The former can be exogenous (say, originating from external actuators such as an applied magnetic field affecting embedded beads) or endogenous (say, describing myosin molecular motors cross-linking actin filaments). The boundary condition associated with Eq. (2) is

$$\sigma(l_{-}(t), t) = \sigma(l_{+}(t), t),$$

which implies the existence of a stiff spring representing a constraint which connects the two edges of the crawling system $l_+$, $l_-$ and ensures that its total length

$$L = l_{+}(t) - l_{-}(t)$$

remains constant.

**Energy balance.** We consider that the system is isothermal and introduce its free energy in the form

$$\Psi = \int_{l_{-}}^{l_{+}} \rho \psi(A, G) dx, \quad (4)$$

where we assume that elasticity and turnover are uncoupled material properties of the skeleton:

$$\psi(A, G) = \psi_e(A) + \psi_c(G).$$

For simplicity, we consider as variable only the mechanical contribution $\psi_e(A)$ while the chemical contribution $\psi_c(G) = \tilde{\psi}$ is assumed constant. This reflects the fact that only variations of the elastic strain change the stored internal energy while the amount of stress-free swelling due to material turnover does not impact the chemical free energy as it is often formulated in classical morphoelastic theories [45]. In view of the local mass conservation (1), we obtain

$$\frac{d\Psi}{dt} = \int_{l_{-}}^{l_{+}} \left( r\psi + \rho \frac{d\psi}{dt} \right) dx = \int_{l_{-}}^{l_{+}} (\sigma_e \partial_x v + r\mu_c) dx,$$

where we introduced the elastic stress in the skeleton,

$$\sigma_e = \rho_m \frac{\partial \psi_e}{\partial A},$$

and the chemical potential of the skeleton,

$$\mu_c = \psi - A \frac{\partial \psi_e}{\partial A}.$$

Next we need to introduce the exerted power $\Pi$. Since we assumed for simplicity that there are no boundary traction forces (no cargo to carry), such a power is delivered only in the bulk, through both active stresses and active mass exchange.
The power of active forces acting on the system has the standard form of mechanical work per unit time,

$$\Pi^m = \int_{l_{-}}^{l_{+}} f_a v dx.$$

Similarly, the power of generalized forces performing the regulation of mass exchange and acting at the chemical level can be written as

$$\Pi^c = \int_{l_{-}}^{l_{+}} \mu_a r dx.$$

The external mass reservoir is characterized by the chemical potential $\mu_a(x, t)$, which is a direct chemical analog of $f_a(x, t)$, while the rate of the exchange reaction $r(x, t)$ is the analog of the local velocity $v(x, t)$.

Due to the fact that our system is active, the integral energy inequality representing the second law of thermodynamics can be interpreted in two equivalent forms.

First, considering $f_a$ and $\mu_a$ as internal drivings, we can write the integral energy balance as $-d\Psi/dt = \hat{R}$, where the function $\hat{R} = R - (\Pi^m + \Pi^c)$ describes all energy exchanges due to the interaction between the material system described by the energy [Eq. (4)] and the machinery operating to maintain the active terms. The latter includes the rate of supply of mechanical and chemical energy due to the microscopic active agents $\Pi^m + \Pi^c$ and the irreversible dissipation $R \ge 0$ characterizing the interaction of the system with a thermal reservoir. Note that if the active agents are designed to always operate in the regime where $\Pi^m + \Pi^c \ge 0$, as will be the case in the rest of the paper, the corresponding entry in $\hat{R}$ can be considered as representing “antidissipation.”

The second approach, which we effectively use in what follows, is based on the assumption that the power $\Pi^m + \Pi^c$ represents an external work. We can then write the integral energy balance in the form

$$\begin{aligned} R &= \Pi^m + \Pi^c - \frac{d\Psi}{dt} \\ &= \int_{l_{-}}^{l_{+}} [(f_a - f_e)v + (\mu_a - \mu_c)r]. \end{aligned} \quad (5)$$

Here the dissipation,

$$R \ge 0,$$

is interpreted as the difference between the work done per unit time by the macroscopic active agents $\Pi^m + \Pi^c$ and the concurrent rate of change of the macroscopic free energy of the system, $d\Psi/dt$ (see, e.g., Refs. [35,46]).

To ensure that the term $R$ representing the dissipated energy is non-negative as required by the second law of thermodynamics, we make the simple assumption that the linear Onsager close-to-equilibrium theory [47] is operative. Specifically, we assume that the thermodynamic fluxes and forces are linearly related so that

$$\begin{aligned} v &= \Delta f / \xi + \lambda_{12} \Delta \mu, \\ r &= \lambda_{21} \Delta f + \Delta \mu / v, \end{aligned} \quad (6)$$

where $\Delta f = f_a - f_e$ and $\Delta \mu = \mu_a - \mu_c$. Note that in a more general tensorial framework, the cross term should involve a

director field corresponding to the preferred direction of the material turnover. In the rest of the paper, we will neglect this chemomechanical cross term as it is usually done in the context of cell motility [30,49–57]. While it is in fact a possibility in view of the presence of various interconnecting pathways, to our knowledge, a direct material coupling between the dynamics of the cell adhesion and that of cytoskeleton turnover has not been so far well documented [58].

In Eqs. (6), we have introduced the purely mechanical friction coefficient $\xi > 0$ and the purely chemical kinetic coefficient $v > 0$. Under these assumptions, friction with the substrate, mimicking the interaction with the mechanical reservoir, and kinetic turnover of the material, mimicking the interaction with the chemical reservoir, represent two sources of the quadratic dissipation in the system. As we show below, the two parameters, $\xi$ and $v$, set the timescales of the corresponding mechanical and chemical relaxation processes and their dimensionless ratio is an important control parameter of the problem.

In Eqs. (6), we have also implicitly assumed only weak deviation from equilibrium, which implies a weak level of activation. This is a common assumption in similar conceptual studies (see, for instance, the assumptions behind the theory of active gels [59] and some elementary models of molecular motors [60]). While such choice mainly aims at achieving analytical transparency, it is also motivated by the fact that it is still not clear how strong activation leading to considerable deviation from equilibrium should be modeled.

**Specification of the model.** We further impose for simplicity that the free energy represents only entropic elasticity and set

$$\psi_e(A) = \psi_0(A \log(A) - A), \quad (7)$$

where $\psi_0$ characterizes the elastic stiffness of the material. The chosen expression (7) can be viewed as reflecting a skeleton containing almost ideal polymer chains. Note that $\partial_A \psi(1) = 0$, which agrees with our assumption that after the mass exchange the reference configuration with density $\rho_m$ remains stress-free. Under assumption (7), the elastic stress-strain relation takes the form

$$\sigma_e = -E \log \left( \frac{\rho}{\rho_m} \right), \quad (8)$$

where $E = \rho_m \psi_0$ is the Young modulus at the stress-free state. Expression (8) penalizes both infinite polymer chains compression and extension. The resulting Eq. (6)$_1$ can be written as

$$f_a - E \frac{\partial_x \rho}{\rho} = \xi v. \quad (9)$$

Analogously we can now write an explicit expression for the thermodynamic chemical potential

$$\mu_c = \tilde{\psi} + \frac{E}{\rho}, \quad (10)$$

which allows us to rewrite the resulting Eq. (6)$_2$ in the form

$$\mu_a - \tilde{\psi} + \frac{E}{\rho} = vr. \quad (11)$$

Equation (9) effectively describes the kinetics of the mechanical relaxation of the internal flow $v$ while Eq. (11)
describes the kinetics of the chemical relaxation of the material turnover; the two equations are coupled through the mass balance equation (1).

### III. NONDIMENSIONALIZATION

To analyze the resulting system of equations, it is convenient to introduce the change of variable:

$$y = \frac{x - S(t)}{L} \in [-\frac{1}{2}, \frac{1}{2}], \quad \text{where } S(t) = \frac{l_{-}(t) + l_{+}(t)}{2}$$

is the geometric center of the segment. We use L to nondimensionalize space, $\xi L^2/E$ to nondimensionalize time, $\rho_m$ to nondimensionalize density, and $E/L$ to nondimensionalize force. Note that to normalize the density we could have also used the natural scale $\xi^2 L^2/E$. In this case our system acquires an additional dimensionless parameter $K = (E\rho_m)/(\xi^2 L^2)$, which, however, can be completely eliminated by simply renormalizing the density $\rho$ to $\rho/K$. The parameter K would, of course, have remained essential if we had taken inertia into account.

Combining Eqs. (1), (9), and (11), the problem reduces to a single dimensionless reaction-drift-diffusion equation for the skeleton density field $\rho(x, t)$:

$$\partial_t \rho + \partial_y(\rho(f - V) - \partial_y \rho) = \frac{1}{\tau} \left( \frac{1}{\rho} - \frac{1}{\tilde{\rho}} + \mu \right), \quad (12)$$

with the boundary conditions

$$\rho|_{-1/2} = \rho|_{1/2} \quad \text{and} \quad V(t) = \left(f - \frac{\partial_y \rho}{\rho}\right)|_{\pm 1/2}. \quad (13)$$

In Eqs. (12) and (13),

$$V = \dot{I}_{\pm}$$

is the macroscopic velocity of the crawler, while $f = f_a/(E/L)$ and $\mu = \mu_a/(E/\rho_m)$ are dimensionless active spatiotemporal controls.

We also introduced two nondimensional parameters. The first one

$$\tau = \frac{\nu \rho_m^2}{\xi L^2}, \quad (14)$$

represents the ratio of the two characteristic timescales: the one characterizing the chemical reaction controlling the skeleton turnover and the other one describing mechanical relaxation due to the sliding friction against the rigid background. The second dimensionless parameter compares the stored elastic energy compared to the chemical stored energy:

$$\tilde{\rho} = \frac{E}{\rho_m \tilde{\psi}}. \quad (15)$$

In Sec. VI, we show that, as we will be interested in small material density variations only, the only important dimensionless parameter, controlling the choice of the optimal crawling strategy, will be the combination of the two parameters (14) and (15):

$$\lambda = \frac{1}{\tau \tilde{\rho}^2} = \frac{\xi \tilde{\psi}^2 L^2}{\nu E^2}. \quad (16)$$

Thus, when $\lambda = 0$, material turnover is absent and sliding friction with the substrate is the only source of dissipation and therefore the main rate-limiting process. Instead, when $\lambda \gg 1$, turnover is much faster than frictional relaxation and the energy entering the system is predominantly dissipating due to the chemical reaction.

### IV. ACTUATION

The free boundary problem (12), with the boundary conditions (13), also contains two dimensionless active fields, which are still not specified. One is the active force distribution $f(y, t)$ and the other one is the active chemical potential distribution $\mu(y, t)$. Both of them represent nonequilibrium reservoirs driving the system, a mechanical one and a chemical one, respectively.

Given that intrinsic mechanical action is exerted by force couples, it can be convenient to write

$$f = -\partial_y s,$$

where the potential $s(y, t)$ is the active stress [48]. In the context of ATP-driven actomyosin systems, it describes the field of force dipoles generating mechanical contraction. Instead, the field $\mu(y, t)$ represents the dynamic target of an out-of-equilibrium chemical reaction responsible for the internal turnover of matter.

Observe that when both fields vanish and, therefore, the driving is absent, the solution of Eqs. (12) and (13) is $\rho \equiv \tilde{\rho}$ and $V \equiv 0$. Two nonequilibrium limiting cases are of interest: $\mu = 0$, corresponding to a purely mechanical driving, and $s = 0$, when the actuation is purely chemical.

One possibility to set $s(y, t)$ and $\mu(y, t)$ is to impose that they are controlled by the presence of chemomechanical feedbacks so that the driving “follows” the responding system as the two are described by a coupled system of equations. For instance, in some models of contraction-driven cell motility, the actuation process is coupled with the dynamics of molecular motors which, in turn, is linked to the flow of cytoskeleton [35]. Other types of coupling involving a feedback relating material flow with the active stress describing various regulatory pathways have been considered, for instance, in Refs. [61–63]. Physically motivated systems involving the coupling between active turnover and an intracellular chemical messenger affected by the material flow have been considered as well [56].

### V. OPTIMAL DRIVING

In this paper we are not assuming any of the aforementioned feedbacks and instead ask the question of how to reverse engineer the crawling machinery to reach an optimal performance. One can think of our system as describing a soft crawling body driven either by internal (cell motility) or external (robotics) agents. The task is to understand how an optimal performance can be achieved when it can be driven using both mechanical—$f(y, t)$—and chemical—$\mu(y, t)$—actuation. We leave the question of whether such controls can be indeed implemented in a self-consistent and physically meaningful manner to a separate study.
Some works have previously addressed the issue of cell motility in the framework of control theory by resorting to chemomechanical controls in specific situations [33,51] but without considering the energetic cost of such actuation. The present paper complements and further develops this previous work.

The first step is to restrict the class of admissible controls $f(y, t)$ and $\mu(y, t)$. To reflect a cyclic nature of the actuation, we assume that these controls are $T$-periodic in time where $T$ is a positive number to be found in the process of finding the solution of the optimization problem. Since we do not take into account any resultant part in either mechanical and chemical drive, it is natural to assume that the time and space averages of the controls $f(y, t)$ and $\mu(y, t)$ are equal to zero:

$$ \langle f \rangle = \langle \mu \rangle = \bar{f} = \bar{\mu} = 0. \qquad (17) $$

Here

$$ \langle (\cdot) \rangle = \lim_{T \to \infty} \frac{1}{T} \int_{0}^{T} (\cdot) dt \text{ and } \overline{(\cdot)} = \int_{-1/2}^{1/2} (\cdot) dy $$

denote the time averaging and the space averaging, respectively.

First recall that the total time-averaged power injected into the system as a result of both mechanical and chemical activity is

$$ C = \overline{\langle \Pi^m + \Pi^c \rangle} = \overline{\langle fv + \mu r \rangle}. \qquad (18) $$

By averaging Eq. (5) over time, we obtain that $C \ge 0$, showing that the total actuation cost, expressed in this way, is always non-negative. This also shows that the power of the active driving is indeed antidissipative over each stroke. Substituting the expressions for $v$ and $r$ from the Onsager constitutive relations (9)–(11) into Eq. (18), we can rewrite the expression for the measure of the energetic cost in the form

$$ C = \left\langle f \left( f - \frac{\partial_y \rho}{\rho} \right) \right\rangle + \frac{1}{\tau} \left\langle \mu \left( \mu + \frac{1}{\rho} - \frac{1}{\tilde{\rho}} \right) \right\rangle. $$

The definition of the functional performance is less straightforward since the system can move even in the absence of cargo. Naturally, in this case one would like to associate a nonzero functionality even if the resultant applied force is equal to zero. To circumvent this classical problem, several proposals have been made. Among them we consider the Stokes performance [64]. In our case this means choosing the rate of frictional dissipation $\xi L \langle V \rangle^2$ necessary to advance the system as a rigid object as the functionality measure. When both the length $L$ of the system and the sliding friction coefficient $\xi$ are fixed, this choice is equivalent to associating the performance of the crawling mechanism with the averaged velocity:

$$ P = \langle V \rangle. $$

Other choices are possible as well, accounting, for instance, for the metabolic expenses required to maintain at a certain level the distributed active stresses (see, for instance, the discussion in Ref. [65]).

In this paper we have chosen to set the problem of maximizing the performance $P$ at a fixed energetic cost $C$. This optimal control problem, also considered in Refs. [14,40],

reflects the desire to compare crawler designs with an imposed availability of metabolic resources and to select the design that achieves the best performance. Other choices would be possible such as minimizing the cost at a fixed performance, which would correspond to selecting, among the crawlers that achieve a given performance, the design that consumes the less resources. It is also possible to maximize the efficiency $P/C$ (see Ref. [35]), which corresponds to an optimal trade-off between the cost and the performance. The choice between these and other seemingly arbitrary options is ultimately dictated by the targeted functionality of the crawler.

## VI. WEAK ACTUATION

To get analytic results, we only consider the case where the spatial and temporal inhomogeneities of the driving forces are small. Then if $\epsilon$ is a small parameter characterizing the scale of activation, we can write

$$ \begin{aligned} s(y, t) &= \epsilon s_1(y, t) + \epsilon^2 s_2(y, t) + \dots, \\ \mu(y, t) &= \epsilon \mu_1(y, t) + \epsilon^2 \mu_2(y, t) + \dots. \end{aligned} \qquad (19) $$

The smallness of the parameter $\epsilon$ reflects the fact that the actuation is assumed to be sufficiently weak to produce only small deformations (i.e., density variations) of the material. To quantify the precise range of values of $\epsilon$ where the linear regime is operative, one needs to solve numerically the full nonlinear optimization problem, assuming that the actuation fields can be large. This considerable task is beyond the scope of our exploratory paper and is left for a separate study.

To assess the cost $C$ and the performance $P$, we will need to compute the first two terms in both expansions. To this end, it will be convenient to introduce the new auxiliary variables

$$ u = \mu \tilde{\rho}, $$

characterizing chemical control, and

$$ w = s + \log(\rho/\tilde{\rho}), $$

mixing the mechanical control with the mechanical response. We can then also write the associated expansions:

$$ \begin{aligned} u(y, t) &= \epsilon u_1(y, t) + \epsilon^2 u_2(y, t) + \dots, \\ w(y, t) &= \epsilon w_1(y, t) + \epsilon^2 w_2(y, t) + \dots. \end{aligned} \qquad (20) $$

Inserting Eqs. (20) in Eqs. (12) and (13), the first-order problem takes the form

$$ \partial_t w_1 - \partial_{yy} w_1 + \lambda w_1 = \partial_t s_1 + \lambda (s_1 + u_1) \qquad (21) $$

with the periodic boundary conditions

$$ w_1|_{-1/2} = w_1|_{1/2} \text{ and } V_1(t) = -\partial_y w_1|_{\pm 1/2}. \qquad (22) $$

In the context of weak actuation, our choice of the auxiliary variables $u$ and $w$ leads to the appearance in Eq. (21) of only one nondimensional parameter $\lambda$ introduced in Eq. (16). Since $\langle u_1 \rangle = \langle s_1 \rangle = 0$, the first nonvanishing contribution to the cost is of second order, $C = \epsilon^2 C_2 + \dots$, where

$$ C_2 = \langle \partial_y w_1 \partial_y s_1 \rangle + \lambda \langle u_1 (s_1 + u_1 - w_1) \rangle. \qquad (23) $$

As it involves only first-order terms, $C_2$ can be computed from the solution of Eqs. (21) and (22).
Similarly, since $\langle V_1 \rangle = 0$, the first nonvanishing contribution to the performance is also second order, $P = \epsilon^2 P_2 + \dots$, where

$$P_2 = \langle V_2 \rangle. \qquad (24)$$

To compute this contribution we need to consider the second-order expansion for $w_2$:

$$\begin{aligned} \partial_t w_2 - \partial_{yy} w_2 + \lambda w_2 = & \partial_t s_2 + \lambda(s_2 + u_2) \\ & + (\partial_y w_1 + V_1)(\partial_y w_1 - \partial_y s_1) \\ & + \lambda(s_1 - w_1)(2u_1 - 3w_1 + 3s_1)/2, \end{aligned} \qquad (25)$$

with the periodic boundary conditions

$$w_2|_{-1/2} = w_2|_{1/2} \text{ and } V_2(t) = -\partial_y w_2|_{\pm 1/2}.$$

Then, after performing the temporal averaging of these equations over time, we obtain

$$\begin{aligned} -\partial_{yy} \langle w_2 \rangle + \lambda \langle w_2 \rangle = & \langle (\partial_y w_1 + V_1)(\partial_y w_1 - \partial_y s_1) \rangle \\ & + \frac{\lambda}{2} \langle (s_1 - w_1)(2u_1 - 3w_1 + 3s_1) \rangle, \end{aligned} \qquad (26)$$

with $P_2$ entering the periodic boundary conditions:

$$\langle w_2 \rangle|_{-1/2} = \langle w_2 \rangle|_{1/2} \text{ and } P_2 = -\partial_y \langle w_2 \rangle|_{\pm 1/2}. \qquad (27)$$

Again, as the right-hand side of Eq. (26) only involves first-order terms, it is sufficient to consider the first-order contributions of the controls $s_1$ and $u_1$ in order to obtain both $C_2$ and $P_2$.

Our problem is to maximize the functional $P_2$ at a given value of the functional $C_2$. To identify admissible solutions of this problem we also need to impose that

$$\langle f_1^2 + u_1^2 \rangle \le 1, \qquad (28)$$

which ensures that the spatiotemporal variations of the obtained solutions are in agreement with the asymptotic expansion (19). The upper bound in Eq. (28), which we chose to be equal to 1, controls the magnitude of the active fields $f_1$ and $u_1$, up to a renormalization of $\epsilon$.

Note that since we maximize $P_2$, we focus exclusively on positive values of the velocity $V_2$. In view of the spatial symmetry of the system of Eqs. (21) and (22), in addition to an actuation protocol ($u_1(y)$, $s_1(y)$) that gives the performance $P_2$, there is always an actuation protocol ($u_1(-y)$, $s_1(-y)$) delivering the performance $-P_2$ at the same cost. So the solution for the negative velocity corresponding to minimization of $P_2$ can be obtained from our results.

As we have already mentioned, our main goal is to investigate the role of the parameter $\lambda$, characterizing the relative importance of the chemical versus the mechanical activity, in the choice of the optimal crawling gait. More precisely, we ask how the optimal actuation strategy adjusts as $\lambda$ increases from zero to infinity.

## VII. PURELY MECHANICAL DRIVING

Our starting point is the purely mechanical (elastic) limit where the mass exchange with the chemical reservoir main-

taining the turnover is absent. In other words, we assume that $\lambda \to 0$ and study the limit when the reaction regulating material turnover is much slower than the mechanical flow of matter. In this approximation the chemical driving $u_1$ becomes irrelevant.

We will use this limiting case as a benchmark for the rest of the analysis: we will not consider any energetic cost constraint in this section and only maximize the performance regardless of the cost. The obtained maximal performance will then be associated as a benchmark to a certain cost value which will be maintained constant when we study the other cases where $\lambda > 0$. In this way, all the actuation protocols considered in this paper will be comparable as they will be characterized by the same level of injected power.

Using the convenient variables $j_1 = -\partial_y w_1$ and $f_1 = -\partial_y s_1$ we can rewrite Eqs. (21) and (22) in the form

$$\begin{aligned} \partial_t j_1 - \partial_{yy} j_1 &= \partial_t f_1, \\ V_1 &= j_1|_{\pm 1/2} \text{ and } \partial_y j_1|_{-1/2} = \partial_y j_1|_{1/2}. \end{aligned} \qquad (29)$$

The corresponding expressions for the cost [Eq. (23)] and the performance [Eq. (24)] simplify accordingly:

$$C_2^0 = \langle j_1 f_1 \rangle, \quad P_2^0 = \langle \bar{\rho}_1 j_1 \rangle, \qquad (30)$$

where $\partial_y \rho_1 = f_1 - j_1$ and $\bar{\rho}_1 = 0$. To justify the expression for $P_2^0$ in Eq. (30), we introduce the variable $j_2 = -\partial_y w_2$ and write

$$\partial_y \langle j_2 \rangle = -\langle \partial_y \rho_1 (j_1 - V_1) \rangle = -\langle \partial_y (\rho_1 (j_1 - V_1)) \rangle,$$

where we used that, in view of Eqs. (29), $\partial_t \rho_1 = -\partial_y j_1$. Hence, given that $P_2^0 = \langle j_2 \rangle|_{\pm 1/2}$, the expression in Eq. (30) follows.

To solve the linear heat equation (29) we use a standard approach and represent the mechanical driving $f_1$ in Fourier series,

$$f_1(y, t) = \sum_{l=1}^{\infty} f_1^{2l}(t) v_{2l}(y) + f_1^{2l-1}(t) v_{2l-1}(y), \qquad (31)$$

where we separated the terms containing spatially even and spatially odd modes, $v_{2l}(y) = \sqrt{2} \cos(2l\pi y)$ and $v_{2l-1}(y) = \sqrt{2} \sin(2l\pi y)$, respectively. In view of the time-periodic nature of the driving, the time-dependent coefficients $f_1^{2l, 2l-1}$ are T periodic. We can also write a similar representation for the solution of Eqs. (29),

$$j_1(y, t) = \sum_{l=1}^{\infty} j_1^{2l}(t) v_{2l}(y) + j_1^{2l-1}(t) v_{2l-1}(y),$$

whose time-dependent coefficients satisfy the equation

$$\partial_t j_1^{2l, 2l-1}(t) + \alpha_l j_1^{2l, 2l-1}(t) = \partial_t f_1^{2l, 2l-1}(t),$$

where $\alpha_l = 4\pi^2 l^2$.

At large times, the solution of the above ordinary differential equation takes the form

$$j_1^{2l, 2l-1}(t) = f_1^{2l, 2l-1}(t) - \int_0^\infty K_l(t-u) f_1^{2l, 2l-1}(u) du,$$

where the kernel is $K_l(t) = \alpha_l \exp(-\alpha_l t) H(t)$ while $H(t)$ is the standard Heaviside function. In terms of physical variables
it means that

$$ \rho_1^{2l, 2l-1}(t) = \frac{(-1)^{2l, 2l-1}}{2l\pi} \int_0^{\infty} K_l(t-u) f_1^{2l-1, 2l}(u) du. $$

We can now compute the spatial average of interest:

$$ \overline{\rho_1 j_1} = \sum_{l=1}^{\infty} \frac{1}{2l\pi} \int_{0}^{\infty} \int_{0}^{\infty} G_l(t, u, v) [f_1^{2l}(u) f_1^{2l-1}(v) - f_1^{2l-1}(u) f_1^{2l}(v)] du dv, $$

where $G_l(t, u, v) = K_l(t-u)[\delta(t-v) - K_l(t-v)]$. Finally, performing the time averaging of $G_l$, we obtain

$$ P_2^0 = \lim_{M \to \infty} \sum_{l=1}^{\infty} \int_0^M \int_0^M Q_l(u-v) \frac{f_1^{2l}(u) f_1^{2l-1}(v)}{Ml\pi} du dv, \quad (32) $$

where $Q_l(t) = -\alpha_l \text{sign}(t) \exp(-\alpha_l |t|)/2$. Using the same approach we can also compute the cost:

$$ C_2^0 = \lim_{M \to \infty} \sum_{l=1}^{\infty} \left[ \int_0^M (f_{2l}(u)^2 + f_{2l-1}(u)^2) du - \int_0^M \int_0^M L_l(u-v)(f_{2l}(u), f_{2l}(v)) + f_{2l-1}(u)f_{2l-1}(v) du dv \right], \quad (33) $$

where $L_l(t) = \alpha_l \exp(-\alpha_l |t|)/2$.

The next step is to express the time-dependent coefficients $f_1^{2l, 2l-1}$ in temporal Fourier series:

$$ f_1^{2l, 2l-1}(t) = \sum_{k=1}^{\infty} A_{2l, 2l-1}^{k} v_{2k} \left(\frac{t}{T}\right) + B_{2l, 2l-1}^{k} v_{2k-1} \left(\frac{t}{T}\right). \quad (34) $$

Substituting these expressions into Eqs. (32) and (33) and performing integration, we finally obtain

$$ P_2^0 = \sum_{l=1}^{\infty} \sum_{k=1}^{\infty} \frac{2kT \alpha_l (A_{2l}^k B_{2l-1}^k - A_{2l-1}^k B_{2l}^k)}{l(\alpha_k + \alpha_l^2 T^2)} \quad (35) $$

and

$$ C_2^0 = \sum_{l=1}^{\infty} \sum_{k=1}^{\infty} \frac{\alpha_k [(A_{2l}^k)^2 + (A_{2l-1}^k)^2 + (B_{2l}^k)^2 + (B_{2l-1}^k)^2]}{\alpha_k + \alpha_l^2 T^2}. \quad (36) $$

The constraint (28) reduces to

$$ \sum_{l=1}^{\infty} \sum_{k=1}^{\infty} (A_{2l}^k)^2 + (A_{2l-1}^k)^2 + (B_{2l}^k)^2 + (B_{2l-1}^k)^2 \le 1. \quad (37) $$

Note that, as a result of these manipulations, our original partial differential equation control problem has been reduced to an algebraic optimization problem.

Before moving to the solution of this algebraic problem, we observe that for actuation with time-reversal symmetry $f_1(t) = f_1(-t)$ all coefficients $B^k = 0$ and therefore we have from Eq. (35) that $P_2^0 = 0$. This observation can be viewed as a variant of the scallop theorem [23] in our system: the

time-reversal symmetry of the actuation must be broken for self-propulsion to become possible.

We now focus on the maximization of the performance $P_2^0$ subjected to the inequality constraint (37). To this end we introduce the matrix

$$ \mathbb{Q} = \begin{pmatrix} 0 & 0 & 0 & 1/2 \\ 0 & 0 & -1/2 & 0 \\ 0 & -1/2 & 0 & 0 \\ 1/2 & 0 & 0 & 0 \end{pmatrix} $$

and the vector

$$ U_{k,l} = \begin{pmatrix} A_{2l}^k \\ A_{2l-1}^k \\ B_{2l}^k \\ B_{2l-1}^k \end{pmatrix}. $$

In these notations, the problem is to maximize at each $k$ and $l$ the quadratic form

$$ U_{k,l}^{T} \mathbb{Q} U_{k,l} = A_{2l}^{k} B_{2l-1}^{k} - A_{2l-1}^{k} B_{2l}^{k}, $$

under the constraint that

$$ \sum_{l=1}^{\infty} \sum_{k=1}^{\infty} m_{k,l} \le 1, \quad (38) $$

where $m_{k,l} = U_{k,l}^T U_{k,l}$. The eigenspace of $\mathbb{Q}$ corresponding to its largest eigenvalue is two dimensional. With the inequality constraint (38) taken into account, this eigenspace, containing the optimal actuation modes, can be parametrized by the two sets of coefficients $a_{k,l}$ and $m_{k,l}$ and written in the form

$$ U_{k,l} = \frac{\sqrt{m_{k,l}}}{\sqrt{2}} \begin{pmatrix} \cos(a_{k,l}) \\ -\sin(a_{k,l}) \\ \sin(a_{k,l}) \\ \cos(a_{k,l}) \end{pmatrix}. $$

Substituting this expression of $U_{k,l}$ into Eq. (35), we obtain

$$ P_2^0 = \sum_{l=1}^{\infty} \sum_{k=1}^{\infty} \frac{klm_{k,l}T}{k^2 + \alpha_l(lT)^2}, $$

which does not involve the coefficients $a_{k,l}$. The remaining problem of maximizing $P_2^0$ in $k, l$ can be viewed as a problem of allocating the weights $m_{k,l}$ at integer points $(k, l)$ where the expression $kl/(k^2 + \alpha_l(lT)^2)$ reaches its largest value. There is actually only a single point $(k_0, l_0)$ where such an expression reaches its maximal value and it is then natural to set $m_{k_0, l_0} = 1$. Finding the maximum of the function $(k, l) \to kl/(k^2 + \alpha_l(lT)^2)$ when $k$ and $l$ are considered to be continuous variables and setting $\gamma = T/(2\pi)$, we find

$$ k_0 = \left\lceil \frac{1}{2} \left( \sqrt{4\gamma^2 + 1} - 1 \right) \right\rceil $$

and

$$ l_0 = \left\lceil \frac{\sqrt{3} \sqrt{\gamma (\gamma + 2\sqrt{\gamma^2 + 12}) - 3\gamma}}{6\gamma} \right\rceil, $$

where $\lceil \cdot \rceil$ is the ceiling function. This choice delivers the maximum of $P_2^0$ at a given $T$.
Finally, maximizing $P_2^0$ with respect to this remaining pa- rameter, we find that $T = k/(2\pi)$, where $k \ge 1$ is an arbitrary integer. This means that all such actuation protocols (with $k/T$ fixed) give the same optimal performance value $P_2^0 = 1/(4\pi)$ while the corresponding spatial mode is always the same one with $l = 1$. The optimal actuation protocol is then

$$f_1(y, t) = \sqrt{2} \cos(2\pi(2\pi t - y + a)), \quad (39)$$

where $a$ is an arbitrary phase. Using Eq. (36), we find that such performance level is reached at the cost $C_2^0 = 1/2$.

The resulting optimal actuation can be characterized as a traveling wave propagating from the rear to the front. This actuation strategy is in fact often observed in the motion of limbless crawlers. It has already been analytically shown to be optimal in a similar framework in Refs. [14,40]. In these previous works, which also adopted a one-dimensional setting, the definition of the cost and performance are the same as ours and the authors also consider the maximization of the performance at a fixed imposed cost. However, the friction law is more general, allowing for a strain-dependent dissipation (which would involve the dependence of our parameter $\xi$ on the dynamic variable $\rho$). An important difference with Ref. [14] is that there the active control is imposed on a discretized reference strain rather than our control of active stress. In this sense the approach of Ref. [40], also adopting that actuation is performed by internal forces, is closer to ours. The results in Ref. [40] can be also considered more general since no a priori assumption of time periodicity of the mechanical actuation is made.

Instead of focusing on these already investigated issues, in the rest of the paper we move away from the purely mechanical problem and engage chemical activation while fixing the cost at $C_2 = C_2^0 = 1/2$. Thus, we restrict the available power of actuation based on the cost which emerged from the maximal performance reachable in the purely mechanical problem. This will allow us to compare chemical and mechanical activation strategies at the same level of power delivery.

## VIII. PRIMARILY CHEMICAL DRIVING

We now turn to the opposite limit $\lambda \to \infty$ when the reaction regulating material turnover is much faster than the mechanical flow of matter. In this limit, the reaction source term in Eq. (12) dominates the mechanical drift diffusion. However, mechanics is still playing a role by affecting the chemical potential $\mu_c$ through its strain dependence which makes both $s_1$ and $u_1$ operative in the optimization problem.

In order to avoid the formation of boundary layers, we further assume that both the mechanical driving $s_1$ and the chemical driving $u_1$ satisfy periodic boundary conditions. In this case, we obtain from Eq. (21) that $w_1 - s_1 + u_1 \to 0$, physically corresponding to the equilibrium condition $\mu_c = \mu_a$ (in dimensional form).

Under these assumptions we can write

$$ \begin{aligned} C_2^\infty &= \langle \partial_y w_1 \partial_y s_1 \rangle + \lambda \langle u_1(s_1 + u_1 - w_1) \rangle \\ &= \langle \partial_y(u_1 + s_1) \partial_y s_1 \rangle + \langle u_1(\partial_t w_1 - \partial_{yy} w_1 - \partial_t s_1) \rangle \\ &= \langle \partial_y(u_1 + s_1)^2 \rangle. \end{aligned} \quad (40) $$

For the performance, Eq. (26) reduces to $\langle r_2 \rangle \to \langle u_1^2 \rangle/2$; hence,

$$P_2^\infty = -\frac{1}{2} \partial_y \langle u_1^2 \rangle |_{\pm 1/2}. \quad (41)$$

The goal now is to maximize $P_2^\infty$ with the constraint $C_2^\infty = 1/2$ while respecting the small perturbation inequality (28).

The problem is explicitly expressed in terms of the controls without the use of the auxiliary function $w_1$. Note also that time is playing a transparent role in this problem since it is involved only in the final averaging operation. It is therefore possible to first solve the optimization problem with controls that are only space dependent and then multiply the obtained solutions by any T-periodic function whose square average is equal to 1 (in order to fulfill the constraints). In other words, time and space variables can be separated in the optimal actuation protocols.

We thus expand $u_1$ and $s_1$ in Fourier series using the $v_{2l,2l-1}$ basis

$$ \begin{aligned} u_1(y) &= \sum_{l=1}^{\infty} u_1^{2l} v_{2l}(y) + u_1^{2l-1} v_{2l-1}(y), \\ s_1(y) &= \sum_{l=1}^{\infty} s_1^{2l} v_{2l}(y) + s_1^{2l-1} v_{2l-1}(y), \end{aligned} \quad (42) $$

which gives

$$ \begin{aligned} P_2^\infty &= -4\pi \sum_{l=1}^{\infty} l \langle u_1^{2l} u_1^{2l-1} \rangle, \\ C_2^\infty &= 4\pi^2 \sum_{l=1}^{\infty} l^2 \langle (u_1^{2l} + s_1^{2l})^2 + (u_1^{2l-1} + s_1^{2l-1})^2 \rangle. \end{aligned} \quad (43) $$

The inequality constraint (28) takes the form

$$ \sum_{l=1}^{\infty} 4\pi^2 l^2 \left[ \langle (s_1^{2l})^2 + (s_1^{2l-1})^2 \rangle + \langle (u_1^{2l})^2 + (u_1^{2l-1})^2 \rangle \right] \le 1. \quad (44) $$

In the above formulation, all the spatial modes can be multiplied by a time-dependent function $\phi(t)$, leaving the problem unchanged as soon as $\langle \phi^2 \rangle = 1$.

Following the same approach as in Sec. VII, we first maximize the term $u_1^{2l} u_1^{2l-1}$ under the constraints

$$4\pi^2 l^2 \left[ (u_1^{2l} + s_1^{2l})^2 + (u_1^{2l-1} + s_1^{2l-1})^2 \right] = c_l$$

and

$$4\pi^2 l^2 \left[ (s_1^{2l})^2 + (s_1^{2l-1})^2 \right] + (u_1^{2l})^2 + (u_1^{2l-1})^2 = m_l,$$

where $c_l$ is a non-negative sequence of numbers that sums to $1/2$ and $m_l$ is a similar sequence whose sum is smaller than 1.

The resulting mathematical problem can be again qualified as the maximization of a quadratic form under two quadratic constraints. Maximizing the associated Lagrangian function with two Lagrange multipliers corresponding to the two constraints, we obtain that at given $c_l$ and $m_l$, the maximal performance is

$$ P_2^\infty(l, c_l, m_l) = \frac{2\pi l}{(4\pi^2 l^2 + 1)^2} [c_l(4\pi^2 l^2 - 1) + m_l + 4\pi l(\sqrt{c_l(-c_l + 4\pi^2 l^2 m_l + m_l) + \pi l m_l})]. \quad (45) $$
As $P_2^\infty$ is decreasing as a function of $l$ (with other variables fixed) and increasing as a function of both $c_l$ and $m_l$ (also with other variables fixed), expression (45) is made maximal when the weights $c_l = 1/2$ and $m_l = 1$ are allocated to the first spatial mode at $l = 1$.

We then obtain that the optimal performance is

$$P_2^\infty = \frac{1}{1/\pi + 12\pi - 4\sqrt{1 + 8\pi^2}} \quad (46)$$

and the corresponding optimal actuation protocols are

$$\begin{align}
u_1^\infty &= \phi(t) \frac{\sqrt{1 + 4\pi(\pi + 8\pi^3 - \sqrt{1 + 8\pi^2})}}{\sqrt{2(1 + 4\pi^2)(4\pi - \sqrt{1 + 8\pi^2})}} \nonumber \\
&\quad \times \sin\left(\pi\left(\frac{1}{4} - 2y\right)\right), \\
s_1^\infty &= \phi(t) \frac{\sqrt{1 + 4\pi(\pi + 8\pi^3 - \sqrt{1 + 8\pi^2})}}{-2\sqrt{2\pi}(1 + 4\pi^2)} \nonumber \\
&\quad \times \sin\left(\pi\left(\frac{1}{4} - 2y\right)\right). \tag{47}
\end{align}$$

One can see that the control strategy maximizing the performance while maintaining the cost fixed is in this case a *standing wave* where chemical and mechanical drivings effectively conspire. Note that the maximal performance in the presence of such cooperativity is actually larger than in the case where only mechanical driving is present ($P_2^\infty \simeq 0.44 > P_2^0 \simeq 0.08$). The reason is that the possibility of fast chemically driven mass redistribution (material turnover) complements mechanical deformation in driving the internal flow, which can be then facilitated chemically even in the presence of large frictional forces.

## IX. CHEMOMECHANICAL CROSSOVER

Depending on the relative importance of chemical versus mechanical pathways there may be two very different optimal strategies of actively driving the steady crawling on a rigid substrate. If the kinetics of turnover is much slower than that of sliding friction with the substrate, crawling is optimally driven by a control in the form of a *traveling wave*. If, instead, mechanics is much slower than chemistry, the optimal control takes the form of a *standing wave*. In this section we study the crossover between these limiting control strategies by considering the case when the parameter $\lambda \ge 0$ is finite.

Our approach in this section does not address the implied infinite-dimensional optimal control problem in its full generality. Instead, we only consider a finite-dimensional version of the problem by restricting the space of admissible controls to an eight-parameter subspace generated by a particular subset of the possible controls. The motivation behind such a choice is simple: the chosen subset of controls interpolates between the optimal solutions of the general optimization problem in the two limiting cases when turnover is absent ($\lambda \to 0$; see Sec. VII) and when turnover is infinitely fast ($\lambda \to \infty$; see Sec. VIII). This choice allows one, after an appropriate finite-

<figure><img src="image_3.png" /><figcaption>FIG. 1. Optimal performance as a function of turnover (logarithmic scale). The dashed lines indicate the analytic limits $\lambda \to 0$ and $\lambda \to \infty$ computed in Secs. VII and VIII.</figcaption></figure>

dimensional optimization, to analytically find the intermediate optimal regimes.

More specifically, we now consider the following parametric form of the controls:

$$\begin{align}
s_1(y, t) &= s_{cc} \cos(4\pi^2 t) \cos(2\pi y) + s_{cs} \cos(4\pi^2 t) \sin(2\pi y) \nonumber \\
&\quad + s_{sc} \sin(4\pi^2 t) \cos(2\pi y) + s_{ss} \sin(4\pi^2 t) \sin(2\pi y), \nonumber \\
u_1(y, t) &= u_{cc} \cos(4\pi^2 t) \cos(2\pi y) + u_{cs} \cos(4\pi^2 t) \sin(2\pi y) \nonumber \\
&\quad + u_{sc} \sin(4\pi^2 t) \cos(2\pi y) + u_{ss} \sin(4\pi^2 t) \sin(2\pi y). \tag{48}
\end{align}$$

The goal now is to find the optimal set of coefficients $s_{i,j}$ and $u_{i,j}$ where the indices take the values either $c$ or $s$.

We use this specific form of the controls (48) to analytically solve the linear problem of Eqs. (21) and (22) and compute the function $w_1(y, t)$. With this information at hand we can directly express the cost (23) and also the performance (24) by explicitly solving Eqs. (26) and (27). We do not include the corresponding straightforward but cumbersome derivations here. The remaining Karush-Kuhn-Tucker problem of maximizing $P_2$ under the constraint $C_2 = 1/2$ (while respecting the inequality constraint (28)) becomes finite dimensional in terms of the height coefficients $s_{i,j}$ and $u_{i,j}$. The structure of the optimal solution can then be studied as a function of the remaining parameter $\lambda$.

In Fig. 1 we show the optimal performance $P_2(\lambda)$ computed by numerically solving the optimization problem with an interior-point method for each value of $\lambda$. We observe a gradual transition between the two limiting regimes. Thus, as $\lambda \to 0$ the function $P_2(\lambda)$ reaches the value $P_2^0 = 1/(4\pi)$ computed analytically in the $\lambda \ll 1$ limit. The corresponding optimal actuation strategy is the purely mechanical driving

$$s_1^0 = -\pi \sqrt{2} \sin(2\pi(2\pi t - y + a))/2 \quad (49)$$

without any chemical driving: $u_1^0 = 0$ found in Eq. (39). Note that the optimal performance first decays to reach a mini-
<figure><img src="image_4.png" /><figcaption>FIG. 2. Evolution of the mechanical (black) and chemical (red) costs as a function of the turnover (logarithmic scale). The sum of the two costs is maintained at the fixed value $C_2 = 1/2$.</figcaption></figure>

mum, but then finally starts to increase, reaching eventually the plateau $P_2^\infty = 1/(1/\pi + 12\pi - 4\sqrt{1+8\pi^2})$ computed analytically in the $\lambda \gg 1$ limit. The corresponding optimal actuation agrees with Eq. (47), but now, in view of the particular structure of the form (48), with an explicitly specified time-dependent multiplier $\phi(t) = \sqrt{2}\sin(4\pi^2t + a)$. In both limits $a$ is an arbitrary phase.

In Fig. 2 we show the contributions of the mechanical $\langle\Pi_2^m\rangle$ and the chemical $\langle\Pi_5^m\rangle$ activity to the total energetic cost of self-propulsion $C_2 = \langle\Pi_2^m\rangle + \langle\Pi_5^m\rangle = 1/2$. We see that in the regimes with $\lambda \ll 1$ the main energy supply comes from the mechanical source and is represented by the work of the active stresses. Instead, in the regimes with $\lambda \gg 1$ the main energy supply comes from the chemical “pumps” ensuring the appropriate target density of the turnover. Moreover, in those regimes the flux of mechanical energy even changes sign such that mechanical actuators actually work to extract energy from the system (effectively corresponding to a brake on the global motion in the absence of chemical activity). In this way more chemical energy can be injected in. This is a consequence of the constraint fixing the total chemomechanical energy input and leaving the system the freedom to self-organize to optimally use this cost.

Another observable feature of the interaction between the mechanical and the chemical actuation is the dip in the performance $P_2(\lambda)$ at small to finite values of $\lambda$ with respect to the performance $P_2^0$ achieved at $\lambda = 0$ (see Fig. 1). This means that, when first activated, the chemical machinery is detrimental as it interferes with the traveling wave mechanical activation and effectively works as a brake. As the corresponding “chemical engine” gets sufficiently strong, it starts to modify the very regime of actuation from a traveling to a standing wave type, and the performance starts to grow, reaching eventually the limit $P_2^\infty$.

To corroborate this explanation, we show in Fig. 3 the evolution of the optimal performance when we use a restricted

<figure><img src="image_5.png" /><figcaption>FIG. 3. Optimal performances as a function of the turnover with a mechanical driving only (red curve) and with a chemical driving only (blue curve). The cost is maintained at the fixed value $C_2 = 1/2$. The dashed line indicates the value of $P_2^0$.</figcaption></figure>

actuation form with either $u_1 = 0$ (purely mechanical actuation) or $s_1 = 0$ (purely chemical actuation). As expected, the purely mechanical actuation performance only deteriorates with growing $\lambda$. Indeed, the increase of material turnover still triggers the transition from a traveling to a standing wave for the optimal actuation protocol. But the latter is associated with a decrease of performance in the absence of any active chemical recycling. Instead, the purely chemical optimal actuation, which always takes the form of a standing wave, can take advantage of fast turnover while it is less effective than the purely mechanical actuation for slow turnover. The nonmonotony of the optimal performance at small to intermediate values of $\lambda$ of the purely chemical driving is due to a switch of the form (in terms of coefficients $u_{cc}$, $u_{cs}$, $u_{sc}$, $u_{ss}$) of the optimal stationary wave at a threshold value $\lambda = 8\pi^2$. It is interesting that, in the limit $\lambda \gg 1$, the purely chemical actuation achieves exactly the same level of performance $P_2^0$ as the purely mechanical actuation achieves at $\lambda \ll 1$ even with a fundamentally different spatiotemporal pattern: a traveling wave for the mechanical actuation and a standing wave for the chemical actuation.

Another illustration of the progressive transition from a traveling-wave-type to a standing-wave-type actuation can be provided if we use Eq. (48) to construct an “order parameter”-type variable $\theta$, normalized to vanish in the standing wave regime. To this end we first rewrite Eq. (48) in the form

$$s_1 = 2s_0 \cos \left( 4\pi^2 t - \frac{\psi_- + \psi_+}{2} \right) \cos \left( 2\pi y - \frac{\psi_+ - \psi_-}{2} \right) \\ + (s_+ - s_0) \cos (2\pi (2\pi t - y) - \psi_-) \\ + (s_- - s_0) \cos (2\pi (2\pi t + y) - \psi_+)$$

where $s_0 = \sqrt{s_{cc}^2 + s_{ss}^2 + s_{cs}^2 + s_{sc}^2}/2$, $\Delta s = (s_{cc}s_{ss} - s_{cs}s_{sc})/2$, $s_+ = \sqrt{s_0^2 + \Delta s}$, $s_- = \sqrt{s_0^2 - \Delta s}$, $\tan(\psi_+) = (s_{ss} + s_{cc})/(s_{sc} + s_{cs})$, and $\tan(\psi_-) = (s_{cc} - s_{ss})/(s_{sc} - s_{cs})$. This representation splits the terms representing a standing
<figure><img src="image_6.png" /><figcaption>FIG. 4. Transition from a traveling wave optimal actuation to a standing wave as $\lambda$ increases. (a) The decay of the order parameter $\theta$ to zero (logarithmic scale). (b) The spatiotemporal density variations $w_1 - s_1$ in the system at specific values of $\lambda$.</figcaption></figure>

wave contribution from the two traveling waves moving, respectively, prograde and retrograde. Note that at $\Delta s = 0$, or equivalently $s_+ = s_- = s_0$, the traveling wave contribution vanishes. Therefore, with the above definitions of $s_-$ and $s_+$, the indicator

$$ \theta = \frac{(s_- - s_0)^2 + (s_+ - s_0)^2}{4s_0^2} = 1 - \frac{s_+ + s_-}{2s_0} $$

compares the magnitude of the traveling-wave-type contributions with the standing-wave-type contribution while being normalized to vanish in the purely standing-wave-type regime. The same type of “order parameter” can be constructed using the $u_1$ part of the form (48) or indeed any linear combination between $u_1$ and $s_1$.

In Fig. 4 we show the behavior of the function $\theta(\lambda)$ in the regimes when both mechanical and chemical actuations are present. In the same figure we illustrate how the spatiotemporal density profile evolves as the parameter $\lambda$ changes from zero to infinity and the traveling wave becomes progressively arrested, opening the way to the formation of quiescent nodes separating periodically “breathing” sectors.

As we are dealing only with the simple form (48), the exact evolution of the optimal actuation regime at a finite timescale of the mechanical friction vis-à-vis the timescale of the chemical relaxation would still have to be obtained for general time-periodic controls. For instance, it cannot be excluded that the optimal transition between the two limiting regimes is abrupt rather than continuous, as it is suggested by our finite-dimensional approximation. Also, some other spatiotemporal patterns may emerge along the way as one moves from one limiting regime to another. The rigorous clarification of all these issues which require more intense numerical approaches is left for future studies.

## X. CONCLUSIONS

We have proposed a 1D prototypical model of a chemomechanically driven system which can crawl on a solid substrate. While it can be interpreted as a paradigmatic approach to the understanding of how active stress and meshwork remodeling can conspire in living cells, it can also be used as a set of guiding principles for the design of biomimetic self-propelling soft robots.

The motion is actuated by two time-periodic active controls which are fully intrinsic as they exert zero average action on the system. Those two controls are physically fundamentally different, with one being mechanical and the other chemical. The active mechanical force field is deforming the elastic scaffold while the out-of-equilibrium mass reservoir actively controls the availability of building blocks of the scaffold.

For our analysis we have chosen a simple close-to-equilibrium framework which enabled us to introduce the turnover kinetic timescale and the characteristic timescale of mechanical relaxation. When the turnover kinetics is slow, the best performance at a fixed energetic cost is achieved by a purely mechanical actuation represented by a traveling wave deformation propagating from the rear to the front of the crawler. Instead, when the turnover kinetics is fast, the best performance is reached when both mechanical and chemical drivings cooperate and form a standing wave.

Our approach also allowed us to study the continuous crossover between the limiting “mechanics-dominated” and “chemistry-dominated” regimes for a specific actuation form. In fact, one can expect that in more elaborate models various other intermediate regimes of actuation, characterizing alternative optimal crawling gaits, can become possible with continuous as well as discontinuous transitions between them. The potential complexity of this issue is already suggested by our observation that when the activity of the chemical reservoir is still weak, the associated material turnover represents itself only as a dissipative process which lowers the performance. However, when the chemical activity becomes sufficiently strong, the turnover enhances the performance by offering the possibility to recycle matter without creating a frictional counterflow. This simple example shows that the optimization of the metabolic actuation may involve a complex interplay between mechanical and chemical active agents and suggests that it is cooperativity of these two mechanisms that ultimately ensures optimality of the self-propulsion machinery.

An important remaining open question is the very possibility to separate the active controls, operating in real living systems, into a purely mechanical actuation and a purely chemical actuation, that are only connected to each other by the constraint that they should operate at a fixed total metabolic cost. In crawling cells both the mechanical activity of molecular motors exerting contractile forces on the polymer network and the chemical activity regulating the turnover of the meshwork through its polymerization and depolymerization are ultimately driven by the same chemical process: the out-of-equilibrium reaction of ATP hydrolysis. The chemical and mechanical actuations are also tightly dynamically cross-regulated through enzyme-coupled receptors [58]. In this way
the internal driving mechanisms are then coupled to the dynamics through a system of feedback loops. While all these important processes are left outside the present study, they should become a subject of future work.

ACKNOWLEDGMENTS

P.R. acknowledges the support from the French Grant No. ANR-19-CE13-0028. L.T. acknowledges the support from the French Grant No. ANR-10-IDEX-0001-02 PSL.

[1] B. Flaherty, J. McGarry, and P. McHugh, Cell Biochem. Biophys. **49**, 14 (2007).

[2] A. Mogilner, J. Math. Biol. **58**, 105 (2009).

[3] I. S. Aranson, *Physical Models of Cell Motility* (Springer, Cham, 2016).

[4] T. Svitkina, Cold Spring Harbor Perspect. Biol. **10**, a018267 (2018).

[5] C. Ahn, X. Liang, and S. Cai, Adv. Mater. Technol. **4**, 1900185 (2019).

[6] S. Chen, Y. Cao, M. Sarparast, H. Yuan, L. Dong, X. Tan, and C. Cao, Adv. Mater. Technol. **5**, 1900837(2020).

[7] Q. Ze, S. Wu, J. Nishikawa, J. Dai, Y. Sun, S. Leanza, C. Zemelka, L. S. Novelino, G. H. Paulino, and R. R. Zhao, Sci. Adv. **8**, eabm7834 (2022).

[8] F. Ahmed, M. Waqas, B. Jawed, A. M. Soomro, S. Kumar, A. Hina, U. Khan, K. H. Kim, and K. H. Choi, Smart Mater. Struct. **31**, 073002 (2022).

[9] S. Wu, Y. Hong, Y. Zhao, J. Yin, and Y. Zhu, Sci. Adv. **9**, adf8014 (2023).

[10] D. R. Yao, I. Kim, S. Yin, and W. Gao, Adv. Mater. **36**, 2308829 (2024).

[11] M. A. McEvoy and N. Correll, Science **347**, 1261689 (2015).

[12] X. Zhou, C. Majidi, and O. M. O'Reilly, Int. J. Non Linear Mech. **74**, 7 (2015).

[13] T. Shen, M. G. Font, S. Jung, M. L. Gabriel, M. P. Stoykovich, and F. J. Vernerey, Sci. Rep. **7**, 16178 (2017).

[14] D. Agostinelli, F. Alouges, and A. DeSimone, Front. Robot. AI **5**, 99 (2018).

[15] Z. Liu, Y. Wang, J. Wang, and Y. Fei, Robotica **40**, 3995 (2022).

[16] L. Giraldi and F. Jean, SIAM J. Control Optim. **58**, 1700 (2020).

[17] M. Yu, W. Yang, Y. Yu, X. Cheng, and Z. Jiao, Actuators **9**, 26 (2020).

[18] D. K. Patel, X. Huang, Y. Luo, M. Mungekar, M. K. Jawed, L. Yao, C. Majidi, D. K. Patel, L. Yao, X. Huang, Y. Luo, C. Majidi, M. Mungekar, and M. K. Jawed, Adv. Mater. Technol. **8**, 2201259 (2023).

[19] G. I. Taylor, Proc. R. Soc. London A **209**, 447 (1951).

[20] M. Lighthill, J. Fluid Mech. **9**, 305 (1960).

[21] E. Lauga and T. R. Powers, Rep. Prog. Phys. **72**, 096601 (2009).

[22] A. DeSimone and A. Tatone, Eur. Phys. J. E **35**, 85 (2012).

[23] E. M. Purcell, Am. J. Phys. **45**, 3 (1977).

[24] A. Shapere and F. Wilczek, J. Fluid Mech. **198**, 557 (1989).

[25] J. E. Avron and O. Raz, New J. Phys. **10**, 063016 (2008).

[26] R. Golestanian and A. Ajdari, Phys. Rev. E **77**, 036308 (2008).

[27] R. Golestanian and A. Ajdari, Phys. Rev. Lett. **100**, 038101 (2008).

[28] F. Alouges, A. DeSimone, and A. Lefebvre, Eur. Phys. J. E **28**, 279 (2009).

[29] F. Alouges, A. DeSimone, L. Giraldi, and M. Zoppello, IFAC-PapersOnLine **50**, 4120 (2017).

[30] F. Julicher, K. Kruse, J. Prost, and J.-F. Joanny, Phys. Rep. **449**, 3 (2007).

[31] A. C. Callan-Jones and R. Voituriez, Curr. Opin. Cell Biol. **38**, 12 (2016).

[32] P. Recho, T. Putelat, and L. Truskinovsky, Phys. Rev. Lett. **111**, 108102 (2013).

[33] P. Recho and L. Truskinovsky, Math. Mech. Solids **21**, 263 (2016).

[34] P. Recho, T. Putelat, and L. Truskinovsky, J. Mech. Phys. Solids **84**, 469 (2015).

[35] P. Recho, J.-F. Joanny, and L. Truskinovsky, Phys. Rev. Lett. **112**, 218101 (2014).

[36] A. Farutin, J. Étienne, C. Misbah, and P. Recho, Phys. Rev. Lett. **123**, 118101 (2019).

[37] Y. I. Wu, D. Frey, O. I. Lungu, A. Jaehrig, I. Schlichting, B. Kuhlman, and K. M. Hahn, Nature (London) **461**, 104 (2009).

[38] L. Valon, A. Marín-Llauradó, T. Wyatt, G. Charras, and X. Trepat, Nat. Commun. **8**, 14396 (2017).

[39] H. Fang, S. Li, K. Wang, and J. Xu, Bioinspir. Biomim. **10**, 066006 (2015).

[40] S. Santhosh and M. Serra, Phys. Rev. E **106**, 024610 (2022).

[41] E. K. Rodriguez, A. Hoger, and A. D. McCulloch, J. Biomech. **27**, 455 (1994).

[42] D. Ambrosi, G. A. Ateshian, E. M. Arruda, S. Cowin, J. Dumais, A. Goriely, G. A. Holzapfel, J. D. Humphrey, R. Kemkemer, E. Kuhl et al., J. Mech. Phys. Solids **59**, 863 (2011).

[43] A. Goriely, *The Mathematics and Mechanics of Biological Growth* (Springer, Berlin, 2017), Vol. 45.

[44] T. D. Pollard, W. C. Earnshaw, J. Lippincott-Schwartz, and G. Johnson, *Cell Biology E-book* (Elsevier Health Sciences, 2022).

[45] A. Erlich and P. Recho, J. Mech. Phys. Solids **178**, 105342 (2023).

[46] V. Deshpande, A. DeSimone, R. McMeeking, and P. Recho, J. Mech. Phys. Solids **151**, 104381 (2021).

[47] S. R. De Groot and P. Mazur, *Non-equilibrium Thermodynamics* (Courier Corporation, 2013).

[48] K. Kruse, J.-F. Joanny, F. Jülicher, J. Prost, and K. Sekimoto, Eur. Phys. J. E **16**, 5 (2005).

[49] B. Rubinstein, M. F. Fournier, K. Jacobson, A. B. Verkhovsky, and A. Mogilner, Biophys. J. **97**, 1853 (2009).

[50] R. J. Hawkins, M. Piel, A.-M. Lennon-Dumenil, J.-F. Joanny, J. Prost, and R. Voituriez, Biophys. J. **96**, 197a (2009).

[51] A. Carlsson, New J. Phys. **13**, 073009 (2011).

[52] K. Doubrovinski and K. Kruse, Phys. Rev. Lett. **107**, 258103 (2011).

[53] D. Shao, H. Levine, and W.-J. Rappel, Proc. Natl. Acad. Sci. USA **109**, 6851 (2012).

[54] A. Callan-Jones and R. Voituriez, New J. Phys. **15**, 025022 (2013).

[55] E. Tjhung, A. Tiribocchi, D. Marenduzzo, and M. E. Cates, Nat. Commun. **6**, 5420 (2015).

[56] C. Giverso and L. Preziosi, Phys. Rev. E **98**, 062402 (2018).

[57] B. Amiri, J. C. Heyn, C. Schreiber, J. O. Rädler, and M. Falcke, Biophys. J. **122**, 753 (2023).
[58] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter, *Molecular Biology of the Cell*, 4th ed. (Garland Science, New York, 2002).

[59] J. Prost, F. Jülicher, and J.-F. Joanny, *Nat. Phys.* **11**, 111 (2015).

[60] F. Jülicher, A. Ajdari, and J. Prost, *Rev. Mod. Phys.* **69**, 1269 (1997).

[61] M. Nishikawa, S. R. Naganathan, F. Jülicher, and S. W. Grill, *eLife* **6**, e19595 (2017).

[62] X. Qin, E. Hannezo, T. Mangeat, C. Liu, P. Majumder, J. Liu, V. Choesmel-Cadamuro, J. A. McDonald, Y. Liu, B. Yi *et al.*, *Nat. Commun.* **9**, 1210 (2018).

[63] G. B. Blanchard, J. Étienne, and N. Gorfinkiel, *Curr. Opin. Genet. Dev.* **51**, 78 (2018).

[64] M. J. Lighthill, *Commun. Pure Appl. Math.* **5**, 109 (1952).

[65] R. García-García, P. Collet, and L. Truskinovsky, *Phys. Rev. E* **100**, 042608 (2019).