Motility and morphodynamics of conﬁned cells

We introduce a minimal hydrodynamic model of polarization, migration, and deformation of a biological cell conﬁned between two parallel surfaces. In our model, the cell is driven out of equilibrium by an active cytsokeleton force that acts on the membrane. The cell cytoplasm, described as a viscous droplet in the Darcy ﬂow regime, contains a diffusive solute that actively transduces the applied cytoskeleton force. While fairly simple and analytically tractable, this quasi-two-dimensional model predicts a range of compelling dynamic behaviours. A linear stability analysis of the system reveals that solute activity ﬁrst destabilizes a global polarization-translation mode, prompting cell motility through spontaneous symmetry breaking. At higher activity, the system crosses a series of Hopf bifurcations leading to coupled oscillations of droplet shape and solute concentration proﬁles. At the nonlinear level, we ﬁnd traveling-wave solutions associated with unique polarized shapes that resemble experimental observations. Altogether, this model offers an analytical paradigm of active deformable systems in which viscous hydrodynamics are coupled to diffusive force transducers. DOI:


I. INTRODUCTION
Biological cells provide a striking example of far-fromequilibrium, highly deformable microscopic systems. During embryonic development, immune responses, or diseases such as cancer, cells undergo vital functions involving morphological changes. These include cell division and cell migration, two hallmarks of life that are driven by a perplexing array of coordinated mechano-chemical processes.
Individual cell motility takes various forms depending on the cell type and the environment in which the cell is migrating. The specifics of crawling motility are influenced by factors such as adhesion strength, the type of substratum (tissues and gels), external migratory cues (both chemical and mechanical), and the organization of the cellular cytoskeleton [1]. The cytoskeleton, a dynamic network of interlinking protein filaments, is responsible for generating the molecular forces necessary for active (ATP-fueled) translocation. By interacting with the external environment, the cytoskeleton produces locomotion and cell shape changes that correlate with internal flows and shifting distributions of intracellular biomolecules [2]. While the key force-generating components in the cell have long been identified, it remains unclear how their local activity and transport dynamics determine the macroscopic cell polarity and morphology.
In light of the many cell motility phenotypes that have been observed, two main migration modes are put forward in the literature: (1) the mesenchymal mode, characterized by strong specific adhesion to the environment and critically powered by actin polymerisation, and (2) the amoeboid mode, which in extreme cases involves only weak nonspecific interactions * Corresponding author: idolavi1@gmail.com with the environment and is characterised by high contractility levels of the cytoskeleton [1,3]. The former is typically associated with actin-rich structures that protrude form the leading edge of the cell, while the latter is associated with smoother, more rounded morphologies. Historically, most in vitro experimental studies have focused on mesenchymal cell migration over flat two-dimensional (2D) surfaces. More recently, it has been recognised that in vivo cell migration mainly occurs in three-dimensional (3D) micro-environments that impose stringent spatial constraints [4]. In fact, it was shown over the last decade that environmental cues, and in particular geometric confinement, could trigger dramatic phenotypic changes and induce prototypical amoeboid motility in several cell types [4][5][6][7][8].
For the most part, previous droplet-based models of cells have described aqueous suspensions of cytoskeleton filaments with (or without) their associated motors [16,17,21,24,37,41,[51][52][53]. To the best of our knowledge and understanding, the existing 2D models have lacked an explicit hydrodynamic account of the internal cytoplasmic fluid, meaning the cytosol. Yet this viscous fluid, which typically constitutes the bulk of the cell mass, imposes important mechanical constraints on the overall dynamics, namely, (1) in certain situations, the fluid viscosity could be a considerable source of kinetic energy dissipation, (2) the interaction of the fluid with the external substratum may lead to a marked transfer of momentum, and (3) the incompressibility of the fluid means that cell volume is controlled by water fluxes through the semipermeable cell membrane. Although a strict conservation of volume may not apply to all migrating cells, there exists experimental evidence of area conservation during keratocyte migration on substrates [46] and HL-60 (leukemia cells) migration in confinement [54].
In addition, one should be mindful of the fact that the cytoplasm is mechanically linked to the active cell cortex. Indeed, cell deformations induced by cortex activity naturally produce internal cytoplasmic flows [32,[55][56][57]. Moreover, it has been shown that cytoplasmic flows can shift the distribution of proteins that regulate cell cortex activity [58]. With regard to cell motility, such flows are thought to play an important role due to the hydrodynamic forces they apply, as well as to their influence on the transport of cytoskeleton components [55].
Here we introduce a minimal active droplet model of motility under quasi-2D confinement. Unlike previous descriptions, we tackle the problem from the perspective of the passive cytoplasmic fluid. More specifically, we describe the cell cytoplasm as a mass-conserving viscous fluid that is driven on its boundary (the cell membrane) by an active traction force. This force, which could be powered either by actin polymerisation or by actomyosin contraction, is controlled by a diffusive cytoplasmic component (a force-transducing solute). Closing the coupling, we account for the advection of the solute by the internal cytoplasmic flow.
This paper is structured as follows. Section II outlines the formulation of the model, its general properties, and a reduced dimensionless representation. Section III provides a rigorous linear stability analysis characterizing the reststate solution. This analysis shows that a global polarizationtranslation (motility) mode becomes unstable beyond a critical threshold of solute activity. Upon increasing activity further, the system crosses a series of Hopf bifurcations of multipolar modes, destabilizing coupled shape-concentration waves. In Sec. IV we introduce a nonlinear semi-analytical method for computing self-consistent traveling solutions. The results obtained in this section reveal how the steady-state speed and shape of traveling cells vary as functions of model parameters. In Sec. V we conclude the analysis and discuss the significance of our theoretical findings in light of experimental observations in vitro.

II. MODEL
We consider a fluid droplet of viscosity μ confined between two parallel plates separated by a gap h (a Hele-Shaw cell), as illustrated in Fig. 1. Let u = u(t, x, y) be the gap-averaged planar flow and p = p(t, x, y) be the internal fluid pressure FIG. 1. Model illustration. An active actomyosin force (dark arrows) drives the motility and morphodynamics of the confined cytoplasmic droplet (blue). This force acts on the droplet freeboundary in the normal direction and is modulated locally by a diffusive solute (green circles). Inset: Cross section highlighting the fluid flow and solute transport within the droplet. Thin arrows mark the parabolic flow profile, induced by the pressure gradient (blue-level background) and averaged by u (thick blue arrow). The solute binds on and off the plates with rates k on , k off and is advected by the fluid flow only in the unbound state.
(up to a constant). The 2D moving-boundary-value problem for the droplet is given by In Eq. (1), M = h 2 /12μ represents the mobility in the effective Darcy's law, which averages the Stokes momentum balance in thin films. In Eq. (2) we impose the fluid incompressibility, which [along with Eq. (1)] mandates that the pressure is Laplacian in the droplet domain. Equation (3) defines the normal force balance on the droplet free-boundary ∂ (t ) (discarding exterior pressure gradients by assuming negligible viscosity of the surrounding fluid). Here σ denotes the surface tension and κ denotes the curvature. The Young-Laplace condition is perturbed in our model by an active traction force, F act (c)n, where n is the unit normal pointing outward (see Fig. 1). This force is defined per unit length and is controlled locally by the gap-integrated concentration of an internal solute, c = c(t, x, y). We stress that F act (c) can be either positive (pushing outwards) or negative (pulling inwards). Due to incompressibility, any uniform term F 0 ∈ R added to this force would merely offset the pressure p by a constant (−F 0 /h) and thus be irrelevant to the dynamics. For the sake of generality, we do not specify an explicit form of F act (c) at this stage. In Eq. (4) we present the kinematic condition, stating that the normal velocity of the sharp interface, V n , is given by the normal velocity of the fluid on ∂ (t ).
The active force in our model represents a substrate reaction to cytoskeleton activity. By taking only a boundary force, we essentially neglect friction between the cytoskeleton and the cytosol in the bulk. This is justified if, for instance, (1) the actin filaments are relatively stationary in the laboratory reference frame or (2) the actin flow is rapid but its density is sufficiently low as to not modify Eq. (1). Importantly, our formulation does not imply that active cytoskeleton stresses and traction forces are not at play in (t ). Rather, we suggest that the net effect of bulk cytoskeleton activity on the cytosol is transmitted through the membrane. In this regard, we note that Eq. (2) constrains the problem in such a way that if one introduces an inhomogeneous active stress in (t ), which is isotropic, our model would remain equivalent [by simply redefining the pressure in Eqs. (1) and (3)].
In order to produce motion and/or deformation of the droplet, the active force must be graded along the boundary. The simplest choice is then F act (c)n, meaning a normal force that is controlled locally by one chemical species. As we focus on the dynamics of the cytosol, we do not describe the force-generation and force-transmission mechanisms at the molecular level. We postulate that our coarse-grained F act (c)n is driven either by actin polymerization against the membrane or by the contraction of actomyosin filaments. In order to transmit the external force on the cell membrane, the actin filaments must also interact mechanically with the substrate, e.g., via cortex-substrate adhesion.
To this end, the solute can be any cytoplasmic protein controlling the force-generation or adhesion machinery. In our model, we assume that c either induces an inwards pulling force or inhibits an outwards pushing force (see Fig. 1). In line with this principal assumption, we define the negative linear response of F act (c) to solute deviations about the mean planar concentration c 0 , namely, F act (c 0 ) ≡ −χ , where χ > 0.
In the bulk, we assume fast solute adsorption on the top and bottom plates (or onto an adhered cortex), as shown in Fig. 1. With rapid on and off rates, k on and k off , the quasi-2D transport of the gap-integrated solute is given by D∇c · n + aV n c = 0 on∂ (t ), where a = k on /(k on + k off ) is the steady fraction of adsorbed molecules not advected by the average flow and D is an effective diffusion coefficient. In Eq. (5) we present the effective advection-diffusion dynamics in (t ). Note here that the total (advective + diffusive) solute flux is j = (1 − a)uc − D∇c. In Eq. (6) we impose zero solute flux on the moving boundary, i.e., j · n − V n c = 0, where we inserted the kinematic condition, Eq. (4). Simply put, the solute is effectively advected at a slower velocity than that of the fluid. Hence, its concentration decreases (increases) towards an advancing (retracting) front.
Mass conservation. Our equations of motion explicitly conserve both the fluid volume (and hence the projected droplet area A ≡ da) and the total solute (C tot ≡ c da). Specifically,Ȧ = 0 results from the fluid incompressibility and the kinematic condition, Eqs. (2) and (4), whileĊ tot = 0 is due to the lack of unbalanced reactions (sources and sinks) in (t ) along with the no-flux condition on ∂ (t ), Eqs. (5) and (6).
External force balance. By calculating the time-derivatives of the droplet geometric moments (through a complex-number mapping, detailed in Appendix A), we obtain via the first moment the fluid center-of-mass velocity, u cm = A −1 ∂ xV n dl, in terms of the solute distribution on the boundary It is easy to recognize that the line integral on the RHS is precisely the net active traction force that acts on ∂ . This driving force is balanced by the net effective friction that acts over . It can be shown that the term −M −1 hAu cm equals the total viscous shear force applied externally on the two fluid layers in contact with the no-slip plates (see Appendix A). Nondimensionalization. We reduce the number of independent model parameters by formulating our physical equations in terms of rescaled dimensionless variables. In this representation, length will be given in units of R 0 = √ A/π (the droplet radius), time in units of R 2 0 /D (solute diffusion time over the droplet), solute concentration in units of c 0 = C tot /A, and pressure in units of D/M. Accordingly, the curvature (κ) will be given in units of 1/R 0 and velocities (u, V n ) in units of D/R 0 . In addition, F act will be given in units of Dh/M. We then make the transformation F act → F act /h to simplify the notation in Eq. (3).
Our reduced dimensionless equations are given explicitly in Appendix B. They are defined with three dimensionless parameters: For brevity, we shall omit the primes hereinafter, such that σ, χ will denote the dimensionless parameters in Eq. (8), unless stated otherwise.

III. LINEAR STABILITY ANALYSIS
In this section, we examine the shape-concentration dynamics close to the circular homogeneous rest state, which is a straightforward solution to our problem for all parameter values. In this state, the shape is defined by R = 1 and the internal solute concentration is c = 1. It follows that κ = 1 on ∂ and the resulting pressure is constant: p = σ − F act (1) in , meaning that u = −∇ p = 0. We now perturb both the shape of the droplet and the solute concentration such that R = 1 + δR(t, θ ) and c = 1 + δc(t, r, θ ), where δR 1 and δc 1. Let us expand these perturbations in normal cosine modes, where m = 0, 1, 2, . . . , and δR m 1, δc m 1. Note that the sine modes are omitted here since orthogonal perturbations are uncoupled at the linear level.
From the kinematic condition, Eq. (4), it follows that V r −∂ r δ p| r=1 , and thus where we used Eqs. (9) and (10) and the pressure expansion. Note that for χ = 0 we recover the cubic dispersion relation that characterizes the morphological stability of the passive droplet.
We proceed by expanding the solute transport problem, Eqs. (5) and (6), to first order in {δR m , δc m }: Note that, while the bulk equation [Eq. (12)] neglects the quadratic advection term, the linearized no-flux condition [Eq. (13)] still couples the solute transport problem with the droplet shape evolution [Eq. (11)]. Taken together, Eqs. (11)-(13) describe a closed dynamical system for the pair of cos(mθ ) perturbations, δR m (t ) and δc m (t, r). Discarding solutions that diverge at r = 0, the kernels of Eq. (12) may be written as J m (−i √ sr)e st , where J m is the Bessel function of the first kind of order m and s is the eigenvalue (constrained by the boundary condition and normally real-negative for pure diffusion problems). We search for the coupled shape-concentration eigenmodes by imposing a shared growth rate s for both degrees of freedom, where s, α ms , β ms ∈ C and we consider the real part of the RHS in both equations. As a final reduction of the linearized system, we substitute the ansatz, Eq. (14), back in Eqs. (11) and (13): The eigenvalues of Eqs. (15) and (16) are then computed as the roots of the complex characteristic function, The implicit equation G m (s) = 0, which determines s and thus governs stability, depends only on two control parameters: σ , the stabilizing surface tension, and aχ , the destabilizing droplet-solute coupling. Note that aχ also has the meaning of a Pe number, as it represents the ratio between the driven advection rate of the solute (relative to the fluid) over its diffusive transport rate. For each eigenvalue s, the associated shape-concentration eigenmode, v ms (r, θ ) ≡ {α ms , β ms J m (−i √ sr)} cos(mθ ), depends on all three independent parameters.
For aχ = 0, we recover both s = −σ m(m 2 − 1) and the infinitely many real-negative zeros of The former characterizes the morphological stability of the viscous droplet and the latter characterize the decaying diffusion modes on a rigid noflux disk. As a function of aχ , the solute diffusion modes become increasingly entwined with the droplet deformation mode (for more details, see Appendix C). In our analysis, we will focus on those coupled eigenmodes that can be destabilized at some critical aχ . These are summarized in the stability phase diagram shown in Fig. 2. We will also take note of the marginally stable eigenmodes (associated with s = 0) as these represent symmetries in our system. Mass modes. Substituting m = 0 in Eq. (17), we find that the roots of G 0 (s) are all real and nonpositive, meaning that there are no azimuthally symmetric instabilities. The eigenvalue s = 0 has two nontrivial eigenmodes: v R 00 (r, θ ) = {1, 0} and v C 00 (r, θ ) = {0, 1}. These correspond to perturbations of the fluid mass (expanding the droplet radius via v R 00 ) and the solute mass (elevating the uniform concentration via v C 00 ). Both are marginally stable due to the conservation of fluid volume and total solute.
Polarization-translation modes. Substituting m = 1 in Eq. (17) and expanding G 1 (s) for small s, we obtain the roots s = 0 and s 1 8(aχ−1) 3−aχ . Here s = 0 is associated with a single nontrivial eigenmode, v 10 (r, θ ) = {1, 0} cos(θ ), describing the infinitesimal translation of the circular homogeneous rest state (a signature of translational invariance). The more interesting eigenvalue, s 1 , changes sign from negative to positive as the control parameter aχ exceeds the value of 1 (s 1 4(aχ − 1) when aχ ≈ 1, note the black line in Fig. 2). Assuming small s 1 , the associated eigenmode is approximated by v 1s 1 (r, θ ) ≈ {χ, −s 1 r} cos(θ ). That is, the global translation of the droplet is coupled to an internal solute polarity, as suggested already by Eq. (7).
Multipolar shape-concentration modes. For m 2, we find marginally stable modes (with s = 0 for any aχ ) only if σ = 0. In fact, in this zero surface tension limit (treated in Appendix C), all normal modes are closely analogous to mode m = 1, with the additional real eigenvalue approximated by s m 4m(m+1)(aχ−1) 2+(1−aχ )m 2m(m + 1)(aχ − 1) for aχ ≈ 1. However, for σ > 0, we find a pair of complex-conjugate eigenvalues whose real part changes sign from negative to positive as aχ increases or as σ decreases (when aχ > 1); see destabilize each m 2 mode through a Hopf bifurcation. To do so, it must work not only against the solute diffusion (as in m = 1), but also against the surface tension. To first order in σ , the mth Hopf bifurcation occurs at aχ c 1 + 3m 2 +m−4 4(m+2) σ with s m ±im(m + 1) √ 2(m − 1)σ (see Appendix C for derivation and the thin translucent lines in Fig. 2). Then, the associated eigenmodes, v ms m (r, θ ), represent shape-concentration standing waves oscillating with frequency ω m = |Im(s m )|. Note that, in terms of the dimensional parameters, the oscillation period scales like T ∼ √ τ solute τ droplet , where τ solute = R 2 0 /D is the solute diffusion time over the droplet and τ droplet = R 3 0 M −1 σ −1 is the characteristic shape relaxation time of the viscous droplet.
From a physical standpoint, the standing waves result from a dynamic interplay between the active driving force and the restoring surface tension. Consider an initial state in which finite solute and curvature gradients produce force variations on the boundary that directly balance each other [ Fig. 2(c), t = 0]. At this instant, the resultant pressure is uniform and thus u = 0. The solute then spreads out via diffusion while the tension persists in countering the curvature gradients (t = T /8). As the boundary is driven by tension, retracting edges amass solute while advancing edges disperse with solute [see Eq. (6)]. By the time the circular shape is recovered (t = T /4), a residual solute gradient on the boundary induces a nonuniform active force that pulls in the retracting edges further. In essence, this overshoot (or memory) comes from the diffusive transport of the solute. As the droplet is deformed in the transverse direction by F act (c), the tension counters the new deformation until the fluid stalls again (t = T /2). The same mechanism then drives the droplet back to the initial state (as seen in Supplemental Movie 1.a [59]). We stress that such oscillations can be stable (damped) or unstable (amplified), depending on the strength of the shape-concentration coupling (see diagram in Fig. 2).
Traveling morphological waves are also supported by the model in the linear regime. Such waves can be constructed in a straightforward manner by superimposing two orthogonal standing waves [cos(mθ ) and sin(mθ )], evolving at the same amplitude and frequency with a temporal phase shift of a quarter period (see Supplemental Movies 1 and 2 [59]). Interestingly, in these traveling waves, the fluid pathlines circulate locally over time while the instantaneous streamlines remain irrotational by definition (being that u = −∇ p). Our model sustains such flow patterns because the inhomogeneous active force [F act (c)n] continuously drives the droplet in a timedependent manner.

IV. STEADILY MOVING STATES
Reverse-engineering problem. In this section we search for self-consistent solutions to our nonlinear model that are characterized by a fixed stationary shape, where we relied on the kinematic condition, Eq. (4). We prove in Appendix D that Eq. (18) mandates that the internal Darcy flow is uniform and given by u = u cm in (t ). Hence, at steady state, we take u = u 1x without loss of generality. In this state, both the pressure p and the solute concentration c are stationary in the moving frame of reference. It follows that Eq. (1) and Eqs. (5) and (6) are solved by the unique (y-independent) forms: wherex = x − u 1 t and p 1 , c 1 are normalization constants.
Hereinafter, we shall work in the moving frame of reference and omit the tilde inx for brevity. The generic steady-state forms in Eq. (19), which follow from the fixed-shape condition, Eq. (18), significantly reduce the degrees of freedom in our problem. Indeed, if a steadily moving solution exists for given set of parameters, it is fully determined by the three constant numbers (u 1 , p 1 , c 1 ). The challenge is to find those combinations of numbers that produce a self-consistent (physically viable and dimensionless) shape that is sustained by Eq. (3).
Active force saturation. Working in the nonlinear regime, one must specify a formula for F act (c). To account for the plausible saturation of this active force, we choose here the simplest Hill function (with Hill coefficient n = 1), as adopted in previous one-dimensional (1D) models [44,45], where β denotes the maximal pulling force induced by the solute and c s is the saturation parameter. In the dimensionless model, the phenomenological parameters β and c s are defined as Mβ Dh and Ac s C tot , respectively (in terms of the dimensional parameters). The linear response parameter χ , which determines the stability of the rest state, is then given by Optimization procedure. Substituting Eqs. (19) and (20) back in the normal force balance, Eq. (3), gives the curvature on ∂ as a single-valued function of x: This equation translates to an ODE problem for the boundary line that we solve numerically. To be physically viable, the resultant contour ∂ must be a simple (nonintersecting) closed curve without cusps. In addition, the enclosed domain must satisfy the dimensionless normalization conditions, namely, A = π and C tot = π . Together, these constraints translate to three nonlinear equations in the three numbers that make up our reverse-engineered solution (u 1 , p 1 , c 1 ). To solve this problem efficiently, we built a computational optimization procedure (for details, see Appendix D).
Computational solutions. In Fig. 3 we fix a moderate finite tension (σ = 1) and trace the dimensionless steadily moving states as they branch out from the rest state (u 1 = 0) through a pitchfork bifurcation that occurs at aχ = 1 (the threshold of the motility-mode instability). As shown in Fig. 3, and also analytically in Appendix D, this bifurcation is supercritical if c s 0.5 and subcritical if c s < 0.5. For c s < 0.5, there is also a saddle-node bifurcation of finite velocity occurring at some aχ * < 1. This bifurcation structure admits three motility phases: (1) rest phase, where only the stable rest state exists (aχ < 1 and c s 0.5 or aχ < aχ * and c s < 0.5), (2) traveling phase, where the rest state is unstable and there exists a stable traveling state (aχ > 1), and (3) bistable phase, where both the rest state and a high-speed traveling state are stable (aχ * < aχ < 1 and c s < 0.5). In this phase, we  Fig. 3. Note here the additional saddle-node bifurcation occurring at high χ (thin gray line). This bifurcation annihilates the stable branch and gives rise to a new unstable branch beneath it. The latter is itself cut off at a pinched state, traced on the diagram by the thin orange curve. As the two saddle-node bifurcations disappear simultaneously for low c s , the subcritical branch merges with the unstable pinching branch. Large symbols mark the representative solutions visualized in panels (a)-(c). In these panels, graphic objects and colors have the same meaning as in Figs. 3(a)-3(c). In panel (b), the second (saddle-node) and third (pinched) states correspond to χ 1.063 and χ 1.058, respectively. All other parameters can be inferred from the diagram.
force β = χ (c s + 1) 2 /c s . This picture is robust for moderate to high tension, extending to the limit σ → ∞ (the rigid circular case, also treated in Appendix D). Note that the same bifurcations and consequent motility phases appear similarly in the analogous 1D model [44,45].
We stress that our results for the deformable droplet do not merely reconstitute a familiar force-speed relationship. Using our computational optimization method for resolving the self-consistent traveling states, we also capture their shapes [Figs. 3(a)-3(c)]. This acquisition allows us to explore the force-shape (or speed-shape) relationship efficiently at the nonlinear level, without relying on heavy numerical simulations of the coupled moving-boundary dynamics.
In Fig. 4 we fix a low tension (σ = 0.1) and compute the steadily moving solutions in the same manner. Here, the bifurcation diagram exhibits a distinct qualitative change in the force-speed dependence. Strikingly, we find a parametric regime (about c s = 0.5), for which the steady-state speed u 1 (on the stable branch in blue) becomes a nonmonotonic function of the force amplitude β. The puzzling decrease in speed ends at an inverted saddle-node bifurcation, occurring at some aχ * * > min[aχ * , 1] (depending on c s ). This critical point essentially annihilates the stable branch for higher force amplitudes. The newly formed unstable branch is terminated as well, at some χ < χ * * , where a pinching point occurs, implying a topological singularity [see Fig. 4(b)]. In addition, there exists a critical c s (smaller than 0.5), below which the two saddle nodes are annihilated and the stable branch completely disappears. At this stage, the unstable pinching branch merges with the subcritical branch that originates at aχ = 1.
The recovered shapes in Fig. 4 reveal that this intricate bifurcation structure is caused by a particular deformation tendency that takes hold in the low c s (and low σ ) regime. This provides a unique opportunity to gain insights into the nonlinear physics underlying the force-shape-speed dependence.
Interpretation of the steady state morphology. We may infer the qualitative makeup of ∂ from a brief overview of the curvature extrema. Due to reflection symmetry about the axis motion (the polarity axis x), the front (rightmost) and rear (leftmost) limits of the boundary are themselves local curvature extrema. As detailed in Appendix D, both coordinates, denoted x R and x L respectively, are determined implicitly by the constraints which grant a self-consistent solution. Deriving Eq. (21), we also obtain a maximum (x + ) and minimum (x − ) of κ (x), We stress here that u 1 , c 1 (themselves determined by the constraints) also depend nonlinearly on all dimensionless parameters, including σ . We also note that x ± are necessarily real numbers since the minimal aβ for which steady traveling states exist is 4.
The extrema x ± have the meaning of crossover points in the normal force balance, Eq. (21). For x > x + , the curvature decreases with x because the pressure, given in Eq. (19), dominates the equation [F act (c) is negligible at high x due to the exponential drop of c(x), Eq. (19)]. For x ∈ (x − , x + ), the curvature increases with x due to the nonlinear rearwards amplification of the active pulling force [in this regime, For high c s , we generally find that x − < x L < x + < x R . This means that x − is irrelevant (out of bounds) and both x L , x R are local curvature minima. The resulting shape is thus elliptical-shortened along the polarity axis-as seen in Figs. 3(a) and 4(a). A decrease in c s prompts the saturation of F act as a function of x, giving x L < x − < x + < x R . In this case, the front is still a curvature minimum but the rear is now a curvature maximum. The shape is thus triangulated with a bulged rear, as seen in Fig. 3(b). Notice that, about x − , the curvature can also be negative, as seen in Fig. 3(c) (high-speed stable state on the right) and Fig. 4(b). For slow traveling states at low c s , we typically find that x L < x − < x R < x + . In this arrangement, x + is irrelevant and both x L , x R are curvature maxima. The corresponding shape is then necessarily elongated in the direction of motion, as seen in Fig. 3(c) (leftmost steady state) and Fig. 4(c). At low tension and low c s , the force amplitude β amplifies the negative curvature deformation along the stable branch, leading to a peanut-like (cat-tongue) shape, as evident in the saddle-node state [second image in Fig. 4(b)]. This introduces a boundary section, just prior to the rear, on which F act (c) is saturated and yet n x > 0. This adds a meaningful negative contribution, F act (c)n x < 0, to the external force balance, Eq. (7). Hence, the speed u 1 ultimately decreases as a consequence of this (deformation-induced) active drag effect. The same effect also leads to the morphological pinching point that terminates the newly formed unstable branch [bottom state in Fig. 4(b)].

V. CONCLUSION AND DISCUSSION
Our coupled moving-boundary model offers a concise physical description of symmetry breaking, active-capillary waves, and steady-state motility in confined cells. While such phenomena usually entail significantly more complex and analytically intractable modeling, our simple equations of motion capture two types of linear instabilities as well as a range of stationary patterns that connect motion with shape. Despite their intricacy, these results are accordant with stringent constraints imposed by (1) the incompressibility of the internal fluid and (2) the external geometry (or Darcy's law more generally).
The analysis provides both quantitative predictions (i.e., particular speeds and shapes for given parameters) and physical insights, meaning an intuitive understanding of the mechano-chemical self-organization of our system. Indeed, given the simplicity of our formulation, and being that F act (c) is defined phenomenologically, the main value of our work is in these physical insights.
Quantitative predictions for actual cell behaviors remain beyond the scope of the present study. Notwithstanding, the theoretical patterns obtained here do bear some strong qualitative similarities to a number of experimental observations. Specifically, the elongated traveling shape solutions, found in the low tension regime, are reminiscent of in vitro ameoboid motility in quasi-2D confinement [6][7][8]54]. Note the distinct similarity between the polarized progenitor cells in Ref. [7] and the predicted pearlike shape in Fig. 4(b) (first steady state). Our model also captures the formation of a uropodlike structure at the cell rear, which is typical to various motile cell types, particularly under confinement [6,8,54]. The high aspect-ratio shapes, represented in Figs. 3(a), 3(b), and 3(c) (specifically, the high-speed steady state), bear more resemblance to confined dendritic cell migration [44] or "mesenchymal" keratocyte migration on substrates [46,60].
In addition, standing and traveling normal waves have been observed in various systems, including Dictyostelium cells [61], suspended fibroblasts [62] (previously modeled in Ref. [63]), developing embryonic cells [64], migrating micogrlia cells [65], and synthetic membrane vesicles containing the Min protein system [66]. Although these experiments were not set in 2D confinement, the similarity to our coupled multipolar oscillations suggest the involvement of a similar mechanism, that is, a fluid-mediated interplay between the restoring surface tension force and an active pushing (pulling) force controlled by a diffusive inhibitor (activator).
To this end, there are still open questions pertaining to the nonlinear dynamics and limit behaviors in our system. Notably, the global stability of the resolved steadily moving states is not guaranteed. Such states may have only a limited basin of attraction within the vast shape-concentration phase space. To contemplate this point, imagine a randomly perturbed rest state that is linearly unstable with respect to numerous normal modes. In this scenario, the growing multipolar waves would ultimately couple to each other and to the motility mode at the nonlinear level. Plausible additional attractors such as shape-concentration limit cycles could conceivably compete with the steadily moving state for selection. Furthermore, the low tension regime presents a particularly puzzling conundrum: when the rest state is unstable and our branch of stable traveling solutions no longer exists, the dynamic attractor of the motility mode is completely unknown. The possibility of incurring asymmetric fragmentation via a finite-time topological singularity (pinch-off point) is in itself an exciting avenue for further investigation. Such unsettled questions may be answered with an appropriate numerical simulation of the problem.

APPENDIX A: MOMENTS
Here we derive exact relations characterizing the evolution of complex geometric moments associated with (t ). Similarly to Ref. [21], the (x, y) coordinate is mapped to the complex number z = x + iy and the kth moment is given by It is easy to recognize that M 0 = A and M 1 = AR A , where A is the droplet area and R A is the center-of-mass coordinate over the complex plane. We compute the time derivative of M k (t ),

External force balance (k = 1)
Substituting k = 1 in Eq. (A2) giveṡ Over the complex plane,Ṙ A represents the center-of-mass velocity and ∇z · n = n x + in y is the unit normal pointing outwards on ∂ (t ). Multiplying Eq. (A3) by h/M, and going back to R 2 , we may restate Eq. (A3) as where u cm denotes the droplet center-of-mass velocity over R 2 . Here we used the fact that, over a simple closed curve, the line integral of the curvature vector vanishes: κn dl = 0. We recognize that Eq. (A4) represents the external force balance on the droplet. The second term on the LHS equals the net active traction force acting on the droplet boundary in the normal direction. We will now show that the first term equals the net viscous shearing force applied externally on the two fluid layers in contact with the no-slip plates. This dissipative force is given by where τ = μ[∇u + (∇u) T ] is the 3D viscous shear stress tensor and u is the 3D flow field. Given the parabolic profile of the flow, u(x, y, z) = − z(h−z) 2μ ∇ p(x, y), we obtain τẑ = 2z−h 2 ∇ p(x, y), and thus (A6)

APPENDIX B: DIMENSIONLESS MODEL
Here we summarize the equations of motion in our reduced dimensionless model. Upon rescaling variables and parameters as discussed in Sec. II (and taking F act → F act /h for brevity), the dimensional Eqs. (1)-(6) are now given by ∇c · n + aV n c = 0 on∂ (t ). (B6) The linear response of F act (c) to solute deviations about c = 1 (the dimensionless mean concentration) is The system of Eqs. (B1)-(B7) is defined by the three dimensionless parameters given in Eq. (8).
The dimensionless droplet area A and the total solute C tot are then and the dimensionless form of the external force balance, Eq. (A4), is

APPENDIX C: LINEAR STABILITY ANALYSIS
Here we expand on our linear stability analysis of the circular homogeneous rest state discussed in Sec. III.

Special cases
As a series of sanity checks, we analyze Eqs. (15)- (17) in the special cases: a = 0, χ = 0, and σ = 0. First, it is clear that if the droplet-solute coupling is not closed (meaning aχ = 0), the roots of Eq. (17) are the eigenvalues that characterize the completely uncoupled system, where the eigenvalue s mσ 0 depicts the dispersion relation for shape perturbations of the passive droplet, and s m,n 0 are the infinitely many eigenvalues characterizing the decaying diffusion modes on a rigid disk with no-flux condition. More specifically, λ m,n is the nth real-positive root of J m (λ) = We stress that even if aχ = 0, our linear droplet-solute system may still be coupled via χ = 0 or a = 0. To understand why the growth rates in Eqs. (C1) and (C2) remain unaffected by any "one-way" coupling, we examine the eigenmodes associated with s mσ , s m,n in the two special cases: a = 0 and χ = 0.
No solute adsorption (a = 0). Solving Eqs. (15) and ( Let us first make sense of v mσ , Eq. (C3). Without adsorption, the solute is advected at the same speed of the fluid. This means that the moving boundary does not induce solute gradients; see Eq. (16). If the solute profile is not perturbed (δc = 0), the active force on the boundary is kept uniform. It follows that the pure shape perturbation is stabilized by the tension alone [i.e., with s mσ , Eq. (C1)].
The modes v m,n , Eq. (C4), can be understood as follows. Being that a = 0, the solute diffusion problem is unaffected by any small motion of the boundary. Thus, one expects to recover the classical diffusion modes on a rigid no-flux disk [associated with s m,n , Eq. (C2)]. The enslaved shape component in v m,n comes from the fact that small deviations in solute still produce a nonuniform active force on the boundary (via χ = 0). Note that in the special scenarios where σ m(m 2 − 1) λ 2 m,n , the mode v m,n effectively converges to the mode v mσ , Eq. (C3).
No force-transduction (χ = 0). Solving Eqs. (15) and (16) Let us first interpret v mσ , Eq. (C5). Being that χ = 0, the active force on the boundary is kept uniform and independent of small solute deviations. Hence, the shape perturbation must decay with the classical growth rate s mσ , Eq. (C1). Assuming a = 0, the movement of the free boundary still feeds into the solute transport problem [via Eq. The eigenmodes v m,n , Eq. (C6), are easier to understand. Since χ = 0, the shape is unaffected by small solute deviations, Eq. (15). Hence, if the shape is not perturbed (δR = 0), one simply recovers the classical diffusion modes on a rigid no-flux disk [associated with s m,n , Eq. (C2)].
Zero surface tension (σ = 0). We consider the limit σ = 0 and focus only on the multipolar modes m 2. Note that the modes m = 0 and m = 1, which are independent of σ , are discussed separately.
We begin by substituting σ = 0 in Eq. (17), It is clear that s = 0 is a root of Eq. (C7) for any aχ . Substituting σ = 0 in Eqs. (15) and (16), we find that s = 0 is associated with the eigenmode, v zst m0 (r, θ ) = 1 0 cos(mθ ), which describes the pure deformation of the droplet. With no surface tension, and with δc = 0, there is no nonuniform force acting on the boundary that could possibly counter or amplify v zst m0 . Hence, this perturbation must be marginally stable (accordingly, s = 0).
Since G zst m (s) is highly nonlinear, we cannot find its additional roots analytically. As we are interested in instabilities, we expand Eq. (C7) about s = 0 and obtain This function has the additional real root s zst m 4m(m+1)(aχ−1) 2+m(1−aχ ) , which changes sign from negative to positive as aχ exceeds the critical value of 1. We stress that s zst m is a valid approximation of a true eigenvalue so long as it is is small, i.e., for aχ ≈ 1, in which case, we may write s zst m ≈ 2m(m + 1)(aχ − 1). The associated eigenmode is then given by where we expanded the Bessel functions for low s zst m . It is instructive to note that as aχ passes 1, and hence v zst ms m becomes unstable, the solute gradient component in the eigenmode also changes sign with respect to the shape component. This instability is essentially a multipolar analog of the polarization-translation bifurcation (in mode m = 1; see following section).

Generic coupled case
Mass modes (m = 0). Let us consider the azimuthally symmetric perturbations of droplet shape and solute concentration. Substituting m = 0 in Eq. (17) gives As expected, G 0 (s) is independent of both σ and aχ . It has only real-negative roots: s = 0 and s 0,n = − j 2 1,n , where j 1,n is the nth zero of J 1 (x) (the Bessel-J function of order 1).
Substituting m = 0 back in Eqs. (15) and (16), we find that the root s = 0 is associated with two nontrivial eigenmodes, where v R 00 represents a uniform change in the droplet radius and v c 00 represents a uniform elevation in the solute concentration. Due to the conservation of total fluid and total solute, these mass perturbations are both marginally stable (accordingly, s = 0).

022404-10
Polarization-translation modes (m = 1). Substituting m = 1 in Eq. (17) gives We find that s = 0 is a root of G 1 (s) for all parameter values. Substituting m = 1 in Eqs. (15) and (16), and noting that J 1 (0) = 0, one finds that s = 0 has a single nontrivial eigenmode This mode merely describes the infinitesimal translation of the droplet. At the linear level, the shape remains circular and c remains uniform, as in the circular homogeneous rest state. Due to translational invariance, this perturbation is marginally stable (accordingly, s = 0). Since G 1 (s) is highly nonlinear, we cannot find its additional roots analytically. As we are interested in instabilities, we expand Eq. (C14) about s = 0: This function has the additional real root s 1 = 8(aχ−1) 3−aχ , which changes sign from negative to positive as aχ exceeds 1. Note that s 1 approximates a true eigenvalue of G 1 (s), Eq. (C14), so long as it is small, i.e., for aχ ≈ 1. In which case, we may write s 1 ≈ 4(aχ − 1). The associated eigenmode is then given by where we expanded the Bessel functions for small s 1 . It is instructive to note that as aχ passes 1, and v 1s 1 becomes unstable, the solute gradient component in v 1s 1 (r, θ ) also changes sign.
In the stable case (aχ < 1), the solute diffusion dominates over its relative drift with respect to the fluid. As the solute spreads out over the droplet domain (δc → 0), the polar driving force also vanishes, and hence the droplet slows down and stalls (u → 0). Given an initial gradient δc = r cos θ = x, the stopping distance of the droplet, x, can be computed from the eigenmode, Eq. (C17), which gives x = sign[s 1 ] χ 4(1−aχ ) = − χ 4(1−aχ ) . In the unstable case (aχ > 1), the solute drift with respect to the fluid dominates over diffusion. As the solute accumulates at the rear, it amplifies the forward droplet velocity, which then strengthens the solute asymmetry. Hence, in this case, v 1s 1 , Eq. (C17), represents an unstable perturbation that breaks the front-rear symmetry and leads to motility.
Multipolar modes (m 2). We now examine the roots of G m (s) for m 2 with aχ > 0 and σ > 0. It is easy to see that s = 0 is always a root of Eq. (17). However, taking m 2 and σ > 0, it follows from Eq. (15) that α m0 = 0 (the shape component vanishes), and since J m (0) = 0 for any m 1, the solute component also vanishes. Hence, for each m 2 with finite tension, s = 0 has only the trivial eigenmode (v = 0).
As we are interested in instabilities, we look for an eigenvalue s * [complex root of G m (s)] whose real part changes sign from negative to positive as a function of the destabilizing control parameter aχ . At the critical point, denoted aχ c , the real part of s * is zero and thus s * = iω (where ω ∈ R represents an oscillation frequency of the associated eigenmode). By solving Re[G k (iω)] = Im[G k (iω)] = 0, we may compute both ω and aχ c as functions of m and σ . Assuming small (yet nonzero) frequency, we expand G m (iω), Eq. (17), to up to order (m/2 + 3) in ω and obtain where F m (ω 2 ), H m (ω 2 ) are the following real-valued functions: (C20) Since s * = iω is a root of G m (s), it follows that F m (ω 2 ) 0 and H m (ω 2 ) 0. This system of two implicit equations can be solved in the variables ω 2 and aχ . Using Mathematica, we find one explicit solution that agrees with our result for the zerosurface-tension limit, i.e., for σ → 0 one obtains aχ c → 1 and ω → 0. We stress that this explicit solution can be considered a valid approximation of the bifurcation point so long as ω is small. Hence, we expand it here up to leading order in σ To recapitulate, we find a pair of complex-conjugate eigenvalues whose real part crosses 0 as aχ exceeds the critical point aχ c [approximated by Eq. (C21)]. At this Hopf-bifurcation point, the frequency ω is approximated by Eq. (C22). Using Eq. (15), we can express the coupled eigenmode associated with s * = iω. This eigenmode represents a shape-concentration standing wave, where we expanded the Bessel functions for small ω, and substituted Eq. (C22). Using Eq. (C23), we find that φ ≈ π ± arctan { √ 2/[σ (m − 1)]} approximates the phase shift between the shape and solute components in the standing wave. For low σ , we obtain φ ≈ ∓π/2.
The oscillation period is then where τ droplet = R 3 0 /(Mσ ) is the characteristic shape relaxation time of the viscous droplet. The fact that T is proportional to the geometric mean of τ solute and τ droplet implies that the oscillatory mechanism crucially depends on both relaxation processes.

Computational eigenvalues
We complete our linear stability analysis by finding the numerical roots of G m (s), Eq. (17). First, we trace these roots as functions of σ and aχ through a continuous extension procedure. Second, we describe our continuation method for tracing the Hopf-bifurcation lines over the aχ -σ diagram.
Continuous extension of the eigenvalues. For any m 1, one can compute the real-negative eigenvalues corresponding to aχ = 0, see Eqs. (C1) and (C2) and Table I. For any arbitrary σ , we increase aχ incrementally (starting from 0) and look for the numerical roots of G m (s) using the Find-Root function of Mathematica. This function implements a computational root-finding algorithm about an initial guess (s g ∈ C). To facilitate convergence, and to prevent inadvertent confusion between the multiple roots of G m (s), we define s g as the numerical root of interest obtained in the previous aχ iteration.
In Fig. 5 we demonstrate this analysis for m = 1-5. As there are infinitely many roots for each m, we focus on the droplet eigenvalue [s mσ = −σ m(m 2 − 1)] and the two leastnegative diffusion eigenvalues (s m,1 and s m,2 , given in Table I).
These roots are plotted as functions of aχ for different values of σ . The graphs essentially validate our main conclusions regarding the linear stability of each mode m, inferred previously through Taylor expansions of G m (s). Namely, the plots demonstrate that (1) there exists a single steady bifurcation in m = 1 (occurring at aχ = 1), (2) for σ = 0, all modes m 2 are closely analogous to m = 1, and (3) for σ > 0 (including large values of σ ), each m 2 is destabilized once through a Hopf bifurcation. Interestingly, the plots also reveal that the eigenvalue s mσ tends to merge with its closest neighboring s m,n as aχ is increased from zero. Upon merger, these two eigenvalues become a complex-conjugate pair, signifying a transition from a stable node to a stable focus. At high aχ , beyond the Hopf bifurcation, there is an additional point at which the complex conjugates separate into two distinct realpositive eigenvalues, signifying a transition from an unstable focus to an unstable node. We stress that these transitions do not imply an overall change in stability. We find computationally that all diffusion eigenvalues that do not merge with s mσ remain real-negative for any aχ .
Tracing the critical lines on the stability phase diagram.-For each m 2, we look for the critical line that marks the Hopf-bifurcation over the aχ -σ plane. Precisely on the bifurcation point, G m (s) has two complex-conjugate roots whose real part crosses zero. Hence, one must solve two implicit equations, Re[G m (iω)] = 0 and Im[G m (iω)] = 0, in the variable ω ∈ R and one of the control parameters (either aχ or σ ). Previously, we addressed this problem through a Taylor expansion of G m (iω), Eq. (C18). This allowed us to derive an explicit low-σ approximation for aχ c and ω at the critical point [see Eqs. (C21) and (C22)]. Here, we approach the same problem numerically, without expansions, using the FindRoot function of Mathematica directly on Eq. (17). This function requires a good initial guess of the two variables in order to converge on the solution of interest. We chose to work with small iterations of aχ (starting from aχ = 1) and solve Re[G m (iω)] = Im[G m (iω)] = 0 in the variables ω and σ . About aχ = 1, the critical tension σ and the frequency ω are both small, so we use Eqs. (C21) and (C22) to compute the initial guess: σ g = 4(m+2)(aχ−1) 3m 2 +m−4 and ω g = m(m + 1) √ 2σ g (m − 1). At higher values of aχ , the guess is computed by polynomial continuation of the previously registered numerical solutions. We use this procedure to recover the critical lines in the linear-stability phase diagram (Fig. 2, full colored lines). As expected, these lines are tangent to a a a a  Table I. In each plot, full lines mark the real part of the eigenvalues, while dashed lines mark the imaginary part. Bifurcations (dark diamonds) correspond to a critical aχ at which the real part of at least one eigenvalue turns from negative to positive. For m = 1 (and also for m 2 if σ = 0) there exists one purely real eigenvalue that changes sign precisely at aχ = 1. For m 2 (and σ > 0) there exist two complex-conjugate eigenvalues whose real part changes sign at a Hopf-bifurcation (marked by "H"). For each m, we find numerically that the infinitely many diffusion eigenvalues that do not merge with s mσ remain real-negative for any aχ . our explicit low-σ approximation, Eq. (C21) (thin translucent lines in Fig. 2), at the critical point (aχ, σ ) = (1, 0). . In all panels, the color density map in the bulk represents δc (negative in blue to positive in light green). The color-coded boundary represents δκ (negative in dark red to positive in yellow). The blue vector field in the bulk represent the instantaneous fluid flow, u = −∇δ p. On the boundary of the traveling wave (c), we also trace several fluid path lines over time.

APPENDIX D: STEADILY MOVING SOLUTIONS
Here we expand on our derivation of the self-consistent traveling-wave solutions discussed in Sec. IV.

Uniform flow from fixed shape assumption
The stationary shape assumption corresponds to vanishing normal velocity of the sharp interface in the moving frame of reference: where we used Eq. (B4). Lemma D.1. Let u satisfy Given Eq. (D1), u is uniform and given by where u cm is a constant in R 2 defined as u cm = A −1 ∂ (t ) x(u · n) dl. Proof. We define the flow in the moving frame of reference asũ = u − u cm . Since u cm is a constant and u is both incompressible and irrotational, it follows that ∇ ·ũ = 0 and ∇ ×ũ = 0. Thus, one can define as a Laplacian flow potential forũ:ũ Note that, for Hele-Shaw flow (Darcy's law), one has = −p − u cm · x.
From Green's identity, one has (D2) In Eq. (D2), it is easy to see that both the LHS and the boundary integral on the RHS vanish. Therefore, It follows thatũ = 0 in and thus u = u cm in (t ).

Pressure and solute profiles at steady state
We look for traveling-wave solutions, p(t, x, y) = p(x, y) and c(t, x, y) = c(x, y), wherex = x − u 1 t. Being that u 1 = −∇ p, it is straightforward that the pressure has the form p = p 1 − u 1x , where p 1 is a normalization constant for the droplet area. As for the solute, we substitute c = c(x, y) in Eqs. (B5) and (B6) and obtain [n x (au 1 + ∂ x ) + n y ∂ y ]c(x, y) = 0 on ∂ , where we used the fact that ∂ t c(x, y) = −u 1 ∂ x c(x, y). Lemma D.2.
For any smooth open subset ⊂ R 2 (with n = (n x , n y ) being the outward normal on ∂ ), the solution to the boundary-value problem in Eqs. (D4) and (D5) is unique and given in the following y-independent form: where c 1 is constant (normalization factor for the total solute).
Proof. Without loss of generality, let us substitute c(x, y) = f (x, y)e −au 1x back in Eqs. (D4) and (D5). One obtains the following boundary-value problem for f This elliptic Neumann problem is clearly solved by any constant f = c 1 . All that is left is to is show that nonconstant solutions do not exist. Assume that f ∈ C 2 ( ) ∩ C 0 (¯ ) is a nonconstant solution to Eq. (D7), we denote its maximum by N = max¯ f > 0. Being that the Laplacian is positive semidefinite, Eq. (D7) is uniformly elliptic, and thus the strong maximum principle holds. If there exists a point x ∈ such that f (x) = N, then by strong maximum principle f is constant. It follows that f < N in and there exists a point y on ∂ for which f (y) = N. This means that ∂ n f (y) > 0, in contradiction with the Neumann condition, Eq. (D8).

Rigid circular cell
We first work under the crude assumption that the steadily moving shape is the rigid unit disk, as in the rest state. This simplifying assumption is physically valid at the limit of high tension and low steady-state velocity, σ /(χ a 2 u 2 1 ) → ∞. Since the unit disk is physically viable and dimensionless, there is no need to explicitly resolve the constant p 1 .
We find the constant c 1 analytically (in terms of u 1 ) by imposing the dimensionless normalization of solute; see Eq. (B8). We consider the unit disk geometry and substitute the steady-state form of the concentration, Eq. (D6), in polar coordinates, c(r, θ ) = c 1 e −au 1 r cos θ , where I 1 denotes the modified Bessel function of the first kind of order 1.
To find u 1 , we substitute u cm = u 1x and c(r, θ ) back in the global external force balance, Eq. (B9). Taking thex component (n x = cos(θ )), we obtain We briefly interpret the significance of each term in this expansion. The sign of the first term controls the stability of the trivial solution u 1 = 0 (i.e., the circular homogeneous rest state). Clearly, this term is negative for aχ < 1 and positive for aχ > 1, in agreement with our linear stability analysis. The second term determines the second-order nature of the bifurcation into traveling solutions. This term changes sign at c s = 0.5. The saturation of the Hill-type force, Eq. (20), is expected to ultimately dampen the steady-state velocity at the nonlinear level. This is captured by our expansion so long as the last term is negative (that is, for 0.091273 < c s < 0.927885). Outside of this parametric range, Eq. (D10) is best solved numerically.
We obtain three branches of symmetric solutions to Eq. (D11), where g = c s [c s (2c s + 57) − 60] + 5 (strictly negative in the regimes of interest). The traveling solution branches, u 1s± and u 1u± , are realvalued over the parametric regimes, 1 and 2 , respectively. These are given by  bifurcation. In the supercritical case (c s 0.5), the stable moving states u 1s± branch out continuously from the stable rest state as aχ exceeds 1. In the subcritical case (c s < 0.5), the unstable moving states u 1u± branch out continuously from the unstable rest state as aχ falls behind 1. The latter scenario typically implies entrance into a bistable parametric regime. Indeed, over 2 , Eq. (D16), we find coexistence of both the stable rest state u 1 = 0 and the stable traveling states u 1s± . In this regime, |u 1s± | |u 1u± | and u 1s± = u 1u± = ± 2(c s +1) a √ 6(2c s − 1)/g on a saddle-node bifurcation, occurring at aχ * = 1 − The explicit Eqs. (D13) and (D14) provide useful analytical insights. However, we stress that these are strictly low-speed approximations of the true solutions to the exact Eq. (D10) (the steady-state force balance on the rigid cell).
Computational solutions. To complete the nonlinear bifurcation picture associated with the rigid problem, we solve Eq. (D10) numerically. As stated below this equation, we may set a = 1 without loss of generality. The idea is to trace the stable and unstable solution branches starting from the bifurcation point at χ = 1. To do so, we employ a continuousextension procedure in which we vary χ incrementally. Working with Mathematica, we use (1) the NIntegrate function to compute the integral in Eq. (D10), and (2) the FindRoot function to obtain a numerical root of this equation about an initial guess u g 1 . To facilitate convergence and to ensure continuity of the solution branch, we define u g 1 as the numerical root obtained in the preceding parametric iteration.
Remark D.1. About the bifurcation point (χ = 1), the true u 1 is small, and so we use Eqs. (D13) and (D14) to define the initial guess. Specifically, we take u g = u 1s+ for c s 0.5 and u g = u 1u+ for c s < 0.5.
In Fig. 6 (left panels), we compare the results of our numerical continuous-extension procedure (thick lines) to the approximations given in Eqs. (D13) and (D14) (thin lines). The numerical bifurcation structure behaves as expected in both the low and high c s regimes, including ranges over which the analytical approximations, Eqs. (D13) and (D14), are no longer valid. Due to the transition from a super-to a subcritical bifurcation at c s = 0.5, and the emergence of a saddle node in the subcritical regime, it is clear that Eq. (D10) admits three motility phases: (1) rest phase, where only the stable rest state exists (aχ < 1 and c s 0.5 or aχ < aχ * N and c s < 0.5, with aχ * N denoting the numerical saddle-node bifurcation point), (2) traveling phase, where the rest state is unstable and there exists a stable traveling state (aχ > 1), and (3) bistable phase, where both the rest state and a high-speed traveling state are stable (aχ * N < aχ < 1 and c s < 0.5). In this phase, we find an additional low-speed steady state that is unstable.
In Fig. 7 we show different representations of the motility phase diagram. Figure 7(a) is the χ -c s diagram (essentially a top view of the 3D bifurcation diagram in Fig. 6). Figure 7(b) shows the motility phases on the β-c s plane [the original parameters defining Eq. (20)]. Note here the close correspondence with the analogous 1D model [44] (with its deterministic version analyzed in the Supplemental Information of Ref. [45]). Since the pitchfork bifurcation occurs at β = (c s + 1) 2 /c s = 2[1 + cosh(log c s )], we found it useful to also present the phase diagrams in terms of log c s [Figs. 7(c) and 7(d)]. Interestingly, Fig. 7(d) reveals that β * N (the numerical saddle node bifurcation) is closely approximated by the tangent to the critical line at c s = 0.5.

Deformable cell
We wish to obtain the steadily moving states for a deformable cell of arbitrary surface tension σ . This problem consists of finding both the speed and the shape of the steady state in a self-consistent manner. It is therefore more complicated than the rigid circular case, which amounted to one nonlinear equation in one variable (u 1 ). To deal with this challenge, we must first derive three necessary and sufficient equations for finding the three numbers that make up our reverse-engineered solution (u 1 , p 1 , c 1 ). Shown are different parametric representations of the motility phases prescribed by Eq. (D10). In each plot, the thick black line marks the motility pitchfork bifurcation occurring at χ = βc s /(c s + 1) 2 = 1. To the right of this line, we find the traveling phase (dark blue region), in which the rest state is unstable and there exists a stable traveling state. The rest state is always stable to the left of this line. For c s < 0.5 and aχ * N < χ < 1 (bistable phase, light cyan region) there exists in addition a stable (high-speed) traveling state and an unstable (low-speed) traveling state. The leftmost limit of this phase (thin line with points) corresponds to the numerical saddle-node bifurcation. The rest phase (white region) spans the parametric regime over which there are no steadily moving solutions.
Working in the reference frame of the moving cell, our aim is to recover the shape defined by Eq. (21). Let us reproduce this equation here for convenience: Note that κ (x) depends on the three solution variables (u 1 , p 1 , c 1 ) and the four dimensionless parameters (σ , a, β, c s ). This can be simplified further by multiplying Eq. (D19) by a = 0. It is then easy to see that the same κ (x) depends on three variables (au 1 , ap 1 , c 1 ) and only three contracted parameters (aσ , aβ, c s , where aβ can be replaced by aχ (c s + 1) 2 /c s ). Hence, like the rigid case, we can set a = 1 without loss of generality.
The fact that the steady-state curvature is strictly a singlevalued function of x implies reflection symmetry about the polarity axisx. It is therefore natural to define the boundary line with the curve ±y(x), where +y(x) 0 represents the top half of ∂ and −y(x) represents the mirrored bottom half. The curvature on +y(x) is given by κ = −y (1+y 2 ) 3/2 . Changing variables to Y = −σ y / 1 + y 2 , such that Y = σ κ, reduces Eq. (D19) to a first-order ODE in Y (x). We solve this equation analytically along with the condition Y (0) = 0, which corresponds to y (0) = 0: Remark D.2. The condition y (0) = 0 is chosen arbitrarily for convenience. It essentially aligns the top and bottom poles (meaning the symmetric poles that are transverse to the direction of motion) at x = 0. While p 1 and c 1 vary due to shifts along the axis of motion, the translational invariance of the problem implies that the shape itself and the speed u 1 will be unaffected by this arbitrary choice.
Next, we find the leftmost and rightmost limits of the curve y(x). Strictly speaking, these are the singular points x L < 0 and x R > 0 at which y (x L ) = +∞ and y (x R ) = −∞. We may use Eq. (D20) to compute x L and x R numerically, substituting Y (x L ) = −σ and Y (x R ) = +σ .
Finally, we find the curve y(x) over the domain x ∈ (x L , x R ) by integrating numerically the first-order ODE , y(x R ) = 0.
Here the boundary condition y(x R ) = 0 ensures that the curves ±y(x) join continuously at the front end [the point (x R , 0)] without forming a cusp. In other words, this condition, along with the fact that y (x R ) = −∞, guaranties the smoothness of ∂ at the front. For ∂ to be physically viable, an equivalent condition must hold at the rear. The freedom to impose this additional boundary condition is implicit in our choice of (u 1 , p 1 , c 1 ). Three constraints. For a given set of dimensionless model parameters, we look for solutions to Eq. (D21) that satisfy the following three constraints: (1) The boundary defined by ±y(x) represents a simple (nonintersecting) closed contour without cusps. To exclude intersections, we must take care that y(x) 0 for x ∈ (x L , x R ). To unsure the smoothness of ∂ at the rear, we look for a solution that satisfies y(x L ) = 0.
(2) The dimensionless area is conserved, such that (3) The dimensionless solute is conserved, such that Together, these constraints describe a system of three implicit nonlinear equations in the three variables (u 1 , p 1 , c 1 ). If a dimensionless steadily moving solution exists for a given set of parameters, it can in principle be found by varying those three numbers while reintegrating Eq. (D21).
Optimization procedure. Rather than resorting to an exhaustive 3D grid search, we employ the following optimization procedure that greatly facilitates convergence.
(1) For a new set of dimensionless model parameters, we begin by finding a good quantitative guess of the solution (u g 1 , p g 1 , c g 1 ). Normally, the computation of (u g 1 , p g 1 , c g 1 ) relies on polynomial interpolations of (u 1 , p 1 , c 1 ) as functions of the varied parameter. Technically, these interpolation functions are constructed using the registered numerical solutions obtained in preceding parametric iterations.
(2) We define a matrix of phase-space coordinates {(p 1 , c 1 )} about (p g 1 , c g 1 ). For each such coordinate, we run the following computational element.
(3) With fixed (p 1 , c 1 ), we resolve constraint no. 1 (minimization of |y(x L )|) by varying u 1 about u g 1 while reintegrating Eq. (D21). This automated shooting-solution method converges rapidly on a velocity u 1 > 0 that grants an acceptable error of |y(x L )| < 10 −7 . The resultant contours ±y(x) then bound a shape that is physically viable (∂ is closed without cusps), but generally not consistent with the dimensionless normalization conditions: constraints no. 2 and no. 3. We compute the area A and the total solute C tot associated with the physically viable solution.
(4) We use the registered numerical data from steps 2 and 3 to generate the polynomial interpolation functions, A(p 1 , c 1 ) and C tot (p 1 , c 1 ). We then solve two implicit functions in two variables, namely, A(p 1 , c 1 ) − π = 0 and C tot (p 1 , c 1 ) − π = 0. The solution, (p * 1 , c * 1 ), is then substituted back in element 3, giving a speed u * 1 > 0 and a corresponding curve y * (x) that resolve constraint no. 1.
By using this powerful continuation + minimization procedure we are able to converge with high precision on the coveted steady states. Difficulties in convergence arise in highly nonlinear regimes (e.g., close to bifurcations), where the speed and/or the shape vary strongly as functions of the parameters or as functions of (p 1 , c 1 ). Hence, in these regimes we work with smaller parametric iterations and denser {(p 1 , c 1 )} grids. While these practices are generally time-consuming, they tend to improve the accuracy of both the initial guess and the interpolation functions, which together promote convergence. Remark D. 3. As we span the parameter space, we fix σ and c s arbitrarily and perform incremental iterations of χ [substituting β = χ (c s + 1) 2 /c s in Eq. (D20)]. We always begin at χ = 1 (the motility-mode instability), where the steadily moving states are known to branch out from the circular rest state via a super-or subcritical pitchfork bifurcation. About this critical point, both the speed u 1 and the shape deviation from the circle are small, so we can use our rigid-droplet results to compute the initial guess. In more detail, Eqs. (D13) and (D14) are used for computing u g 1 . Then u g 1 is substituted in Eq. (D9) for computing c g 1 . Finally, c g 1 is used to compute p g 1 = σ − F act (c g 1 ). Remark D.4. Any physically viable solution (u 1 , p 1 , c 1 ), which resolves constraint no. 1 but is not consistent with constraints no. 2 and no. 3 (obtained, e.g., via step 3 in our procedure), can be mapped directly to a self-consistent dimensionless solution. To perform this mapping correctly, one must rescale the solution variables and the model parameters in accordance with our nondimensionalization scheme, where R 0 = √ A/π . Note that when working with χ instead of β, i.e., with β = χ (c s + 1) 2 /c s , one should rescale χ like χ → C tot (c s +1) 2 A(c s +C tot /A) 2 χ . For the new set of three dimensionless parameters, the rescaled solution automatically satisfies all three constraints (discussed above) and is therefore a dimensionless steadily moving state. In this appealing shortcut to constraints no. 2 and no. 3, we stress that the dimensionless parameters (σ and c s specifically) may not be fixed a priori. Hence, unlike our complete optimization procedure, this direct rescaling approach is not well-suited for spanning the parameter space in a controlled, structured manner.