Crawling in a fluid

TL;DR

This paper develops a coupled actin-myosin and fluid flow model to understand how mammalian cells initiate motility in fluids, revealing bifurcations and spontaneous oscillations in cell movement.

Contribution

It introduces a novel analytical model linking actin-myosin kinetics with fluid dynamics to explain cell motility onset in suspension.

Findings

Identification of super- and subcritical bifurcations leading to motility.

Discovery of spontaneous polarity oscillations via Hopf bifurcation.

Predicted cell shape trends beyond spherical assumptions.

Abstract

There is increasing evidence that mammalian cells not only crawl on substrates but can also swim in fluids. To elucidate the mechanisms of the onset of motility of cells in suspension, a model which couples actin and myosin kinetics to fluid flow is proposed and solved for a spherical shape. The swimming speed is extracted in terms of key parameters. We analytically find super- and subcritical bifurcations from a non-motile to a motile state and also spontaneous polarity oscillations that arise from a Hopf bifurcation. Relaxing the spherical assumption, the obtained shapes show appealing trends.

Tables

Table 1. Table 1: Integrals of Green’s kernels for a spherical cell.

$j$	1	2	3
$g_{j,l}$	$\frac{l}{\left(2l+1\right)\left(2l+3\right)}$	$\frac{l+1}{\left(2l-1\right)\left(2l+1\right)}$	$\frac{1}{2l+1}$
$k_{j,l}$	$\frac{3}{2\left(2l+1\right)\left(2l+3\right)}$	$-\frac{3}{2\left(2l-1\right)\left(2l+1\right)}$	$-\frac{3}{2\left(2l+1\right)}$

Table 2. Table 2: Estimates of material coefficients and non dimensional parameters definitions.

name	symbol	typical value
cortical thickness	$h$	${10}^{-7}$ m Clark et al. (2013); Turlier et al. (2014)
cortical viscosity	${\eta }_s$	$h\times \left({10}^3-{10}^6\right)$ Pa m s Turlier et al. (2014); Bergert et al. (2015b)
myosin contractility	$\chi c_0^{\mu }$	$h\times \left({10}^2-{10}^3\right)$ Pa m Bergert et al. (2015b); Turlier et al. (2014)
F-actin compressibility	$\alpha c_0^a$	$h\times {10}^3$ Pa.m Hawkins et al. (2011)
myosin diffusion coefficient	$D^{\mu }$	${10}^{-13}-{10}^{-12}$ ${\text{m}}^2{\text{s}}^{-1}$ Uehara et al. (2010); Hawkins et al. (2011)
cell size	$R$	${10}^{-5}$ m
F-actin turnover	$\beta$	${10}^{-2}-{10}^{-1}$  ${\text{s}}^{-1}$ Hawkins et al. (2011); Fritzsche et al. (2013); Turlier et al. (2014)
characteristic length	$l_0=R$	${10}^{-5}$ m
characteristic time	$t_0=R^2/D^{\mu }$	${10}^2-{10}^3$ s
characteristic surface stress	${\sigma }_0=D^{\mu }{\eta }_s/R^2$	$1-{10}^4$ Pa m
contractility parameter	$\overline{\chi }=\chi c_0^{\mu }/{\sigma }_0$	${10}^{-2}-{10}^3$
compressibility parameter	$\overline{\alpha }=\alpha c_0^a/{\sigma }_0$	${10}^{-1}-{10}^3$
turnover parameter	$\overline{\beta }=\beta t_0$	$1-{10}^2$

Equations

$${\sigma }^s={\eta }_s{ϵ}^s+{\eta }_b{\mathsf{I}}^s({\nabla }^s\cdot {\mathbf{u}}^c)+{\mathsf{I}}^s(\chi c^{\mu }-\alpha c^a),$$

$$\frac{{\eta }_{in}+{\eta }_{out}}{2}{u}_i^c(\mathbf{x})=\oint G_{ij}(\mathbf{x},{\mathbf{x}}^{\prime })f_j({\mathbf{x}}^{\prime })dA({\mathbf{x}}^{\prime })$$

$$+$$

$${\dot{c}}^a+{\nabla }^s\cdot \left[c^a({\mathbf{u}}^c-{\mathbf{v}}_s)\right]=\beta (c_0^a-c^a)$$

$${\dot{c}}^{\mu }+{\nabla }^s\cdot \left[c^{\mu }({\mathbf{u}}^c-{\mathbf{v}}_s)\right]=D^{\mu }{\Delta }^s{c}^{\mu },$$

$$v_s=\frac{2}{3}R\frac{\alpha c_1^a-\chi c_1^{\mu }}{2({\eta }_s+{\eta }_b)+R(3{\eta }_{in}+2{\eta }_{out})},$$

$$c^{a,\mu }(\theta )=1+\sum _{l=1}^{\infty }{c}_l^{a,\mu }{P}_l(\cos \theta ),$$

$${\dot{c}}_l^a=J_l^{aa}{c}_l^a+J_l^{a\mu }{c}_l^{\mu },{\dot{c}}_l^{\mu }=J_l^{\mu a}{c}_l^a+J_l^{\mu \mu }{c}_l^{\mu },$$

$${\bar{\chi }}_1=2+\frac{2\bar{\alpha }}{\bar{\beta }},{\bar{\chi }}_2=2+\bar{\alpha }+\bar{\beta }.$$

$${\dot{C}}_1=\frac{{\bar{\beta }}^2}{{\bar{\beta }}^2+\bar{\alpha }\bar{\beta }-2\bar{\alpha }}\left[(\bar{\chi }-{\bar{\chi }}_1)C_1+\nu C_1^3\right],$$

$$\nu =-\frac{6{\bar{\alpha }}^2\bar{\beta }-12{\bar{\alpha }}^2+\bar{\alpha }{\bar{\beta }}^3+6\bar{\alpha }{\bar{\beta }}^2-24\bar{\alpha }\bar{\beta }+{\bar{\beta }}^4}{\bar{\beta }(\bar{\alpha }+2\bar{\beta })({\bar{\beta }}^2+\bar{\alpha }\bar{\beta }-2\bar{\alpha }{)}^2}.$$

$$v_s=\frac{\bar{\beta }}{3({\bar{\beta }}^2+\bar{\alpha }\bar{\beta }-2\bar{\alpha })}\sqrt{\frac{\bar{\chi }-{\bar{\chi }}_1}{-\nu }}.$$

$$G_{ij}(\mathbf{x},{\mathbf{x}}^{\prime })$$

$$K_{ijk}(\mathbf{x},{\mathbf{x}}^{\prime })$$

$$Y_{l,m}(\mathbf{x})=\sqrt{\frac{2l+1}{4\pi }\frac{(l-m)!}{(l+m)!}}{P}_l^m\left(\frac{x_3}{x}\right){\left(\frac{x_1+i{x}_2}{|x_1+i{x}_2|}\right)}^m,$$

$$c^a(\mathbf{x})=\sum _{l=0}^{\infty }{c}_l^a(\mathbf{x})$$

$$c_l^a(\mathbf{x})=\sum _{m=-l}^l{c}_{l,m}^a{Y}_{l,m}(\mathbf{x})$$

$$c^{\mu }(\mathbf{x})=\sum _{l=0}^{\infty }{c}_l^{\mu }(\mathbf{x})$$

$$c_l^{\mu }(\mathbf{x})=\sum _{m=-l}^l{c}_{l,m}^{\mu }{Y}_{l,m}(\mathbf{x})$$

$$f^n(\mathbf{x})=\sum _{l=2}^{\infty }{f}_l^n(\mathbf{x})$$

$$f_l^n(\mathbf{x})=\sum _{m=-l}^l{f}_{l,m}^n{Y}_{l,m}(\mathbf{x})$$

$$U(\mathbf{x})=\sum _{l=1}^{\infty }{U}_l(\mathbf{x})$$

$$U_l(\mathbf{x})=\sum _{m=-l}^l{U}_{l,m}{Y}_{l,m}(\mathbf{x})$$

$${\mathbf{Y}}_{1,l,m}(\mathbf{x})=[{\nabla }^s-(l+1)\mathbf{x}]Y_{l,m}(\mathbf{x})$$

$${\mathbf{Y}}_{2,l,m}(\mathbf{x})=[{\nabla }^s+l\mathbf{x}]Y_{l,m}(\mathbf{x})$$

$${\mathbf{Y}}_{3,l,m}(\mathbf{x})=\mathbf{x}\times {\nabla }^s{Y}_{l,m}(\mathbf{x})$$

$${\mathbf{u}}^c=\sum _{j=1}^3\sum _{l=0}^{\infty }\sum _{m=-l}^l{u}_{j,l,m}^c{\mathbf{Y}}_{j,l,m}(\mathbf{x})$$

$$\mathbf{f}=\sum _{j=1}^3\sum _{l=0}^{\infty }\sum _{m=-l}^l{f}_{j,l,m}{\mathbf{Y}}_{j,l,m}(\mathbf{x})$$

$$f_{1,l,m}^v=-\frac{(l+2)[2{\eta }_s(l^2+l+1)+{\eta }_b(l+1)(l+2)]}{2l+1}\frac{u_{1,l,m}^c}{R^2}-\frac{l(l-1)(l+2)(2{\eta }_s+{\eta }_b)}{2l+1}\frac{u_{2,l,m}^c}{R^2},$$

$$f_{2,l,m}^v=-\frac{(l+1)(l+2)(l-1)(2{\eta }_s+{\eta }_b)}{2l+1}\frac{u_{1,l,m}^c}{R^2}-\frac{(l-1)[2{\eta }_s(l^2+l+1)+{\eta }_b(l-1)l]}{2l+1}\frac{u_{2,l,m}^c}{R^2},$$

$$f_{3,l,m}^v=-{\eta }_s(l+2)(l-1)\frac{u_{3,l,m}^c}{R^2},$$

$$f_{1,l,m}^e=\frac{l+2}{2l+1}\frac{{\left[\chi c^{\mu }-\alpha c^a\right]}_{l,m}}{R}+\frac{f_{l,m}^n}{2l+1},$$

$$f_{2,l,m}^e=\frac{l-1}{2l+1}\frac{{\left[\chi c^{\mu }-\alpha c^a\right]}_{l,m}}{R}-\frac{f_{l,m}^n}{2l+1},$$

$$f_{3,l,m}^e=0.$$

$$\frac{{\eta }_{in}+{\eta }_{out}}{2}{u}_{j,l,m}^c=R{g}_{j,l}{f}_{j,l,m}+({\eta }_{out}-{\eta }_{in})k_{j,l}{u}_{j,l,m}^c,$$

$$u_{1,l,m}^c=\frac{U_{l,m}}{R}\frac{l}{2l+1},u_{2,l,m}^c=\frac{U_{l,m}}{R}\frac{l+1}{2l+1},u_{3,l,m}^c=0\text{ for }l>1.$$

$$U_{l,m}=\frac{R(\chi c^{\mu }-\alpha c^a{)}_{l,m}}{{\lambda }_l}.$$

Full text

Introduction

Cell motility plays an essential role during embryogenesis, tissue renewal and repair, response of the immune system to infection as well as pathological processes such as cancer cell migration. It is also important from an engineering point of view to conceive biomimetic robots.

A longstanding dogma in biology is that mammalian cells need mechanical interaction with a solid substrate in order to move forward. The commonly adopted idea about the source of motion is that actin polymerization and actin contraction assisted by molecular motors (myosin) create a tangential flow of the cell cortex, the retrograde flow. Thrust is then achieved by the transmission of momentum to the outer environment by a dynamic adhesive coupling through transmembrane proteins such as integrins Mogilner (2009). Several modes of such crawling motion are distinguished in the literature, most of which involve ample shape changes Petrie and Yamada (2015); Zhu and Mogilner (2016); Aroush et al. (2017); Paluch and Raz (2013) but crawling in a shape-preserving manner is also possible, as notably observed for fish keratocytes Lee et al. (1994); Rubinstein et al. (2009); Tjhung et al. (2015); Ziebert et al. (2011). A major step forward has been achieved by identifying that leukocytes can migrate in an extracellular matrix without the assistance of integrins Lämmermann et al. (2008), raising the important question about the role of specific adhesion in crawling. Indeed, it is reported that crawling can be mediated by friction only Bergert et al. (2015a) when the cell is confined by a 3D environment. Increasing attention has been recently paid to the understanding of such adhesion-independent crawling Paluch et al. (2016); Renkawitz et al. (2009); Hawkins et al. (2009).

However, it is also reported that some cells can swim in an unconfined fluid Barry and Bretscher (2010); Arroyo et al. (2012); Aoun et al. (2019). Swimming relies on periodic distortions of the cell surface which generate a reaction force from the ambient fluid Lauga and Powers (2009). Beside the special case of ciliated squirmers Stone and Samuel (1996), the mainstream notion of swimming relies on ample shape changes, the ones of flagellas Lauga and Powers (2009) or of the cell body Brennen and Winet (1977); Arroyo et al. (2012); Fackler and Grosse (2008); Farutin et al. (2013); Campbell and Bagchi (2017).

In this letter, we show that swimming may be generically possible also for mammalian cells through the same growth and contraction driven cortical flows operating during crawling and in the absence shape changes. Mechanical interaction with a substrate or extracellular matrix, and ample shape changes are therefore not necessary for crawling. This actuation mechanism necessitates spatial symmetry breaking in order to generate directed motion, in the same way as in the case of crawling on a solid substrate Verkhovsky et al. (1999); Recho et al. (2013); Barnhart et al. (2015). Recently, symmetry breaking in cortex dynamics and the resulting cell shapes have been analyzed for quasi-spherical cells Callan-Jones et al. (2016) but the direct relation between the symmetry breaking and the swimming motion has not been described. By coupling the acto-myosin dynamics to the internal, external, and cortex flow, we provide analytically the functional dependence of the swimming speed on the relevant parameters.

A second finding is the identification of a rich panel of instabilities leading to cell polarization. Depending on three key parameters of the model (actin stiffness, rate of actin turn-over, myosin contractility) we find either a supercritical bifurcation, leading to a continuous transition between static and motile configurations, or a subcritical bifurcation implying a metastable coexistence between the static and motile configurations for a finite range of parameters. For a small enough actin turn-over rate, the system exhibits a Hopf bifurcation resulting in a permanent oscillation of the polarization. Such behaviour is reminiscent of recent experimental evidences of actomyosin oscillations Lavi et al. (2016); Godeau (2016). Finally, we show the shapes of polarized cells if the fixed shape assumption is relaxed. The shapes are obtained in a quasi-spherical limit and show reasonable trends.

Model

We consider a neutrally buoyant cell in a fluid environment. Swimming is achieved by transmission of shear stress by the cortex flow to the surrounding fluid. We leave the question of the exact nature of the transmission mechanism and its efficiency outside the scope of this study, and assume perfect no-slip conditions between the cortex and the fluids inside and outside the cell.

We first focus on a fixed spherical shape and show how spontaneous polarization of acto-myosin dynamics and the resulting cortex flow can drive the deformation-independent swimming of the cell. The cortex of the cell is modeled as a compressible two-dimensional (2D) Newtonian fluid along the surface of the sphere. The surface density of actin meshwork is denoted $c^{a}.$ The distribution of myosin that crosslinks actin is characterized by the concentration $c^{\mu}.$ The velocity field (in the laboratory frame) of the cortex is denoted as $\boldsymbol{u}^{c}$. The surface strain rate tensor reads $\mathsf{\epsilon}^{s}=\boldsymbol{\nabla}^{s}\otimes\boldsymbol{u}^{c}\boldsymbol{\cdot}\mathsf{I}^{s}+\mathsf{I}^{s}\boldsymbol{\cdot}(\boldsymbol{\nabla}^{s}\otimes\boldsymbol{u}^{c})^{T}$ where $\mathsf{I}^{s}$ is the surface projection operator and $\nabla^{s}\equiv\mathsf{I}^{s}\boldsymbol{\cdot}\boldsymbol{\nabla}$ is the surface gradient. In the framework of the active gel theory Kruse et al. (2004), we use the following constitutive equation for the surface stress

[TABLE]

where $\eta_{s}$ and $\eta_{b}$ are 2D shear and bulk viscosities, $\chi$ is the myosin-induced contractility and $\alpha$ is the 2D bulk modulus of the actin meshwork.

The surface forces are calculated as $\boldsymbol{f}=\boldsymbol{\nabla}^{s}\boldsymbol{\cdot}\mathsf{\sigma}^{s}-f^{n}\boldsymbol{n}$, where $\boldsymbol{n}$ is the outward normal and $f^{n}$ is a Lagrange multiplier that enforces the fixed shape Farutin et al. (2012). It satisfies the zero total force condition $\oint f^{n}\boldsymbol{n}dA=0$, where $dA$ is the area element on the boundary of the cell. The forces $\boldsymbol{f}$ are balanced by viscous stresses of the cytoplasm (the fluid occupying the cell interior) and the suspending fluid at the cell boundary. Both are considered incompressible Newtonian fluids of viscosities $\eta_{in}$ and $\eta_{out}$, respectively. With this assumption and at zero Reynolds number, 3D fluid dynamics is governed by Stokes equations, which can be combined into a boundary integral equation Pozrikidis (1992)

[TABLE]

using the boundary conditions at infinity, as well as stress balance and continuity of velocity at the cell surface. Here $\boldsymbol{x}$ and $\boldsymbol{x}^{\prime}$ are points on the cell surface. The Green’s kernels $G_{ij}$ and $K_{ijk}$ are listed in AF et al. (2019).

The swimming speed is defined as the volume-averaged velocity of the cytoplasm, which can be converted ino a surface integral $\boldsymbol{v}_{s}=(\frac{4}{3}\pi R^{3})^{-1}\oint\boldsymbol{x}(\boldsymbol{u}^{c}\boldsymbol{\cdot}\boldsymbol{n})dA$, where $R$ is the cell radius. The fixed shape assumption implies that the cortex velocity must be tangential to the cell surface in the comoving frame: $(\boldsymbol{u}^{c}-\boldsymbol{v}_{s})\boldsymbol{\cdot}\boldsymbol{n}=0$. This condition fixes $f^{n}$.

To close the description we express the conservation equations on the actin and myosin fields $c^{a}$ and $c^{\mu}$:

[TABLE]

where dot denotes time derivative. The term $\beta(c^{a}_{0}-c^{a})$ represents actin turnover, with $\beta$ being the turnover rate and $c^{a}_{0}$ the equilibrium concentration in the cortex. The term $D^{\mu}\Delta^{s}c^{\mu}$ represents the surface diffusion of myosin, $D^{\mu}$ is the diffusion coefficient. The average myosin concentration $c^{\mu}_{0}$ is conserved by eq. (4).

This model contains two different active drivings of cell motility: molecular motors are pulling agents generating positive stresses in (1) while actin turnover in (3) can be associated with pushing agents generating negative stress in (1). The interplay between these two types of agents in cell crawling was analyzed in Carlsson (2011); Recho and Truskinovsky (2015).

Results

System (1)-(4) constitutes a closed set of equations for determining the cortex velocity field and the actomyosin dynamics. The solution strategy is as follows: It can be shown (see AF et al. (2019)) that for a spherical cell the cortex flow is potential, $\boldsymbol{u}^{c}=\boldsymbol{\nabla}^{s}U+\boldsymbol{v}_{s}$, where $U$ is the flow potential. The mechanical part of the problem (1)-(2) is then solved analytically by expansion in scalar or vector spherical harmonics AF et al. (2019). This allows us to express $U$ for a spherical shape as a linear combination of the harmonic coefficients of $c^{a}$ and $c^{\mu}$ and compute the swimming velocity

[TABLE]

where $c_{1}^{a}$ and $c_{1}^{\mu}$ are the first harmonics of the actin and myosin concentration, defined precisely below.

Equations (3) and (4) are nonlinear but still can be tackled analytically by perturbation expansion. We validated the analytical results by numerical solution of eqs. (3) and (4). The details of the numerical method are given in AF et al. (2019).

An important observation is that the dissipation mechanism highlighted in Eq. (26) arises from a combination of cortex, internal and external fluid viscosities. They act as dashpots in parallel, as shown by viscosity additivity (up to numerical prefactors). Interestingly, the speed is finite even if $\eta_{out}=0$: while the suspending fluid is essential for swimming, its viscosity is not decisive for setting the value of the swimming speed. We therefore consider the limit of $\eta_{in}=\eta_{out}=0$ in the rest of this study, motivated by the fact that in physiological conditions most of the dissipation occurs in the cortex. We also set $\eta_{b}=0$ for simplicity. Along with these assumptions, we use from now on as characteristic scales $R$ for space, $R^{2}/D^{\mu}$ for time, $D^{\mu}\eta_{s}/R^{2}$ for surface stresses, $c_{0}^{a}$ and $c^{\mu}_{a}$ for actin and myosin concentrations but keep the same symbols to avoid new notations. Three non-dimensional parameters fully define the problem, the motor contractility $\bar{\chi}\equiv\chi c_{0}^{\mu}R^{2}/(D^{\mu}\eta_{s})$, the compressibility of actin $\bar{\alpha}\equiv\alpha c_{0}^{a}R^{2}/(D^{\mu}\eta_{s})$ and the turnover rate of actin $\bar{\beta}\equiv R^{2}\beta/D^{\mu}$.

A homogeneous solution for actin and myosin ($c^{a}=c^{\mu}=1$) always exits. This solution yields a zero cortical flow field and the cell is at rest. This solution can become unstable leading to a polarized cortex dynamics. Such symmetry breaking results in a cortical flow and, in turn, in spontaneous cell motion.

To quantify such process, we expand scalar fields in axisymmetric spherical harmonics as

[TABLE]

where $P_{l}$ are Legendre polynomials and $\theta$ is the polar angle between the swimming direction and the vector position on the sphere counted from the center of mass. Linearizing eqs. (3) and (4) about the homogeneous solution, yields

[TABLE]

where $J_{l}^{aa}=-l(l+1)\bar{\alpha}/\lambda_{l}-\bar{\beta},$ $J_{l}^{a\mu}=l(l+1)\bar{\chi}/\lambda_{l},$ $J_{l}^{\mu a}=-l(l+1)\bar{\alpha}/\lambda_{l},$ $J_{l}^{\mu\mu}=l(l+1)\bar{\chi}/\lambda_{l}-l(l+1),$ and $\lambda_{l}=2(l^{2}+l-1)$. The growth rate $\omega$ of a perturbation in (6) therefore obeys a quadratic equation.

When the real part of $\omega$ vanishes (onset of instability), the imaginary part can either be zero or non zero. The first case defines a steady bifurcation, while the second one defines a Hopf bifurcation. A simple analysis shows that the minimum contractility at which the instability occurs is controlled by the first harmonic $l=1$. The instability takes place if $\bar{\chi}>\bar{\chi}_{c}=\min(\bar{\chi}_{1},\bar{\chi}_{2}),$ where

[TABLE]

We get a steady bifurcation for $\bar{\chi}_{1}<\bar{\chi}_{2}$ and a Hopf one otherwise. Once the contractility exceeds the critical value given by (7) the perturbations grow exponentially with time and nonlinear terms can no longer be neglected. In general resorting to numerical resolution is necessary. However, in the vicinity of the bifurcation point a perturbation analysis is possible.

First, let us consider the case of a steady bifurcation. A nonzero first harmonic (concentration polarity) induces a spontaneous cell motion (See (26)). This propulsion mode is driven by the appearance of a retrograde flow of actin, directed from a front region to a rear region and compensated by actin creation at the front ($\bar{\beta}(1-c^{a})>0$, polymerization) and degradation at the rear ($\bar{\beta}(1-c^{a})<0$, depolymerization). Figure 1 shows the acto-myosin and the cortex flow patterns. The nature of this instability is similar to the one presented in Recho et al. (2013) in a simplified picture.

The unstable eigenmode is given by $C_{1}=2(\bar{\alpha}+\bar{\beta})c_{1}^{\mu}-2\bar{\alpha}c_{1}^{a}$. To obtain the velocity close to the bifurcation threshold, we must also consider the second harmonic. We set $\dot{c}^{a}_{2}=\dot{c}^{\mu}_{2}=0$, because relaxation times for these modes are much smaller than that of $C_{1}$. This yields the expressions of $c^{a}_{2}$ and $c^{\mu}_{2}$ as a function of $C_{1},$ which upon substitution into the expression of $\dot{C}_{1}$ leads to

[TABLE]

The sign of $\nu$ dictates the nature of the steady bifurcation with $\nu<0$ corresponding to a supercritical bifurcation (nonlinearity saturates linear growth) and $\nu>0$ corresponding to a subcritical one (nonlinearity amplifies linear growth).

In the case of a supercritical bifurcation, (8) has a stable fixed point given by $C_{1}=\sqrt{(\bar{\chi}-\bar{\chi}_{1})/\nu}$. The swimming speed then reads

[TABLE]

The bifurcation diagram in Fig. 2A shows both the analytical results and the full numerical simulations.

When the steady bifurcation is subcritical (i.e. first order transition), velocity abruptly switches to a finite value. Due to the finite jump, the perturbative expansion leading to eq.(8) is illegitimate. One has thus to resort to a numerical solution. Figure 2B shows the results. There is a region of metastable coexistence between a non-polarized static state and a polarized motile one between a turning point located at $\bar{\chi}=\bar{\chi}_{3}$ and the linear stability threshold $\bar{\chi}=\bar{\chi}_{1}$. A finite perturbation at a constant contractility can therefore be sufficient to initiate or arrest motion in this range.

When the bifurcation is of Hopf type, it is always supercritical with periodic oscillations of velocity emerging continuously from the static trivial solution at the critical value of contractility $\bar{\chi}_{2}$. A second critical contractility $\bar{\chi}_{4}$ exists at which the periodic oscillations disappear and the system abruptly jumps to a stable motile steady-state. The steady-state branch can be continued for $\bar{\chi}\geq\bar{\chi}_{3}.$ An unstable branch bent backwards emerges continuously from the uniform solution at critical contractility $\bar{\chi}_{1},$ and annihilates with the stable steady-state branch. The full bifurcation diagram shown in Fig.3 is reminiscent of a heteroclinic bifurcation. In the range $[\bar{\chi}_{3},\bar{\chi}_{2}]$, we again observe the coexistence of static and motile solutions while in the interval $[\bar{\chi}_{3},\bar{\chi}_{4}]$, oscillatory and motile solutions are metastable potentially giving rise to complex cell gaits in the presence of biological noise.

We summarize these results on the bifurcation nature in a phase diagram in Fig. 4 where the critical value of $\bar{\chi}_{c}$ is shown in color code.

Relaxing the fixed shape assumption, we have determined the trend of the cell shape. For that purpose we have assumed the shape of the cell to be preserved by surface tension (see AF et al. (2019)). A typical shape is shown in Fig.2C having a satisfactory tendency (mushroom-like flattened shape) compared to some real cell shapes Ruprecht et al. (2015).

Discussion

A rich panel of patterns and dynamics due to a spontaneous symmetry breaking of the acto-myosin kinetics is revealed. The cell thus swims thanks to a retrograde flow of the cortex grasping on the external fluid. Orders of magnitudes of the three key parameters $\bar{\alpha}$, $\bar{\beta}$ and $\bar{\chi}$ can be estimated from available experimental data $\bar{\alpha}\sim 10^{-1}-10^{3}$, $\bar{\beta}\sim 1-10^{2}$ and $\bar{\chi}\sim 10^{-2}-10^{3}$ (see AF et al. (2019)). The wide ranges are due to disparate values of $\eta_{s}$ in the literature Wu et al. (2018). The rectangle in Fig. 4 shows the ranges of estimated parameters. It is interesting to see that this corresponds to a rich region of dynamics including the various bifurcations encountered in this study. The swimming speed (9) yields (in the parameter range shown by rectangle in Fig. 4) in physical units about $v_{s}\sim 1-10D^{\mu}/R\sim 0.1-1\;\mu\text{m}/\text{min}$, which is a reasonable value for many mammalian cells Jilkine and Edelstein-Keshet (2011).

How the cortex flow is transmitted to the external fluid is still unclear. Three pathways are possible (i) through the membrane of the cell via recirculation of phospholipids, (ii) through transmembrane proteins (iii) through bumps that are advected backward by cortex in the form of waves (or shape changes). For axisymmetric flow (as considered here) the first scenario is impossible (it can be shown that the membrane is at rest in the swimming frame), except if endo/exocytosis is operating as reported recently Jones et al. (2006). The second scenario has been recently successfully applied to the T-lymphocyte swimming Aoun et al. (2019). The third scenario would require a cooperative motion of bumps along the membrane but numerical evaluation seems to rule out this possibility Aoun et al. (2019).

Our results show that shape-invariant actomyosin-based motility, a hallmark of crawling, can occur away from any solid substrate, provided that the retrograde flow of actin can transmit momentum to the surrounding fluid. This is consistent with recent experimental observations Aoun et al. (2019) and a theoretical model of active nematic droplets Giomi and DeSimone (2014). The rich bifurcation structure that allows actomyosin to break symmetry is also reminiscent of the one known for cells crawling on a solid substrate Recho et al. (2013); Tjhung et al. (2015); Lavi et al. (2016) pointing to the interpretation that the interaction of the actomyosin contractility with its own viscous resistance to deformation and turn-over is sufficient to propel the cell robustly on a solid and in a fluid.

Acknowledgments

We are grateful to Karin John for her role in the initial discussions leading to this work. AF and CM thank CNES (Centre National d’Etudes Spatiales) and the French-German University Programme “Living Fluids” (Grant CFDA-Q1-14) for financial support. JE and PR are supported by a CNRS Momentum grant, ANR-11-LABX-0030 “Tec21” and IRS “AnisoTiss” of Idex Univ. Grenoble Alpes. All are members of GDR 3570 MecaBio and GDR 3070 CellTiss of CNRS. The computations were performed using the Cactus cluster of the CIMENT infrastructure, supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA) and the authors thank Philippe Beys who manages the cluster.

Supplementary information

The supplementary information contains the explicit expressions for the Green’s kernels of the Stokes equations, the calculation of the cortex flow in terms of the concentration fields, the description of the numerical procedure, the technique of shape reconstruction in the quasi-spherical limit and the physical values used to estimate non-dimensional parameters.

I Green’s kernels

The following kernels are to be used in eq. (2) of the main text:

[TABLE]

II Full solution

The results in this section are presented in the most general form, without any assumptions about the relative values of $\eta_{in}$, $\eta_{out}$, $\eta_{s}$, and $\eta_{b}$.

II.1 Spherical Harmonics

All fields on the cell surface are expanded in spherical harmonics of the vector pointing from the center of the cell to a given point of its surface. The spherical harmonics are defined as

[TABLE]

where $P_{l}^{m}$ are associated Legendre polynomials. The following expansions are used

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Note that $f^{n}_{1}=0$ is zero, as required by the condition of the total force acting on the cell being equal to zero, $f^{n}_{0}$ is irrelevant because it just shifts the osmotic pressure drop across the membrane, and $U_{0}$ is irrelevant because only gradients of $U$ enter equations.

We define the vector spherical harmonics as

[TABLE]

[TABLE]

[TABLE]

The following expansions are used

[TABLE]

[TABLE]

II.2 Force calculation

The force $\boldsymbol{f}$ can be represented as a sum of two contributions $f_{j,l,m}=f^{v}_{j,l,m}+f^{e}_{j,l,m},$ where

[TABLE]

The amplitudes $f^{v}_{j,l,m}$ contain the contribution of the surface viscosity terms in eq. (1) of the main text, while the amplitudes $f^{e}_{j,l,m}$ contain contributions of myosin contractility, cortex compressibility, and Lagrange multiplier $f^{n}\boldsymbol{n}$.

II.3 Fluid dynamics

The integrals in eq. (2) of the main text can be calculated analytically for a spherical cell. The results are

[TABLE]

where the coefficients $g_{j,l}$ and $k_{j,l}$ are listed in table 1.

II.4 Explicit solution

We note that the Green’s kernels (I) are diagonal in the basis of vector spherical harmonics for spherical shape. Furthermore, we see that the amplitudes of the surface viscosity force $f^{v}_{3,l,m}$ depend only on $u^{c}_{3,l,m}.$ This implies that the vector spherical harmonics with first index 3 are completely decoupled from the two other types. Since $f^{e}_{3,l,m}=0$ for all $l$ and $m$, we conclude that $u^{c}_{3,l,m}=0$ for all $l$ and $m$ as well. Adding the fixed shape condition $(\boldsymbol{u}^{c}-\boldsymbol{v}_{s})\boldsymbol{\cdot}\boldsymbol{n}=0$ yields that $\boldsymbol{u}^{c}-\boldsymbol{v}_{s}$ can be written as a surface gradient of some surface potential $U$, as used in the main text. Or, in spherical harmonics,

[TABLE]

Solving the equations (21), (22), and (23) for $U$ and $f^{n}$ yields

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Using the expressions (24) and (26), the retrograde flow and the swimming velocity can be expressed as a function of the concentration fields. The time evolution equations of the concentration fields are obtained by substituting $U$ into eqs. (3) and (4) of the main text. The following equation can be used to reduce all calculations to scalar spherical harmonics

[TABLE]

III Numerical procedure

The numerical procedure consists in representing $c^{a}$, $c^{\mu}$ and $U$ by the amplitudes of the spherical harmonics for all values of $l<l_{max}$ and $|m|\leq l$, where $l_{max}$ is the cut-off value. We take $l_{max}=64$ for such calculations. We observed that regardless of the initial conditions, the dynamics relaxed to an axisymmetric solution. We therefore also performed simulations with the shape assumed axisymmetric from the beginning, which is achieved by setting all amplitudes for $m\neq 0$ to zero. With this assumption, $l_{max}=1024$ was used, which proved to be necessary for strongly polarized cells. The eqs. (3) and (4) of the main text were solved by an explicit Euler scheme by truncating the harmonic expansion of the advection terms to $l<l_{max}$. The time step was chosen small enough to avoid the instability due to the stiffness of the diffusion equation (typically $10^{-4}$ in non-dimensional units). In some cases, a small diffusion of actin (diffusion coefficient $10^{-3}$ in non-dimensional units) was added to enhance the stability of the actin advection equation. The steady-state branches in Figs. 2 and 3 of the main text were obtained by solving eqs. (3) and (4) of the main text with $\dot{c}^{a}=\dot{c}^{\mu}=0$ using Newton’s method.

IV Model for the cell shape

The calculation of the shape follows the method used in Ref. [28] of the main text. Spherical shape of a cell in suspension can be physically achieved by a combination of high osmotic pressure $\Delta P$ inside the cell and the inextensibility of the membrane. Further in this section, we allow the shape of the cell to deviate from a sphere, albeit weakly, taking the leading terms in the small-deformation expansion. We parametrize the cell shape by a shape function $\rho$

[TABLE]

where $|\boldsymbol{x}|$ is the distance from the center of the cell to a given point on its boundary. $\rho_{0}=0$ to the leading order in deformation because of the conservation of the membrane area. $\rho_{1}=0$ because this term corresponds to a translation of the cell to the leading order. We show below that $\rho$ scales as $\Delta P^{-1},$ which justifies an expansion in powers of $\rho$. We consider the quasi-spherical limit, taking the leading terms in such expansions. The function $\rho(\boldsymbol{x})$ is expanded in spherical harmonics of $\boldsymbol{x}$ to be used below

[TABLE]

Assuming the tension of the membrane $\zeta_{0}$ to be homogeneous (unaffected by the cortex flow) we can write for the tension force $\boldsymbol{f}^{m}$

[TABLE]

where $H$ is the mean curvature of the membrane (sum of the principal curvatures) and $H_{0}=2/R$ is the value of $H$ for a perfectly spherical cell. The term $-\zeta_{0}H_{0}\boldsymbol{n}$ in eq. (31) corresponds to an isotropic compression of the fluid inside the cell, which is balanced by the osmotic pressure. This relates the tension of the membrane $\zeta_{0}$ to the pressure difference by the Laplace law:

[TABLE]

The term $(H-H_{0})\boldsymbol{n}\zeta_{0}$ in eq. (31) corresponds to a position-dependent normal force, which we identify with the Lagrange multiplier $f^{n}\boldsymbol{n}$ used to maintain the shape of the cell. This justifies that $H-H_{0}$ scales as $\zeta_{0}^{-1}$ or, equivalently, as $\Delta P^{-1}$ for fixed $f^{n}$. Since $f^{n}$ is governed by the actomyosin dynamics, as follows from eq. (27), we obtain that the shape of the cell can be indeed made as close to a sphere as necessary by choosing $\Delta P$ large enough. This shows that all calculations made for perfectly spherical cells remain valid to the leading order in the limit of large $\Delta P$ even if the spherical-shape condition is relaxed.

The mean curvature can be related to the shape function by

[TABLE]

This yields the final relation between the shape harmonics $\rho_{l,m}$ and the Lagrange multiplier $f^{n}$:

[TABLE]

V Physical parameters

We list in table. 2 the physical data we have considered to obtain rough estimates of the three non-dimensional parameters entering in the model.