Modes of mechanical guidance of adhesion-independent cell migration

Abstract

Adhesion-independent migration is a prominent mode of cell motility in confined environments, yet the physical principles that guide such movement remain incompletely understood. We present a phase-field model for simulating the motility of deformable, non-adherent cells driven by contractile surface instabilities of the cell cortex. This model couples surface and bulk hydrodynamics, accommodates large shape deformations and incorporates a diffusible contraction-generating molecule (myosin) that drives cortical flows. These capabilities enable a systematic exploration of how mechanical cues direct cell polarization and migration. We first demonstrate that spontaneous symmetry breaking of cortical activity can lead to persistent and directed movement in channels. We then investigate how various physical cues – including gradients in friction, viscosity, and channel width as well as external flows and hydrodynamic interactions between cells – steer migration. Our results show that active surface dynamics can generate stimulus-specific cell behaviors, such as migration up friction gradients or escape from narrow regions. Beyond cell migration, the model offers a versatile platform for exploring the mechanics of active surfaces in biological systems.

---

1 Introduction

Adhesion-independent migration, commonly referred to as amoeboid migration, plays a key role in development, immune surveillance, and cancer invasion.^1–7 In vitro studies have shown that this migratory mode can be induced in a wide range of cell types when subjected to three-dimensional confinement.^6–8 Although there are alternative mechanisms,⁹ amoeboid motility is often driven by actomyosin accumulation at the cell rear, which produces retrograde flows of the cell cortex. This cortical flow then propels the cell through friction with the environment.^3,5–7,10 Although several theoretical frameworks have been proposed to explain how this migration is initiated and sustained, the exact processes governing the symmetry breaking that triggers polarization remain incompletely understood.^{1,2,6,8,10,11}

Minimal models^{2,3,6,7,10,12–19} have suggested that the cell cortex can be described as an active fluid surface that self-organizes tension-generating molecules – such as myosin – into a localized concentration peak through the interplay of advective transport, contractility, and surface flows. If myosin accumulates towards one side of the cell, this side becomes the cell rear and it initiates retrograde flows that propel the cell forward, consistent with amoeboid migration dynamics.²⁰ While most prior studies have focused on such self-organization on fixed geometries,^18,21–24 experiments and theory have shown that shape changes can actively feed back on cortical flows,^7,25 underscoring the need to incorporate deformable surfaces into modeling frameworks.

Deforming active surfaces were considered in ref. 17, 19, 26 and 27, coupled to hydrodynamics and material flow. Significant shape deformations including strong constrictions could be reproduced for tubular surfaces¹⁷ as well as for ellipsoidal and spherical surfaces.^19,27 However, all previous methods operate with a grid-based representation of the surface. Correspondingly, the surrounding medium was either neglected^17,26 or limited to a simple homogeneous fluid without walls or obstacles.^19,27 Recent work has addressed some of these limitations by modeling the effects of viscous and viscoelastic resistance in confined migration scenarios, providing new predictions on polarization thresholds.¹⁰ In general, the grid-based Lagrangian surface representation of these approaches makes it difficult to include (wall) contact as the mesh typically deforms and tangles in the contact zone and additional contact detection and repulsion algorithms are required. Also, grid-based methods are not suited for topological transitions like cell division, which is often modeled as a self-organized process on an active surface.^18,27

To overcome these limitations, we propose a phase-field approach to represent the evolving surface. The phase-field description makes the model well-suited for simulating migration through complex environments. The model is not only flexible to deal with complex surrounding geometries, including contact to walls and obstacles, but also enables the simulation of topological cell-shape changes, as encountered during self-organized cell division. We then use this model to explore how polarization and migration of cells are triggered by various mechanical cues, such as gradients in friction or viscosity, surrounding fluid flow, and constrictions. Altogether, we establish a versatile modeling platform for studying active surface dynamics and use it to analyze how cells choose direction in diverse microenvironments.

2 Active surface model

To develop a diffuse-interface model, we first introduce a sharp-interface formulation describing the dynamics of an active gel surface immersed in a fluid medium. To this end, we recast the governing equations given in ref. 17–19 in a form suitable for a later conversion into a diffuse-interface representation. Afterward, these equations are non-dimensionalized and the diffuse-interface formulation is introduced, which finally allows us to handle wall contact such as to simulate pattern formation and cell migration in a channel.

2.1 Sharp-interface model

The spatial domain, denoted by Ω, is divided into the intracellular fluid domain Ω₁ and the exterior fluid medium Ω₀. Both domains are separated by the cell surface Γ, representing the cell membrane and cortex, see Fig. 1 (top).

Fig. 1 Simulation setup of a cell in a cylindrical channel. Top: Sharp-interface model with cell (Ω₁) and fluid (Ω₀) separated by interface Γ. Due to axisymmetry, only the blue rectangle is computed. Bottom: Diffuse-interface representation of the inset region by a phase field ϕ (left) and the corresponding adaptive mesh (right).

We refer to the normal vector pointing from Ω₀ into Ω₁ as n, and use it to define the projection P_Γ onto Γ as P_Γ := I − n ⊗ n, where I denotes the identity matrix. Differential operators on the surface Γ are defined as follows:

∇_Γ := P_Γ∇ surface gradient,

∇_Γ := P_Γ: ∇ surface divergence,

Δ_Γ := ∇_Γ·∇_Γ Laplace–Beltrami operator.

Note that applying these operators to a field variable requires the latter to be defined not only on Γ but in a neighborhood of Γ. Throughout this article, all surface variables are assumed to be defined in this way.

We describe intracellular and extracellular flows in the whole domain Ω by the Stokes equations

∇·[η(∇u + ∇u^T)] − ∇p = f_σ + ∇_Γ·(δ_ΓS_visc),

∇·u = 0,

where η(x) denotes the fluid viscosity, both inside and outside the cell: η(x) = η_i for x ∈ Ω_i, i = 0, 1. Moreover, f_σ and S_visc denote the active surface tension and surface viscous stress, specified below, and δ_Γ is a surface Dirac-delta distribution, which is defined as δ_Γ(x) = 0 if x ∉ Γ and , where |Γ| is the area of any surface Γ⊂Ω. We model the fluid dynamics using a single, continuous velocity field u. The cell surface thus moves with this flow field. As a result, there is neither tangential slip nor fluid flow across the interface, leading to an intrinsically imposed no-slip condition on both sides of the cell surface.

Any flow and shape deformation along the cell surface is limited by the ability of the cortex to deform and remodel itself, such that accounting for the cortex mechanics is indispensable.^28,29 On the relevant time scales of cell division and cell migration, the rheology of the cortex is predominantly viscous,²⁹ which gives rise to a viscous surface stress S_visc, as introduced by Scriven³⁰ as

S_visc = (η_b − η_s)(∇_Γ·u)P_Γ + η_sP_Γ(∇_Γu + ∇_Γu^T)P_Γ,

where η_b and η_s are the bulk and shear viscosity of the surface.

The surface tension force f_σ is produced by force-generating molecules such as myosin motor proteins. As in previous literature,^17–19 we model the tension as a monotonically increasing Hill-function of the myosin concentration c by , where c₀ is the characteristic concentration and σ₀ > 0 a scaling factor. Consequently, the surface tension force can be formulated as

f_σ = ∇_Γ·(δ_Γσ(c)P_Γ) = δ_ΓHσ(c)n + δ_Γ∇_Γσ(c), (1)

where H = ∇·n is the total curvature of the surface. The first term on the right-hand side corresponds to the normal component of surface tension. Respectively, the second term is the tangential component, known as the Marangoni force, which arises from gradients of surface tension along the surface. It provides a force toward regions of high myosin concentrations, and is therefore responsible for the retrograde flows that drive amoeboid migration.

These equations are coupled to an advection-diffusion equation to describe the evolution of the myosin concentration c, on Γ:

where d/dt is the material derivative and D > 0 the diffusion constant. We omit more complex formulations of eqn (2) that incorporate detailed binding/unbinding kinetics.^14,31,32 We chose this simplification because these specific kinetics are not necessary to trigger the symmetry-breaking instability that triggers cell migration. Therefore, we use the simplest possible model, although these terms can be readily incorporated into the present framework.

2.2 Non-dimensionalization

To reduce the number of model parameters, we non-dimensionalize the equations by rescaling length in units of the radius R of the cylindrical channel, time in units of the diffusive timescale D/R², and concentration in units of the characteristic concentration c₀ in the Hill function for the active stress. Then, dividing eqn (2) by c₀R²/D yields the dimensionless concentration equation

For the hydrodynamic equations, we base the scaling on the surface bulk viscosity η_b. Thus, we rescale pressure by Dη_b/R³ and introduce the surface viscosity ratio ν = η_s/η_b and the Péclet number as the ratio between time scales of surface diffusion (R²/D) and active flows (η_b/σ₀). We obtain

∇·u = 0, (5)

with . The dynamics of the system is thus governed by the three parameters Pe, ν and ηR/η_b. The values of the parameters used in our simulations are provided in Table 1 in the results section.

Table 1 Default model and simulation parameters (top) and experimental estimates of physical parameters (bottom)

Time step	Δt	0.01
Bulk grid size	hbulk	2.5 × 10^−1
Interface grid size	hint	5.52 × 10^−3
Mobility	M	10^−3
Interface thickness	ε	0.02
Volume correction	β	0.01
Mass correction	α	0.25
Surface viscosity ratio	ν = ηs/ηb	1
Fluid viscosity	η0R/ηb = η1R/ηb	1
Péclet number	Pe	25
Myosin diffusivity^45	D	100 µm^2 min^−1
Viscosities
– Cortical^29,46	ηb,ηs	1 mN s m^−1
– Cytoplasmic^47,48	η1	10^−2–10^1 Pa s
	η1R/ηb	10^−4–10^−1
– Exterior^49–53	η0	10^−3–10^6 Pa s
	η0R/ηb	10^−5–10^4
Cortical tension^29,54	σ0	1 mN m^−1
Cell/channel radius	R	10 µm

---

2.3 Phase-field Ansatz

To flexibly account for large deformations, topological transitions and wall contact, we introduce a phase-field (i.e. diffuse-interface) version of the above equations. In the phase-field model, an auxiliary field variable ϕ – the phase field – is introduced and used to indicate the bulk phases, which can be arbitrary viscous fluids,³³ as well as viscoelastic or elastic materials.³⁴ The phase-field function varies smoothly between the distinct values for the phases across the interface, resulting in a small but finite interface thickness ε. Depending on the application of interest, phase-field methods may offer advantages over other interface-capturing methods. For example, they allow for unconditionally stable inclusion of surface tension³⁵ and fully-discrete energy-stable schemes, see e.g. ref. 36 and 37.

Here, we describe the geometry of the cell such that ϕ(x) ≈ 1 in Ω₁ and ϕ(x) ≈ 0 in Ω₀ with a smooth transition between the two phases, see Fig. 1 (bottom). The phase field is initialized by , where r is the signed distance to Γ, positive in Ω₁. For complex interface geometries, r can be computed numerically using redistancing algorithms.³⁸ In our simulations, we employ a simple cigar shape (cylinder with two hemispherical caps, Fig. 1) which permits an explicit analytical expression for r.

After initialization, the phase field is advected with the fluid flow to capture changes in the surface geometry. This evolution is governed by the convective Cahn–Hilliard equation, which is solved over the entire computational domain Ω:

μ + εΔϕ − ε⁻¹W′(ϕ) = 0, (7)

where μ denotes the chemical potential and W(ϕ) = ϕ²(1 − ϕ)² a double-well potential. The parameter ε > 0 describes the thickness of the interface region, and the mobility M > 0 governs the thermodynamically driven diffusion of the phase field. In our approach, the diffuse interface merely approximates a sharp cell surface, and the phase-field framework acts only to advect that interface consistently. Accordingly, interface motion is imposed by the velocity field rather than diffusion. We therefore choose M ≪ 1, so that it primarily stabilizes and smooths the interface profile without shifting its position. The term scaled with constant β > 0 serves as a volume-correction mechanism to relax the cell volume back to its initial value V(0). In our simulations β is chosen large enough such that volume loss is limited to roughly 1%. Note that the correction term compensates for the fact that the geometric volume enclosed by the ϕ = 0.5 contour is not conserved under the Cahn–Hilliard dynamics, ensuring a more accurate control over the domain occupied by the cell.

The phase-field representation can not only be used to describe the surface Γ = {x|ϕ(x) = 0.5}, but also to approximate the Dirac-delta distribution δ_Γ ≈ |∇ϕ|. Following the diffuse-interface approach,³⁹ the concentration eqn (3) can be extended from the submanifold Γ onto Ω as

∂_t(|∇ϕ|c) + ∇·(|∇ϕ|cu) − ∇·(|∇ϕ|∇c) = 0 in Ω. (8)

See Appendix A.1 for a derivation. Even though the concentration equation suggests mass conservation, numerical discretization errors may lead to small deviations, which can accumulate to a significant mass loss over long simulation times. Therefore, we introduce the mass on the surface

and add a mass correction term
with some constant α > 0 to the right-hand side of the concentration equation to ensure mass conservation.

Finally, we reformulate the momentum eqn (4) in the phase-field formalism. The fluid viscosity is linearly interpolated between the distinct viscosities in the phases,

As usual for diffuse-interface models of two-phase flow (e.g. ref. 35), the constant surface tension term δ_ΓHσn can be reformulated to , where the numerical factor stems from the chosen double-well potential. Moreover, we use the extended normal vector ñ := ∇ϕ/|∇ϕ| ≈ n and the surface projection P̃_Γ := I − ñ ⊗ ñ ≈ P_Γ to define diffuse-interface approximations to the surface differential operators, e.g. the surface gradient and surface divergence . We obtain a diffuse-interface version of the surface stresses and forces

which enter the diffuse-interface hydrodynamics equation

−∇·[η(ϕ)(∇u + ∇u^T)] + ∇p = ∇·(|∇ϕ|S_visc) + f_σ (10)

in Ω.

Note that the first divergence operator on the right-hand side of eqn (10) is not required to be a surface divergence, since it is applied to a tangential tensor.

2.4 Velocity extension

As we will see in the numerical tests, the active surface tension imposes strong tangential flows leading to regions of large tangential compression (∇_Γ·u < 0) or stretching (∇_Γ·u > 0). Due to incompressibility of the flow field, this goes along with the opposite deformation in the normal direction, i.e. tangentially stretched regions get compressed in the normal direction and vice versa. In numerical tests we find that such strong compressional or extensional flows in the normal direction tend to locally shrink or expand the thickness of the interface region, respectively, unless an extremely high, unphysical mobility is used. To eliminate this perturbation of the interface profile for reasonable mobilities, we instead advect the phase field with an auxiliary velocity field v. The idea is that v is constructed as an extension of the velocity u evaluated at the surface (i.e. at the 0.5-level set), which is constant in the direction n normal to the surface. Therefore, using v instead of u in the advective term, velocity differences across the diffuse interface are leveled out, which leads to a consistent advection of the complete diffuse interface region. We propose to construct v by solving the additional equation for its components (v₁, v₂, v₃),

|∇ϕ|v_i − |∇ϕ|u_i − ∇·[|∇ϕ|ñ ñ·∇v_i] = 0 (11)

on Ω and for i = 1, 2, 3.

In the Appendix A.2 we show by matched asymptotic expansion that this formulation converges to the following sharp-interface limit equations

v = u on Γ, (12)

∇v·n = 0 on Γ. (13)

The obtained velocity field v is used to replace the velocity u in the advection terms of both the Cahn–Hilliard equation and the concentration equation.

2.5 Discretization

To avoid numerical instabilities from capillary time step constraints, the six strongly coupled eqn (5)–(8), (10) and (11) are discretized by a monolithic linear semi-implicit time discretization.³⁵ The proposed time discretization ensures mass conservation of regulating surface species on the discrete level as shown in Appendix A.3.

We solve the resulting system in each time step with a finite element method based on the finite element toolboxes DUNE⁴⁰ and AMDiS.^41–43 An adaptive grid is employed to accurately resolve the phase field and surface forces with more details given in Appendix A.4. To make simulations more efficient, we assume that the cell shape and concentration field are axisymmetric. This assumption, which holds in particular for homogeneous cells in cylindrical channels, reduces computations effectively to a 2D domain from which the full 3D solution can be recovered, see Fig. 1. Further details on the axisymmetric formulation can be found in Appendix A.5.

3 Results

To examine how cells determine their direction in confined geometries – such as microchannels or microcapillaries – we consider a cigar-shaped cell initially positioned inside a cylindrical channel (see Fig. 1 (top)). The cell is initialized with a length of 3.24 and a height of 1.84, placed within a channel of diameter 2. The channel length varies and is chosen sufficiently large to ensure that the cell is unaffected by the channel openings. The interface thickness is selected to be sufficiently small such that the specific choice of ε has little influence on the results; see Appendix A.6 for a more detailed discussion. A nearly homogeneous cortical concentration field is prescribed with a base value of 1 and superimposed uniform noise centered on 0 and a width of 0.1. The channel walls are equipped with no-slip (v = 0) and no-adhesion (ϕ = 0, n·∇ϕ = 0) boundary conditions. Unless otherwise specified the model and simulation parameters are as given in Table 1.

The chosen Péclet number corresponds to active cortical tension σ₀ on the order of 1 mN m⁻¹ as suggested from measurements in cells.⁴⁴ The critical Péclet number at which pattern formation occurs can be very roughly estimated from the linear stability theory in ref. 18, which however neglects the exterior fluid and channel walls. Based on our simulations, we identify a critical Péclet number of Pe ≈ 20 (see Appendix A.7). Starting simulations from random initial concentrations, we consistently find the emergence of dominant polar patterns of mode 1 (i.e. high concentration at one cell end and low at the other).

The chosen ratio between cortical and fluid viscosity ηR/η_b = 1 was selected to specifically study a balanced regime where the viscous contributions of the surrounding fluid and the cell cortex are comparable, clearly exposing the coupling mechanism necessary for cellular propulsion. While physically realistic values of cytoplasmic viscosity (see Table 1) suggest a much lower ratio, our results are expected to be robust as in both regimes single-spot myosin concentration patterns (1-modes) emerge.¹⁸ To verify this, we performed additional representative simulations with a smaller fluid-to-cortical viscosity ratio (ηR/η_b = 0.25). These confirmed that the cell consistently develops the same single-spot polarity, and that the overall influence of external cues on the symmetry-breaking mechanism remains robust.

3.1 Spontaneous symmetry breaking

Studying the polarization process experimentally is challenging due to the small spatial scales and complex interplay of physical forces involved. Our modeling approach allows us to overcome these limitations and gain mechanistic insight into how polarization and migration emerge in non-adherent cells. To begin, we consider a single cell confined within a channel, without any external cues, in order to isolate the intrinsic mechanisms driving polarization and movement.

As shown in Fig. 2, the cell initially rounds up (i.e. reduction in length and increase in width) and adapts to the channel. Due to the properties of the phase-field model – specifically, the boundary condition ϕ = 0 at the channel walls and the resulting phase-field profile – a thin liquid film remains between the cell and the wall. The thickness of this film scales with the interface thickness ε. A systematic validation confirms that, for sufficiently small ε, the results are insensitive to its precise value (see Appendix A.6 for details).

Fig. 2 Cell migration via spontaneous symmetry breaking in a channel. (a) A cell is placed inside a channel and elongated relative to its relaxed shape. The surface concentration is initialized to be uniform with slight perturbations. (b) Snapshot of the simulation at t = 1.36. The initial contraction, driven by surface tension, leads to the formation of two spots of high myosin concentration at the cell ends. Streamlines show the velocity field. In case I, the concentration and velocity are slightly higher on the right end; in case II, higher on the left end. (c) Velocity at the cell ends over time. The blue dotted line corresponds to the left end and the green dashed line to the right end. Initially, the cell ends move towards each other as the cell contracts. However, the initial symmetry breaks spontaneously and the cell ends up migrating in one direction. (d) Concentration profiles along the cell surface at different times. Early asymmetries in contraction lead to a buildup of concentration at one end, ultimately defining the cell rear and the direction of motion. (e) Snapshot of the final state at t = 15. Over time, the asymmetries intensify, leading to full polarization and persistent migration to the left (case I) or right (case II). Note that flow velocity arrows are displayed in the lab frame. Even though the vectors visually cross the cell surface, the no-penetration condition is fulfilled, as also the cell boundary moves at the local fluid velocity.

As the cell contracts, two regions of high myosin concentration emerge at the poles. Although this state is nearly symmetric, small asymmetries – stemming from noise in the initial myosin concentration – gradually amplify over time through a self-reinforcing process driven by the Marangoni force. Because the myosin concentration at one end is slightly higher, there is a net flow of myosin towards that region. This net flow eventually causes the peak at one side to outcompete the other, pulling in all the available myosin and breaking the initial symmetry. This process of spontaneous symmetry breaking forms a single spot of high myosin concentration with equal probability at either end of the cell.

Eventually the cell fully polarizes, Fig. 2(d) and (e). The asymmetric distribution of myosin generates a constant force toward the high-concentration region, inducing a retrograde flow (in the cell's frame of reference) of the cell cortex. This flow interacts frictionally with the confining walls, through shear stresses transmitted by the thin liquid layer between the cell surface and the channel wall, allowing the cell to transmit tangential forces to its surrounding. As the rear contracts and the cortex flows backward, its friction against the channel walls propels the cell body forward. At the same time, the cytoplasm's incompressibility leads to a recirculating flow pattern within the cell, Fig. 2(e).

Note that the frictional interactions with the channel walls can arise through two distinct mechanisms. First, they may result from the close proximity with a no-slip wall, which inhibits tangential motion of the cell surface, either by direct contact or by strong viscous shear forces that penalize velocity differences between the cell and the near boundary. In contrast to these hydrodynamic effects, friction can also be modeled as an explicit force opposing fluid motion. This latter approach is employed to investigate the recently observed phenomenon of friction-guided cell migration in the next section.

3.2 Frictio- and viscotaxis: cells migrate up friction and viscosity gradients

Cells rarely migrate in a purely random manner; rather, their movement is often directed by external clues – for example, in response to nutrient sources or during the wound healing process. While migration guided by chemical gradients – known as chemotaxis – is well studied,^2,55–58 cell migration guided by mechanical cues remains relatively less explored. In particular for adhesion-independent migration in channels, recent work found that cells can follow gradients of friction – a phenomenon that was called frictiotaxis.⁵⁹ Here, we explore this new mode of directed migration with our active surface model. Moreover, we also consider guidance by gradients of viscosity of the surrounding medium.

Using in vitro microchannels, ref. 59 showed that cells tend to migrate persistently up a friction gradient – that is, toward regions of higher wall friction. The authors also proposed a simple one-dimensional model to explain how cells polarize towards higher frictions. Here we extend this analysis by investigating frictiotaxis through our 3D model, which accounts for changes in cell shape and captures the fluid flows not just in the cortex but also in the cytoplasm and the extracellular medium.

In our axisymmetric 3D setup, the friction gradient is incorporated by introducing an additional friction force f_fric = −χ_wallk(z) u·ẑ ẑ to the right-hand side of the force-balance eqn (4). This force is oriented in the axial direction ẑ opposing the axial velocity u·ẑ, scaled with a cell–wall friction coefficient k(z) which varies linearly across the channel ranging from k_left on its left to k_right on its right end. The force is localized by the characteristic function χ_wall, which equals 1 only within a thin friction layer of thickness 0.05 and 0 elsewhere. This ensures that the friction force is integrated solely over a narrow region adjacent to the channel wall.

As shown in Fig. 3 (top), the cell develops regions of high myosin concentration at both ends at early times. However, the symmetry in the flow field is broken by the presence of the friction gradient, see Fig. 3(a I). The lower friction on one side drives two key processes. First, it promotes faster contraction. Second, it accelerates the surface transport of myosin toward that same side. Both of these effects lead to a faster accumulation of myosin on the low-friction side. This process is self-reinforcing due to the Marangoni force, leading to the cell polarization. As a result, the cell migrates up the friction gradient driven by the resulting retrograde flows, Fig. 3(a II). A closer examination of the velocity at both cell ends reveals that this symmetry breaking occurs very early in the contraction phase (Fig. 3(b)). The cell speed continuously increases on longer time scales until a maximum velocity is reached, which then slowly decays due to the constantly increasing friction, Fig. 3(c). Note that this decay does not occur when friction is uniform (Fig. 3 (c, gray)).

Fig. 3 Guidance of cell migration by friction and viscosity gradients. A cell is initialized in a channel with either a friction gradient on the channel wall (top panel) or a viscosity gradient in the surrounding fluid (bottom panel). The internal viscosity of the cell remains uniform throughout. (a I) A snapshot of the cell at t = 1.65 shows an asymmetric formation of high-concentration spots at the cell ends. Streamlines indicate asymmetric velocity fields, with higher velocity and concentration observed on the low-friction or low-viscosity side, respectively. (a II) Outlines of the cell at t = 0, 8, 12, 16, and 20. The cell undergoes asymmetric contraction and progressively migrates towards higher friction or viscosity. At t = 20 the cell is fully polarized with a single spot of high concentration at the rear of the cell. (b) Velocity at the cell ends over time for varying friction and viscosity gradients. In both conditions, contraction and velocity are greater on the low-friction or low-viscosity side. For uniform environments (gray curves), cell end velocities remain nearly symmetric. (c) Cell velocity over time for ten different initial conditions each for three cases: uniform friction/viscosity (gray) and a gradient to the left (green) or right (blue) of the channel. In uniform environments, symmetry is broken spontaneously, and the cell moves either left or right with equal probability. Under gradient conditions, cells consistently migrate toward regions of higher friction or higher viscosity. Parameters: k_right = 1000, k_left = 0, Pe = 22, and η_right = 1, η_left = 0.25, Pe = 15.

Migration up the friction gradient is consistently observed for many different initializations of the initial conditions, Fig. 3(c), illustrating the robustness of this mechanism. Hence, our simulations confirm the previous 1D model results,⁵⁹ highlighting that the frictiotaxis mechanism remains robust when including the effects of cell-shape changes and 3D hydrodynamics.

Depending on the wall coating, previous works reported values of the cell–wall friction coefficient in the range 10³–10⁷ Pa s m⁻¹.^3,59 Moreover, ref. 59 experimentally demonstrated frictiotaxis in channels with a two-fold increase in friction along a channel of ∼100 µm. Translating our simulation units into physical units using the parameter estimates in Table 1, we impose a strong friction gradient of 5 × 10⁸ Pa s m⁻¹ along a channel of 120 µm.

Another emerging area of research is viscotaxis – the directed migration of cells in response to spatial variations in environmental viscosity.^55,60 In this study, we investigate how a viscosity gradient in the surrounding medium influences cell behavior. To incorporate this effect, we modify the previously uniform viscosity, described in eqn (9), by introducing a linear gradient in extracellular viscosity that transitions from η_left at the left end to η_right at the right end of the channel. The intracellular viscosity remains constant throughout and we don't include the additional friction term. Note that the left/right gradient in extracellular viscosity naturally extends to the thin liquid film between cell boundary and channel wall.

Simulation results (Fig. 3 (bottom)) show that cells consistently migrate toward the region of higher viscosity. The underlying mechanism closely resembles that of frictiotaxis: lower extracellular viscosities enable faster cortical flows, leading to faster accumulation of myosin. This initial asymmetry strengthens over time, resulting in cell polarization and directed migration up the viscosity gradient as observed consistently across a range of simulations with different initializations of the initial conditions, Fig. 3(c).

Translating our simulation units to physical units based on the parameter estimates in Table 1, the viscosity gradient imposed in our simulations is of the order of 1 Pa s µm⁻¹. As previously discussed, these viscosity values are considerably higher than those found in typical physiological conditions. However, such high values could be relevant for cell migration within viscous tissues. Tissue viscosities are far greater than those of liquid biological media; they are often in the range 10³–10⁶ Pa s,^49–53 which can give rise to viscosity gradients potentially much larger than the one we utilized.

3.3 Rheo- and barotaxis: cells go with the flow

In many physiological and experimental settings, cells experience fluid flows that influence their polarization and migration. Flow-induced navigation, or rheotaxis, occurs when an externally imposed flow field interacts with intrinsic polarization mechanisms. To simulate this, we impose an inflow velocity u_in at the open channel boundaries. Under no-slip wall conditions – mimicking typical microfluidic experiments – cells consistently migrate in the direction of the flow. This behavior closely resembles that of passive droplets, which we model here by setting a uniform surface tension (σ(c) ≡ 1) and thereby decoupling the surface tension from the myosin concentration.

The imposed flow pushes the cell forward while its surface remains stuck to the channel walls, thus generating a retrograde cortical flow in the cell's frame of reference. This retrograde flow, together with enhanced contraction on the upstream side and reduced contraction downstream, drives myosin accumulation at the rear, reinforcing polarization aligned with the flow (Fig. 4(a)). Active migration thus adds to the passive advection, enabling the cell to move faster than a passive droplet (Fig. 4(d)).

Fig. 4 Guidance of cell migration by inflow and pressure gradients. (a)–(d) Rheotaxis: a cell in a channel with no-slip walls is exposed to an imposed inflow from left to right. (a) A snapshot of the cell and surrounding flow field at t = 3.0 shows the development of a high-concentration spot at the rear, accompanied by retrograde cortical flows. Streamlines indicate flow patterns inside and around the cell. (b) Concentration along the cell surface at different times, showing the emergence and strengthening of polarity. (c) Velocity profiles along the channel centerline at multiple times, revealing asymmetric propulsion in the direction of the imposed flow. (d) Cell velocity over time for two inflow magnitudes (u_in = 0.1 and 0.25), comparing active cells (solid) and passive droplets (dashed); active cells consistently move faster. (e) Barotaxis: cell velocity over time for two pressure drops (p_left = 5 and 10; p_right = 0), comparing active cells (solid) and passive droplets (dashed). Active cells migrate toward the low-pressure end and reach higher speeds than passive droplets.

Using the parameter values in Table 1, the used inflow parameters correspond to velocities on the order of 1 µm min⁻¹. The resulting cell speeds are of the same order, consistent with experimental observations.

In confined tissues and microfluidic setups, pressure gradients are another common physical cue that can bias cell motion – a phenomenon known as barotaxis. Unlike rheotaxis, where velocity is imposed directly, here the driving force is a difference in hydraulic pressure across the channel. In vitro experiments have shown that many cell types preferentially move toward the low-pressure branch in bifurcated channels, with cell polarity strongly correlated to directional bias.^61,62

We simulate barotaxis by imposing a pressure gradient: a positive pressure at the left boundary and zero at the right, thereby inducing flow toward the low-pressure end – analogous to rheotaxis. Again, retrograde flow and rearward contraction both promote myosin accumulation at the rear, producing sustained forward propulsion, closely mirroring the results of rheotaxis. This behavior is robust across different initializations of the initial conditions Fig. 4(e). Using the parameter values in Table 1, the imposed pressure gradient corresponds to approximately 0.07 Pa µm⁻¹. This gradient would result in passive flow speeds below 1 µm s⁻¹, which seem representative of the slow flows observed in confined biological environments.

In rheotaxis, the fluid speed is constrained by the prescribed velocity at the channel openings, which prevents migrating cells from propelling the surrounding incompressible fluid. In contrast, during barotaxis, the fluid velocity at the channel openings is not restricted. This allows migrating cells to accelerate the entire fluid in the channel, resulting in a substantially amplified speed difference between active cells and passive droplets (Fig. 4(e)).

3.4 Topotaxis: cells migrate away from confinement

In addition to sensing biochemical signals, cells are highly responsive to mechanical constraints and spatial confinement. Changes in geometry of the extracellular medium – whether due to tissue structure, extracellular matrix composition, or physical barriers – can significantly influence how cells move.^12,55,58 Given the ubiquity of such conditions in vivo, it is important to examine how varying degrees of confinement shape migratory responses.

Rather than considering complex three-dimensional geometries, we focus here on a simple slanted (conical) channel with an opening angle of approximately 1° to demonstrate that even subtle geometric changes can strongly influence cell behavior.

Initially, a cigar-shaped cell is placed in a slightly slanted channel. Symmetry breaking occurs in two stages: in the first stage, the cell adapts its shape to the channel's geometry and contracts more on the side with greater available space, as it can extend further in the perpendicular direction there. This expansion leads to a localized increase in molecular concentration and flow velocity on the less-constricted side (Fig. 5(a I)), until the cell has adapted to the channel geometry at t ≈ 0.78. This dynamics is also observed in simulations with passive droplets (σ = 1, results not shown).

Fig. 5 Guidance of cell migration by environmental topography (topotaxis). A cell is initialized inside a slanted channel with a wider opening on the right. (a I) Snapshots of the cell at t = 0.75. Initially, surface tension causes the cell to conform to the geometry of the slanted channel (see sketch), resulting in increased concentration and velocity on the right side of the cell. (a II) At later time (t = 8), the cell polarization has reversed and persistent contraction on the left drives net migration to the right. Retrograde flow near the channel wall transports myosin toward the rear, reinforcing directional polarization. (b) Cell conformation shows increasing asymmetry over time (t = 0, 4, 8 and 12), indicating progressive polarization and migration. (c) Velocity at the cell ends over time. After an initial conforming-to-the-channel phase (orange) and contraction phase (red), the cell migrates consistently to the right (blue). (d) Concentration profiles along the cell surface at different times. Initially, two peaks emerge, with the right side being more pronounced (t = 0.78, orange). After t = 2.0 the pattern shifts and polarizes with higher concentration on the left side. (e) Cell velocity over time for ten different initializations of the initial conditions shows that active cells consistently move away from confinement, with a velocity much faster than a passive droplet with constant surface tension (dashed line).

Despite the slight asymmetry in concentration resulting from the cell's initial adaption to the channel geometry, surface tension remains largely homogeneous in active cells at this early stage. Consequently, the cell, much like a passive droplet, seeks to minimize its surface area and moves towards the wider region of the channel. During this movement towards the channel opening, the cell cortex, experiences frictional drag from the nearby channel walls, causing it to move slower than the cell body. The difference in velocity constitutes a retrograde flow in the cell's frame of reference transporting myosin to the more-constricted end. At t ≈ 2, the increasing myosin concentration on the constricted end surpasses the concentration on the less-constricted side (Fig. 5(d)).

This behavior is consistent with simulations of both active cells and passive droplets. However, while passive droplets are driven solely by homogeneous surface tension, active cells migrate significantly faster due to the Marangoni forces created by the myosin concentration gradient. This reinforces the asymmetry and accelerates the migration process to significantly larger speeds (Fig. 5(e)).

These findings suggest that directed migration into constricted environments, as reported in some studies,^12,63 requires additional mechanisms beyond those modeled here.

3.5 Hydrodynamic interactions coordinate motion of nearby cells

Cells are typically not isolated in biological organisms but reside in close proximity to neighboring cells, often forming cooperative groups such as cell trains.⁶⁴ It is therefore of interest to investigate how polarization and migration are influenced by the presence of a neighboring cell. To this end, we consider a minimal system consisting of two identical cells placed side by side within a homogeneous narrow channel. The two cells are not connected by adhesive interactions or signaling pathways but are coupled solely through the shared velocity field in the surrounding medium, Fig. 6(a).

Fig. 6 Coordination of cell migration via hydrodynamic interaction. Two cells are initialized side by side within a channel, each with a nearly uniform surface concentration and elongated shape. There is no direct physical contact between the cells – interaction occurs purely through the surrounding flow field. (a I) At an early time (t = 0.9), both cells contract due to surface tension and behave nearly independently. Each cell induces its own flow, which partially cancel each other out between the cells and reduces velocities on the inner sides. (a II) At t = 30, one cell is slightly more polarized and drives the neighboring cell forward. Both cells exhibit rearward localization of concentration and migrate in the same direction. (b) Velocity profiles along the channel centerline over time. Early profiles show near-zero velocity between the cells (at z = 0) due to mutual cancellation of the opposing flows. (c) Velocity at the cell ends over time. Initially, both cells (blue and green) contract inward asymmetrically. The differences between the cell ends (dotted left, dashed right) of each cell are larger than differences between the two cells. Around t = 20, the cells begin to mutually migrate in one direction. (d) Concentration profiles along the surface of both cells (dotted left, solid right). Early contraction causes accumulation at both cell ends. Over time, the initial concentration differences in each cell decrease, and the profile transitions into a single peak at the rear, indicating polarization. (e) Cell velocity over time for four different initial conditions. In all cases, the two cells ultimately migrate in the same direction, illustrating the robustness of hydrodynamically mediated coordination. Parameters: Pe = 20.

During the initial contraction phase, the cortical flows generated by each cell interfere with one another, leading to a partial cancellation of the flow field in the region between them, Fig. 6(b). Despite this, the cells undergo symmetry breaking, with one of the two cells breaking symmetry earlier than the other. As a result, the cortical flows are slightly stronger in the cell where symmetry broke earlier, which hence achieves a stronger polarity earlier, Fig. 6(d). As this small difference amplifies over time, the more strongly polarized cell begins to move, and – through their coupling via the incompressible fluid between them – it drives the other cell to move in the same direction. The resulting pair migration proceeds at a speed comparable to that of a single cell, suggesting that mutual hydrodynamic interactions facilitate directional coherence, Fig. 6(e).

This minimal two-cell setup provides a foundation for exploring more complex collective behaviors. In particular, it raises intriguing questions about the outcome of interactions between already polarized cells approaching one another from opposite directions – an avenue we propose for future study.

4 Conclusion

We have introduced a phase-field model to study the motility of non-adherent, deformable cells guided by a contractile instability of the cell cortex. The model couples surface and bulk hydrodynamics to surface flow of a diffusible species (myosin), which generates an active contractile force. It accounts for surface viscosity and includes a stabilizing auxiliary velocity field, derived through asymptotic analysis, to maintain robust interface dynamics under strong cortical deformations.

Through simulations, we systematically explored mechanisms guiding non-adherent cell migration. In the absence of external cues, confined cells break symmetry spontaneously via contractile activity, forming a single myosin-rich spot and migrating persistently through retrograde flow. Introducing gradients in the environmental properties revealed distinct taxis behaviors (see Table 2): frictiotaxis and viscotaxis arise from gradients in wall friction and extracellular viscosity, respectively, both directing migration toward higher-resistance areas. Rheotaxis, induced by external flow, aligns polarization with the flow direction, causing active cells to move faster than passive droplets. Barotaxis, driven by pressure differences, directs movement toward low-pressure zones. Topotaxis, triggered by gradients in environmental geometry, leads to migration away from constriction, resembling passive droplet behavior. Finally, we examined cell pairs, showing that hydrodynamic coupling induces collective migration of both cells in the same direction. Because of their lack of focal adhesion, both with the environment and with other cells, non-adherent cells are thought to migrate mostly as single cells. Our findings reveal that hydrodynamic interactions could lead to collective amoeboid migration, which was observed in recent experiments.^64,66

Table 2 Modeling modifications and behaviors in different modes of directed cell migration

Migration mode	Schematic representation	Model modifications	Migration direction	Passive droplets	Ref.
Frictiotaxis		Additional friction force term in the momentum eqn (4)	Towards high friction	✗	59
Viscotaxis		Introduction of a gradient in extracellular viscosity	Towards high viscosity	✗	60
Rheotaxis		Inflow at the channel openings	With the flow	✓	3 and 65
Barotaxis		Additional pressure boundary conditions	Towards low pressure	✓	61 and 62
Topotaxis		Alteration of the channel geometry	Away from confinement	✓	12 and 63

---

We show that mechanical cues, such as friction, viscosity, pressure, flow, and confinement, can each independently guide non-adherent cell migration. Moreover, even weak hydrodynamic interactions enable coordinated motion in multi-cell systems. Notably, passive and active responses differ under specific conditions, emphasizing the role of the contractile instability of the cell cortex in shaping stimulus-specific behaviors.

Translating our simulation units into physical units using the parameter estimates in Table 1, the characteristic time scale of all our simulations is 1 min. Hence, the process of cell polarization in our simulations sets in within seconds and is complete within minutes, which is similar to the duration of ∼10 min of noise-induced cell polarization in experiments.⁴⁵ Therefore, the noise level in our simulations is enough to destabilize an unpolarized cell in a time scale similar to experimental observations. Thus, the gradient values used in our simulations inform the gradients that should be used in experiments to obtain similarly robust guidance of cell migration against noise.

Moreover, while we studied cell propulsion in a simple fluid, our results can also be interpreted in the context of a cell moving through tissue. This broader applicability is supported by our choice of large fluid viscosities (which mimic long-term viscous behavior of tissue) and the observed robustness of our results across various cortical-to-fluid viscosity ratios.

Although simulations were restricted to axisymmetric geometries, the model is fully extensible to general 3D settings. This opens avenues for investigating cell behavior in more complex environments, including branched channels, obstacles, and tissue-like matrices. These situations are relevant for the migration of different cell types, such as immune cells surveying tissues for pathogens or cancer cells invading healthy tissues.⁶⁷ Therefore, the guidance mechanisms identified here may help explain how such cells make directional decisions in their complex microenvironments.

Beyond cell migration, our model provides a foundation for studying other biological processes involving active surfaces – for example, cytokinesis, where a contractile ring emerges in the cell cortex and progressively constricts to divide the cell.^68,69 Minimal models¹⁸ suggest that the underlying mechanics are the same as those explored here. Having demonstrated that our approach is capable of describing cell–cell interaction, it also opens the door to studying active surface dynamics in multicellular systems, for example, those central to early animal embryogenesis,^70,71 an area where computational modeling remains remarkably limited.

Taken together, our results underscore how complex migratory behaviors can emerge from minimal ingredients like cortical contractility, hydrodynamic coupling, and environmental gradients within a unified physical framework. This sets the stage for further theoretical and experimental studies into the principles governing active surface-driven motility.

Author contributions

HLG: formal analysis; software; visualization; writing – original draft; writing – review & editing. CW: software; writing – original draft. PAEH: writing – review & editing. RA: conceptualization; supervision; writing – review & editing. SA: conceptualization; funding acquisition; methodology; supervision; writing – review & editing.

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

The numerical simulation code is available at Zenodo (https://doi.org/10.5281/zenodo.17167539). Installation requires the finite element library AMDiS,⁴¹ which is available online (https://gitlab.math.tu-dresden.de/iwr/amdis).

Appendices

A Supplementary model details

A.1 Extending a surface equation to Ω. In this subsection, we formally justify the diffuse-domain formulation of the surface concentration eqn (8). Therefore, we extend the sharp-interface eqn (3) to the full domain Ω. To this end, we introduce the weak formulation by testing with a suitable test function ψ with compact support in Ω × [0,T]. This condition implies that ψ vanishes on the boundary of Ω as well as at the initial time t = 0 and final time t = T. The weak formulation of eqn (3) can then be written as

Using the Reynolds transport theorem on surfaces yields

Now the domain of integration can be extended to Ω ⊃ Γ by use of the interface distribution δ_Γ.

Using the ordinary Reynolds transport theorem (not on a surface) yields

Going back to the strong formulation gives

0 = ∂_t(δ_Γc) + ∇·(δ_Γcu) − ∇·(δ_ΓP_Γ∇c).

The asymptotic analysis of this equation without the surface projection P_Γ shows that the concentration is constant in normal direction, n·∇c = 0 for ε → 0. In this case P_Γ·∇c = ∇c, hence P_Γ can be omitted. Approximating the surface delta function by |∇ϕ| yields the diffuse interface concentration eqn (8).

A.2 Asymptotic analysis. We will argue in this section, that the solution v(x,ε) of

|∇ϕ|v − |∇ϕ|u − ∇·[|∇ϕ|(n ⊗ n)∇v] = 0 on Ω, (14)

is a constant normal extension of u, i.e. for ε → 0, the equations converge formally to

v = u on Γ, (15)

(n·∇)v = 0 on Γ. (16)

To do so, we utilize matched asymptotic analysis and follow the argumentation of ref. 39 and 72. We start by introducing a new inner coordinate system in a neighborhood of the surface Γ, such that for any x in the neighborhood there is an unique s(x) ∈ Γ with minimal distance to x. Then x can be represented as

x = s(x) + r(x)n = s(x) + εz(x)n,

where r is a signed distance function with r < 0 in Ω₁ and r > 0 in Ω₀. The variable z is a scaled distance function defined by . The phase-field function can be expressed in terms of these coordinates as

We expand u(x,ε) and v(x,ε) outside the interface region in terms of the original coordinate system . This procedure is referred to as the outer expansion. In the neighborhood of Γ we introduce for û(s,z) and v̂(s,z) the inner expansion in terms of the new coordinate system. In an overlapping region, both are valid representations of the same function and thus the following matching conditions hold

Inserting the outer expansions into eqn (14) then yields

0 = |∇ϕ⁽⁰⁾|v⁽⁰⁾ − |∇ϕ⁽⁰⁾|u⁽⁰⁾ − ∇·[|∇ϕ⁽⁰⁾|(n ⊗ n)∇v⁽⁰⁾],

which gives the trivial identity 0 = 0 away from the interface.

In the inner coordinate system, eqn (14) turns to

where we used n·∇_Γ = 0 and ∂_zn = 0. By inserting the inner expansion and comparing powers of ε, we obtain the following condition at order 1/ε²

0 = ∂_z(|∇ϕ|∂_zv̂⁽⁰⁾).

Thus, |∇ϕ|∂_zv̂⁽⁰⁾ is constant in z. This constant value must be zero since |∇ϕ| → 0 for z → ±∞ (see eqn (A.1)). Hence, ∂_zv̂⁽⁰⁾ = 0. Using this result, we obtain at order 1/ε

0 = ∂_z(|∇ϕ|∂_zv̂⁽¹⁾).

Similarly to above, we deduce ∂_zv̂⁽¹⁾ = 0. Using matching condition (A.4) we obtain the desired condition (16), as

(n·∇)v⁽⁰⁾ = 0 on Γ.

Finally, we have at order 1

0 = |∇ϕ|û⁽⁰⁾ − |∇ϕ|v̂⁽⁰⁾ + ∂_z(|∇ϕ|∂_zv̂⁽²⁾) + ∇_Γ·(|∇ϕ|n∂_zv̂⁽¹⁾) = |∇ϕ|û⁽⁰⁾ − |∇ϕ|v̂⁽⁰⁾ + ∂_z(|∇ϕ|∂_zv̂⁽²⁾).

Integrating from z = −∞ to z = +∞ yields

where we used that v̂⁽⁰⁾ is constant in z as shown above and û⁽⁰⁾ is constant in z since u is a continuous function (in fact: û⁽⁰⁾ = u_|Γ). Further, from eqn (A.1) we conclude for z = ±∞ that |∇ϕ| = 0, hence [|∇ϕ|∂_zv̂⁽²⁾]^+∞_−∞ = 0. Dividing by the (non-zero) integrals, we obtain v̂⁰ = û⁰. With conditions (A.2) and (A.3) we conclude v⁰ = u⁰, from which we recover the desired eqn (15).

A.3 Time discretization. After establishing the equation system, it remains to solve the six strongly coupled equations. To reduce the size of the linear equation system that we have to solve in each time step, we decouple them and solve the Stokes–Cahn–Hilliard equations independently from the velocity extension equation and the concentration equation.

We employ the stable linear semi-implicit time discretization from ref. 35 for solving the coupled Stokes–Cahn–Hilliard–Navier system in the n-th time step

where we calculated ϕⁿ⁻¹, vⁿ⁻¹ and cⁿ⁻¹ in the previous time step.

Using the current time steps ϕⁿ and uⁿ, we obtain the projected velocity v from

|∇ϕⁿ|vⁿ − ∇·(|∇ϕⁿ|(ñⁿ ⊗ ñⁿ)∇vⁿ) = |∇ϕⁿ|uⁿ.

Finally, we solve the concentration equation using the previously computed values for ϕⁿ and vⁿ

Note that, even without the term on the right hand side (i.e. in case α = 0), the proposed time discretization ensures mass conservation on the discrete level. This property becomes obvious in the weak from of the equation with a test function φ ∈ H¹(Ω):

This holds especially for φ ≡ 1 and thus

which is the phase-field equivalent of exact mass conservation

However, due to adaptive grid refinement and coarsening, small errors in surface mass may accumulate over time, requiring the mass correction (α > 0) for long simulation times.

A.4 Space discretization. We solve the system of equations in each time step with a finite element method based on the finite element toolboxes DUNE⁴⁰ and AMDiS.^41–43 To decrease the size of the system, we avoid solving the full 3D-problem and assume that the cell shape and myosin concentration distribution are axisymmetric. This holds in particular for the biologically most relevant patterns, which are rings and single spots of increased concentration. The assumption of axisymmetry reduces computations effectively to a 2D domain from which the full 3D-solution can be recovered by rotating the calculated 2D solution 360 degrees. Detailed explanations can be found in the Appendix A.5.

An adaptive grid is employed to accurately resolve the phase field and surface forces, see Fig. 1 (bottom right). Interfacial grid refinement is heuristically chosen, based on the value of the phase field, such that the grid size is h_int where 0.05 < ϕ < 0.95 and h_bulk otherwise.

The numerical approximations ñ and P̃_Γ become less accurate when evaluated farther from the interface. To avoid numerical errors accumulating in the outer areas, we replace |∇ϕ| by

Interchanging |∇ϕ|* with max(|∇ϕ|*, 10⁻¹) in the diffusion terms in eqn (8) and (14), we ensure that the induced linear systems remain regular. For the adjusted problem and the introduced discretization, we use Lagrange-P2 elements for the phase field ϕ, the chemical potential μ, the velocities u and v and the concentration c, only for the pressure p we use a P1-ansatz space.

A.5 Axisymmetric formulation. As hinted before we use the rotationally symmetric setup to our advantage and save computation time by solving a 2D-problem and then rotate the solution around the z-axis turning it into 3D. For this we consider cylinder coordinates (ρ,φ,z) throughout this section. A vector u ∈ [Doublestruck R]³ is then represented by u = u_ρe_ρ + u_φe_φ + u_ze_z and we write u = (u_ρ,u_φ,u_z)^T and use that notation for all vectors and matrices following.

To reduce the complexity of the problem, we utilize that the solution will be axisymmetric, and thus constant in φ-direction. We therefore replace the usual 3D-Cartesian differential operators by their cylindrical equivalent without the φ-direction.

We denote the unit vectors in cylindrical coordinates by e_ρ, e_φ and e_z. It holds

∂_ρe_ρ = 0, ∂_φe_ρ = e_φ, ∂_ze_ρ = 0,

∂_ρe_φ = 0, ∂_φe_φ = −e_ρ, ∂_ze_φ = 0,

∂_ρe_z = 0, ∂_φe_z = 0, ∂_ze_z = 0.

In cylindrical coordinates, the gradient is given by

Specifically, the gradient of a scalar field f is

and the gradient of a vector field u = u_ρe_ρ + u_φe_φ + u_ze_z

The divergence is defined in a similar fashion as the inner product

thus the divergence of a vector field u is given by

We observe for a tensor A ∈ [Doublestruck R]^n×n

Furthermore, the normal vector on the surface Γ is of the form n = n_ρe_ρ + n_ze_z and

Now we utilize that our system is constant in azimuthal direction, in particular all first derivatives in φ-direction vanish and u_φ = 0. With that, the incompressibility condition (5) turns to

For the conservation of momentum, eqn (10),

Note that for shorter notation we write here δ_Γ instead of the |∇ϕ|. Now formulating this in the described axisymmetric setting, yields for the first order derivatives

and for the second order derivatives

Finally, for the concentration eqn (8) the application of rotational operators ∇_R and ∇_R· is straightforward. Similarly, for the Cahn–Hilliard equation, Δϕ is replaced by ∇_R·∇_Rϕ, which is necessary for the chemical potential to approximate the correct 3-dimensional curvature. However, we refrain from using axisymmetric operators for Δμ to better preserve cell volume. Imposing rotational symmetry in this term would effectively increase the extracellular-to-intracellular volume ratio, thereby promoting diffusion of the smaller intracellular phase into the larger extracellular domain – a known artifact in phase-field models.⁷³ Since our objective is to conserve intracellular volume, such rotational symmetry is intentionally avoided in this specific term, thereby lowering droplet loss of mass by decreasing the volume contrast of intracellular to extracellular domains, due to two-dimensional instead of three-dimensional domains. Note that since the primary purpose of the phase field is to track the cellular interface, rotationally symmetric operators are not necessary for Δμ; the equations still approximate the same sharp interface limit.

A.6 Influence of the interface thickness ε. To evaluate the effect of the interface thickness parameter ε, we performed two validation studies. In the first study, we compared the initial cell adaptation to the channel for three different values of ε. The simulation results (Fig. 7 (left)) indicate that the thickness of the fluid film separating the cell from the wall increases with the interface thickness parameter ε, consistent with the expected scaling behavior of phase-field models.

Fig. 7 Influence of interface thickness ε on cell adaptation and cell speed. Left: Cell contour at t = 10.0 for three interface thickness values (ε = 0.04, 0.02, 0.01), showing that the thickness of the thin fluid film between cell and wall is proportional to ε. Right: Cell speed u_cell over time for different interface thicknesses ε. Smaller ε values lead to higher steady-state speeds, while the results converge for sufficiently small ε.

The second validation study investigates the influence of ε on the velocity of a cell that has already been initialized with a prescribed shallow 1-mode deformation, thereby avoiding additional effects from pattern formation. The comparison across different values of ε shows that the cell velocity converges as ε decreases, approaching a constant nonzero value. As expected, cells move slightly faster for smaller ε, since a thinner fluid layer between the cell and the wall results in increased shear interaction and, consequently, a higher net driving force.

Overall, these studies confirm that, for sufficiently small ε, the simulation results are largely independent of the precise choice of the interface thickness, corroborating the robustness of the findings presented in the main text.

A.7 Critical Péclet number and pattern formation. To assess the onset of spontaneous polarization in our model, we performed a systematic set of simulations for different Péclet numbers Pe. Each simulation was initialized with a slight 1-mode perturbation (amplitude 0.002) and we monitor both the spatial concentration profile along the axial coordinate and the concentration difference between the two cell ends over time.

Fig. 8 (left) shows that for Péclet numbers below 15, the concentration profiles remain nearly left-right symmetric. Above Pe ≈ 20, polar patterns consistently emerge, characterized by high concentration at one cell end and low concentration at the other. The slight visible 2-mode pattern results from the presence of the channel walls and the corresponding initial adaption of cell shape. This 2-mode, however, does not further increase over time, indicating that the critical Pe number is larger than 25 for this mode. Higher order modes are not observed to grow in all considered cases.

Fig. 8 Emergence of polar 1-mode concentration patterns depends on the Péclet number. Left: Centered concentration profile along the axial coordinate z for different Péclet numbers at t = 10. For Pe < 20, profiles remain nearly left-right symmetric, while for Pe ≥ 20, distinct polar patterns emerge. The slight visible 2-mode pattern results from the presence of the channel walls and the corresponding initial adaption of cell shape. Right: Time evolution of the left-right concentration difference, c_R − c_L, for varying Péclet numbers. Above the critical value (Pe ≥ 20), a persistent symmetry breaking and increasing polarity between cell ends is observed, while lower Péclet numbers result in negligible polarization.

Fig. 8 (right) quantifies the evolution of the concentration difference between the right and left ends of the cell, c_R − c_L. For Pe < 20, 1-mode polarity remains negligible, while for Pe ≥ 20, a pronounced and persistent increase is observed, signifying spontaneous emergence of polar 1-mode patterns. The establishment of a slight concentration difference for Pe = 20 indicates that this value is very close to the critical Péclet number of the 1-mode.

Acknowledgements

SA acknowledges support from DFG (grants 511509575, 328170591, 417223351). The authors acknowledge computing time on the compute cluster of the Faculty of Mathematics and Computer Science of Technische Universität Bergakademie Freiberg, operated by the computing center (URZ) and funded by the Deutsche Forschungsgemeinschaft (DFG) under DFG grant number 397252409. This research was supported in part by grant NSF PHY-2309135 and the Gordon and Betty Moore Foundation Grant No. 2919.02 to the Kavli Institute for Theoretical Physics (KITP).

Notes and references

1. D. L. Bodor, W. Pönisch, R. G. Endres and E. K. Paluch, Dev. Cell, 2020, 52, 550–562.
2. E. K. Paluch, I. M. Aspalter and M. Sixt, Annu. Rev. Cell Dev. Biol., 2016, 32, 469–490.
3. M. Bergert, A. Erzberger, R. A. Desai, I. M. Aspalter, A. C. Oates, G. Charras, G. Salbreux and E. K. Paluch, Nat. Cell Biol., 2015, 17, 524–529.
4. K. M. Yamada and M. Sixt, Nat. Rev. Mol. Cell Biol., 2019, 20, 738–752.
5. M. Tozluoğlu, A. L. Tournier, R. P. Jenkins, S. Hooper, P. A. Bates and E. Sahai, Nat. Cell Biol., 2013, 15, 751–762.
6. Y.-J. Liu, M. Le Berre, F. Lautenschlaeger, P. Maiuri, A. Callan-Jones, M. Heuzé, T. Takaki, R. Voituriez and M. Piel, Cell, 2015, 160, 659–672.
7. V. Ruprecht, S. Wieser, A. Callan-Jones, M. Smutny, H. Morita, K. Sako, V. Barone, M. Ritsch-Marte, M. Sixt, R. Voituriez and C.-P. Heisenberg, Cell, 2015, 160, 673–685.
8. C. Kang, P. Chen, X. Yi, D. Li, Y. Hu, Y. Yang, H. Cai, B. Li and C. Wu, eLife, 2024, 13, RP96821.
9. K. M. Stroka, H. Jiang, S.-H. Chen, Z. Tong, D. Wirtz, S. X. Sun and K. Konstantopoulos, Cell, 2014, 157, 611–623.
10. T. Singha and P. Sens, Phys. Rev. Lett., 2025, 134, 068401.
11. H. D. Moreau, M. Piel, R. Voituriez and A.-M. Lennon-Duménil, Trends Immunol., 2018, 39, 632–643.
12. A. Reversat, F. Gaertner, J. Merrin, J. Stopp, S. Tasciyan, J. Aguilera, I. de Vries, R. Hauschild, M. Hons, M. Piel, A. Callan-Jones, R. Voituriez and M. Sixt, Nature, 2020, 582, 582–585.
13. P. Recho, T. Putelat and L. Truskinovsky, Phys. Rev. Lett., 2013, 111, 108102.
14. A. C. Callan-Jones and R. Voituriez, New J. Phys., 2013, 15, 025022.
15. A. C. Callan-Jones and R. Voituriez, Curr. Opin. Cell Biol., 2016, 38, 12–17.
16. A. Callan-Jones, Front. Cell Dev. Biol., 2022, 10, 1000071.
17. A. Mietke, F. Jülicher and I. F. Sbalzarini, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 29–34.
18. A. Mietke, V. Jemseena, K. V. Kumar, I. F. Sbalzarini and F. Jülicher, Phys. Rev. Lett., 2019, 123, 188101.
19. L. D. Wittwer and S. Aland, J. Comput. Phys.:X, 2023, 17, 100126.
20. T. Lämmermann and M. Sixt, Curr. Opin. Cell Biol., 2009, 21, 636–644.
21. K. V. Kumar, J. S. Bois, F. Jülicher and S. W. Grill, Phys. Rev. Lett., 2014, 112, 208101.
22. T. Moore, S. K. Wu, M. Michael, A. S. Yap, G. A. Gomez and Z. Neufeld, Biophys. J., 2014, 107, 2652–2661.
23. I. M. Sehring, P. Recho, E. Denker, M. Kourakis, B. Mathiesen, E. Hannezo, B. Dong and D. Jiang, eLife, 2015, 4, e09206.
24. C. A. Weber, C. H. Rycroft and L. Mahadevan, Phys. Rev. Lett., 2018, 120, 248003.
25. A. C. Callan-Jones, V. Ruprecht, S. Wieser, C. P. Heisenberg and R. Voituriez, Phys. Rev. Lett., 2016, 116, 028102.
26. A. Torres-Sánchez, D. Millán and M. Arroyo, J. Fluid Mech., 2019, 872, 218–271.
27. M. Bonati, L. D. Wittwer, S. Aland and E. Fischer-Friedrich, New J. Phys., 2022, 24, 073044.
28. M. Kelkar, P. Bohec and G. Charras, Curr. Opin. Cell Biol., 2020, 66, 69–78.
29. E. Fischer-Friedrich, Y. Toyoda, C. J. Cattin, D. J. Müller, A. A. Hyman and F. Jülicher, Biophys. J., 2016, 111, 589–600.
30. L. E. Scriven, Chem. Eng. Sci., 1960, 12, 98–108.
31. E. Asante-Asamani, M. Dalton, D. Brazill and W. Strychalski, Appl. Math. Mod. Challenges, 2024, 2, 17–37.
32. R. J. Hawkins, R. Poincloux, O. Bénichou, M. Piel, P. Chavrier and R. Voituriez, Biophys. J., 2011, 101, 1041–1045.
33. D. M. Anderson, G. B. McFadden and A. A. Wheeler, Annu. Rev. Fluid Mech., 1998, 30, 139–165.
34. D. Mokbel, H. Abels and S. Aland, J. Comput. Phys., 2018, 372, 823–840.
35. S. Aland, J. Comput. Phys., 2014, 262, 58–71.
36. S. Aland and F. Chen, Int. J. Numer. Methods Fluids, 2016, 81, 657–671.
37. G. Grün and F. Klingbeil, J. Comput. Phys., 2014, 257, 708–725.
38. M. Sussman and E. Fatemi, SIAM J. Sci. Comput., 1999, 20, 1165–1191.
39. A. Rätz and A. Voigt, Commun. Math. Sci., 2006, 4, 575–590.
40. O. Sander, DUNE—The Distributed and Unified Numerics Environment, Springer International Publishing, 2020.
41. S. Praetorius, The Adaptive Multi-Dimensional Simulation Toolbox (AMDiS), a discretization module on top of the Dune framework, 2020, https://gitlab.com/amdis/amdis.
42. S. Vey and A. Voigt, Comput. Visualization Sci., 2006, 10, 57–67.
43. T. Witkowski, S. Ling, S. Praetorius and A. Voigt, Adv. Comput. Math., 2015, 41, 1145–1177.
44. G. Salbreux, G. Charras and E. Paluch, Trends Cell Biol., 2012, 22, 536–545.
45. P. Maiuri, J.-F. Rupprecht, S. Wieser, V. Ruprecht, O. Bénichou, N. Carpi, M. Coppey, S. de Beco, N. Gov, C.-P. Heisenberg, C. Lage Crespo, F. Lautenschlaeger, M. Le Berre, A.-M. Lennon-Dumenil, M. Raab, H.-R. Thiam, M. Piel, M. Sixt and R. Voituriez, Cell, 2015, 161, 374–386.
46. A. Saha, M. Nishikawa, M. Behrndt, C.-P. Heisenberg, F. Jülicher and S. W. Grill, Biophys. J., 2016, 110, 1421–1429.
47. T. Kalwarczyk, N. Ziebacz, A. Bielejewska, E. Zaboklicka, K. Koynov, J. Szymański, A. Wilk, A. Patkowski, J. Gapiński, H.-J. Butt and R. Hołyst, Nano Lett., 2011, 11, 2157–2163.
48. B. R. Daniels, B. C. Masi and D. Wirtz, Biophys. J., 2006, 90, 4712–4719.
49. G. Forgacs, R. A. Foty, Y. Shafrir and M. S. Steinberg, Biophys. J., 1998, 74, 2227–2234.
50. P. Marmottant, A. Mgharbel, J. Käfer, B. Audren, J.-P. Rieu, J.-C. Vial, B. van der Sanden, A. F. M. Marée, F. Graner and H. Delanoë-Ayari, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 17271–17275.
51. K. Guevorkian, M.-J. Colbert, M. Durth, S. Dufour and F. Brochard-Wyart, Phys. Rev. Lett., 2010, 104, 218101.
52. C. Blanch-Mercader, R. Vincent, E. Bazellières, X. Serra-Picamal, X. Trepat and J. Casademunt, Soft Matter, 2017, 13, 1235–1243.
53. C. Pérez-González, R. Alert, C. Blanch-Mercader, M. Gómez-González, T. Kolodziej, E. Bazellieres, J. Casademunt and X. Trepat, Nat. Phys., 2019, 15, 79–88.
54. A. Diz-Muñoz, O. D. Weiner and D. A. Fletcher, Nat. Phys., 2018, 14, 648–652.
55. A. Shellard and R. Mayor, Trends Cell Biol., 2020, 30, 852–868.
56. J. Condeelis, A. Bresnick, M. Demma, S. Dharmawardhane, R. Eddy, A. L. Hall, R. Sauterer and V. Warren, Dev. Genet., 1990, 11, 333–340.
57. C. Shi and P. A. Iglesias, Wiley Interdiscip. Rev.:Syst. Biol. Med., 2013, 5, 631–642.
58. S. SenGupta, C. A. Parent and J. E. Bear, Nat. Rev. Mol. Cell Biol., 2021, 22, 529–547.
59. A. Shellard, K. Weißenbruch, P. A. E. Hampshire, N. R. Stillman, C. L. Dix, R. Thorogate, A. Imbert, G. Charras, R. Alert and R. Mayor, Nat. Commun., 2025, 16, 3811.
60. B. Liebchen, P. Monderkamp, B. ten Hagen and H. Löwen, Phys. Rev. Lett., 2018, 120, 208002.
61. A.-M. Lennon-Duménil and H. D. Moreau, Curr. Opin. Cell Biol., 2021, 72, 131–136.
62. H. D. Moreau, C. Blanch-Mercader, R. Attia, M. Maurin, Z. Alraies, D. Sanséau, O. Malbec, M.-G. Delgado, P. Bousso, J.-F. Joanny, R. Voituriez, M. Piel and A.-M. Lennon-Duménil, Dev. Cell, 2019, 49, 171–188.e5.
63. J. Park, D.-H. Kim and A. Levchenko, Biophys. J., 2018, 114, 1257–1263.
64. I. M. N. Wortel, J. Postat, M. Mihaylova, M. Merino, A. Bhagrath, M. Harris, L. Wouters, L. Wiebke, D. R. Parisi, J. N. Mandl and J. Textor, bioRxiv, 2024, preprint DOI:10.1101/2024.10.21.618803.
65. O. Otto, P. Rosendahl, A. Mietke, S. Golfier, C. Herold, D. Klaue, S. Girardo, S. Pagliara, A. Ekpenyong, A. Jacobi, M. Wobus, N. Töpfner, U. F. Keyser, J. Mansfeld, E. Fischer-Friedrich and J. Guck, Nat. Methods, 2015, 12, 199–202.
66. D.-L. Pagès, E. Dornier, J. de Seze, E. Gontran, A. Maitra, A. Maciejewski, L. Wang, R. Luan, J. Cartry, C. Canet-Jourdan, J. Raingeaud, G. Lemahieu, M. Lebel, M. Ducreux, M. Gelli, J.-Y. Scoazec, M. Coppey, R. Voituriez, M. Piel and F. Jaulin, Sci. Adv., 2022, 8, eabp8416.
67. R. Poincloux, O. Collin, F. Lizárraga, M. Romao, M. Debray, M. Piel and P. Chavrier, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 1943–1948.
68. L. L. Satterwhite and T. D. Pollard, Curr. Opin. Cell Biol., 1992, 4, 43–52.
69. T. D. Pollard, Curr. Opin. Cell Biol., 2010, 22, 50–56.
70. J.-L. Maître, H. Turlier, R. Illukkumbura, B. Eismann, R. Niwayama, F. Nédélec and T. Hiiragi, Nature, 2016, 536, 344–348.
71. J. Firmin, N. Ecker, D. Rivet Danon, Ö. Özgüç, V. Barraud Lange, H. Turlier, C. Patrat and J.-L. Maître, Nature, 2024, 629, 646–651.
72. X. Li, J. Lowengrub, A. Rätz and A. Voigt, Commun. Math. Sci., 2009, 7, 81–107.
73. P. Yue, C. Zhou and J. J. Feng, J. Comput. Phys., 2007, 223, 1–9.

---