Cell motility, synchronization, and cell traction orientational order

TL;DR

This paper develops a theoretical model for cell motility based on traction forces and synchronization, revealing conditions for spontaneous migration and wave-like synchronization patterns, aligning with experimental observations.

Contribution

It introduces a novel tensor-based characterization of cell traction that incorporates synchronization, and demonstrates how non-isochrony leads to spontaneous cell motility.

Findings

Spontaneous transition to motility requires non-isochrony.

Synchronization patterns involve wave propagation.

Model aligns with experimental actin oscillator data.

Abstract

Suspensions of swimming micro-organisms provide examples of coordinated active dynamics. That has stimulated the study of a phenomenological theory combining synchronization and polar order in active matter. Here, we consider another example inspired by the traction forces of migrating cells. The novelty, in this case, is the global force-free nature of the traction force field. Such a constraint is absent in the case where the vector field describes swimming speeds in micro-organisms suspensions. Cell traction is characterized by means of a complex tensor quantity, that generalizes the nematic orientation tensor to incorporate the ability of particles to synchronize, and cell motility depends on this quantity being non-zero. We provide a realization of migrating cell which comprises an assembly of dipolar elements exerting traction on a fluid substrate. This model indicates that spontaneous transition to the motile state is possible but requires (as in the case of synchronization) non-isochrony of oscillations and involves subtle synchronization patterns, associated to propagation of waves. Such results are consistent with recent experimental work relating motility to the synchronization of actin oscillators at the periphery of migrating cells.

Equations

$$c(\mathbf{x},\hat{\mathbf{u}},\varphi ,t)=$$

$$⟨\frac{1}{N}\sum _N\delta (\mathbf{x}-{\mathbf{x}}_n(t))\delta (\hat{\mathbf{u}}-{\hat{\mathbf{u}}}_n(t))\delta (\varphi -{\varphi }_n(t))⟩$$

$$\rho (\mathbf{x},t)=\int d\hat{\mathbf{u}}\int d\varphi c(\mathbf{x},\hat{\mathbf{u}},\varphi ,t)$$

$$\mathbf{p}(\mathbf{x},t)=\int d\hat{\mathbf{u}}\int d\varphi \hat{\mathbf{u}}c(\mathbf{x},\hat{\mathbf{u}},\varphi ,t)$$

$$\mathbb{S}(\mathbf{x},t)=\int d\hat{\mathbf{u}}\int d\varphi [\hat{\mathbf{u}}\otimes \hat{\mathbf{u}}-\frac{\mathbb{I}}{2}]c(\mathbf{x},\hat{\mathbf{u}},\varphi ,t)$$

$$\Phi (\mathbf{x},t)=\int d\hat{\mathbf{u}}\int d\varphi e^{i\varphi }c(\mathbf{x},\hat{\mathbf{u}},\varphi ,t)$$

$$\mathbf{\Pi }(\mathbf{x},t)=\int d\hat{\mathbf{u}}\int d\varphi e^{i\varphi }\hat{\mathbf{u}}c(\mathbf{x},\hat{\mathbf{u}},\varphi ,t)$$

$$\Sigma (\mathbf{x},t)=\int d\hat{\mathbf{u}}\int d\varphi e^{i\varphi }[\hat{\mathbf{u}}\otimes \hat{\mathbf{u}}-\frac{\mathbb{I}}{2}]c(\mathbf{x},\hat{\mathbf{u}},\varphi ,t).$$

$$\mathbf{\Pi }(\mathbf{x},t)=\frac{1}{N}⟨\sum _n{\hat{\mathbf{u}}}_n{e}^{i{\varphi }_n}\delta (\mathbf{x}-{\mathbf{x}}_n)⟩.$$

$$\Sigma (\mathbf{x},t)=\frac{1}{N}⟨\sum _n(\frac{\mathbb{I}}{2}-{\hat{\mathbf{u}}}_n{\hat{\mathbf{u}}}_n)e^{i{\varphi }_n}\delta (\mathbf{x}-{\mathbf{x}}_n)⟩.$$

$$\Sigma (\mathbf{x},t)=\Phi \mathbb{S}$$

$$\mathcal{F}(\mathbf{x})=f_0\sum _{\alpha }\cos {\varphi }_{\alpha }\sum _n{\hat{\mathbf{u}}}_n^{\alpha }\delta (\mathbf{x}-{\mathbf{x}}_n^{\alpha }).$$

$$⟨{\mathcal{F}}_i(\mathbf{x})⟩\propto f_0{l}_{0}Re[{\nabla }_k{\Sigma }_{ik}(\mathbf{x})]$$

$$\eta {\nabla }^2\mathbf{v}-\nabla p=-\mathcal{F}(\mathbf{x});\nabla \cdot \mathbf{v}=0.$$

$$v_i(\mathbf{x})=\int d\mathbf{y}{H}_{ij}(\mathbf{x}-\mathbf{y}){\mathcal{F}}_j(\mathbf{y})$$

$$0\approx -\gamma ({\dot{\mathbf{x}}}_n^{\alpha }-\mathbf{v}({\mathbf{x}}_n^{\alpha }))+{\mathbf{f}}_n^{\alpha }$$

$${\omega }_a^{\alpha }\approx {ϵ}_{abi}{\hat{u}}_b^{\alpha }\sum _{\beta \ne \alpha }\frac{{\mathcal{O}}_i^{\alpha \beta }}{L^{\alpha }}.$$

$${\dot{L}}^{\alpha }=2F^{\alpha }/\gamma +\sum _{\beta \ne \alpha }{\hat{\mathbf{u}}}^{\alpha }\cdot {\mathcal{I}}^{\alpha \beta }$$

$${\dot{F}}^{\alpha }=-\mathcal{K}{d}_{\alpha }+\mathcal{M}{F}^{\alpha }(1-\mathcal{S}{d}_{\alpha }^2)/\gamma +\mathcal{A}{d}_{\alpha }^3$$

$${\dot{R}}_{\alpha }\sim (\lambda -\beta R_{\alpha }^2)R_{\alpha }+l_0^2{\hat{u}}_i^{\alpha }{\hat{u}}_j^{\beta }({\hat{u}}_k^{\alpha }{\hat{u}}_l^{\beta }+{\hat{u}}_l^{\alpha }{\hat{u}}_k^{\beta })\times$$

$$\frac{R_{\beta }\gamma {\omega }_0}{4}\sin ({\varphi }_{\beta }-{\varphi }_{\alpha })T_{ijkl}(x_{\alpha }-x_{\beta })$$

$${\dot{\varphi }}_{\alpha }\sim -\Delta -\chi R_0^2-\frac{l_0^2{\omega }_0\gamma }{4}{\hat{u}}_i^{\alpha }{\hat{u}}_j^{\beta }({\hat{u}}_k^{\alpha }{\hat{u}}_l^{\beta }+{\hat{u}}_l^{\alpha }{\hat{u}}_k^{\beta })\times$$

$$[\cos ({\varphi }_{\beta }-{\varphi }_{\alpha })+\frac{\chi }{\beta }\sin ({\varphi }_{\beta }-{\varphi }_{\alpha })]T_{ijkl}(x_{\alpha }-x_{\beta })$$

$$\dot{\mathbf{C}}=\frac{1}{M}\sum _n\sum _{\alpha }{\dot{\mathbf{x}}}_n^{\alpha }$$

$$⟨{\dot{C}}_i⟩\propto \frac{1}{24i}\int d{x}_{\alpha }\int d{x}_{\beta }{U}_{ijklu}(x_{\alpha }-x_{\beta }){\mathcal{G}}_{jklu}[\Sigma ,\Phi ].$$

$${\mathcal{G}}_{jklu}[\Sigma ,\Phi ]=\{{\Sigma }_{kl}^{\alpha }({\Sigma }_{uj}^{\beta }{)}^{\ast }-c.c.$$

$$+{\Phi }^{\alpha }({\Sigma }_{uj}^{\beta }{)}^{\ast }{\delta }_{kl}/2-c.c\}+\{u\leftrightarrow l\}$$

$${\Sigma }_{kl}^{\alpha }({\Sigma }_{uj}^{\beta }{)}^{\ast }-c.c.=|{\Phi }^0{|}^2[{\mathbb{S}}_{uj}^{\beta }{\mathbb{S}}_{kl}^{\alpha }-{\mathbb{S}}_{uj}^{\beta }{\mathbb{S}}_{kl}^{\alpha }]=0.$$

$${\Phi }^{\alpha }({\Sigma }_{uj}^{\beta }{)}^{\ast }{\delta }_{kl}-c.c.=|{\Phi }^0{|}^2[{\mathbb{S}}_{uj}^{\beta }{\delta }_{kl}-{\mathbb{S}}_{uj}^{\beta }{\delta }_{kl}]=0$$

$${\mathcal{G}}_{jklu}[\Sigma ,\Phi ]=\{{\mathbb{S}}_{uj}^0{\mathbb{S}}_{kl}^0[{\Phi }^{\ast }(x_{\beta })\Phi (x_{\alpha })-c.c.]$$

$$+{\mathbb{S}}_{uj}^0\frac{{\delta }_{kl}}{2}[\Phi (x^{\alpha })-c.c.]\}+\{l\leftrightarrow u\}$$

$$⟨{\dot{C}}_i⟩\sim \frac{{\mathcal{G}}_{jklu}^{(0)}}{24}\times (\mathbf{k}\cdot \hat{\mathbf{x}})\int d{x}_{\alpha }\int d{x}_{\beta }(x_{\beta }-x_{\alpha })U_{ijklu}$$

$${\partial }_t\Phi \sim -D_{\varphi }\Phi -\chi \mathcal{Z}\{[(\frac{{\delta }_{il}}{2}\rho +{\mathbb{S}}_{il})\int d{x}_{\alpha }{T}_{ijkl}\times$$

$$(\frac{{\delta }_{jk}}{2}{\Phi }^{\alpha }+{\Sigma }_{jk}^{\alpha })]+[k\leftrightarrow l]\}.$$

$$\displaystyle\partial_{t}\Sigma_{fg}\sim-2D\Sigma_{fg}-\chi\frac{\mathcal{Z}}{8}\int dx_{\alpha}\Big{(}[\{\rho(T_{fjkg}+T_{gjkf})$$

Full text

.1 Introduction

Motile cells resemble microscopic active (i.e. self-driven) droplets. Thanks to this analogy, model cells Kruse et al. (2006); Ziebert et al. (2012); Tjhung et al. (2012); Ruprecht et al. (2015); Bergert et al. (2015) which rely on soft active matter theory Marchetti et al. (2013); Prost et al. (2015) have deepened our understanding of cell motility. So far, however, most of these theoretical models were focused on paradigmatic cases, like that of fish-keratocytes Danuser et al. (2013), where cells form lamellipodia. The motion of such cells involves regular flows of actin, motor contractility Kozlov and Mogilner (2007); Rubinstein et al. (2009), and the interplay with adhesion which makes it possible to picture, in theoretical models, adhesion and force-generation as time-independent processes. Unsteady oscillatory motions, however, are observed in other types of migrating cells, like neurons Jiang et al. (2015) or amoeba Ehrengruber et al. (1996); Tanimoto and Sano (2014); Hoeller et al. (2016); Negrete et al. (2016). Analogous behavior is observed also when ameboid cells and neutrophils are suspended in fluids Barry and Bretscher (2010). In adherent cells Burnette et al. (2011) forces and adhesion have to be coordinated to enable cell motion Giannone et al. (2007) and molecular determinants for the control of an excitable, oscillatory, system have been recently identified in fibroblasts Corallino et al. (2018).

One promising strategy to account for the spontaneous transition De Magistris et al. (2014) from non-motile behavior to directed migration is based on the following hypothesis: coordinated unsteady spatiotemporal patterns observed in migrating cells are self-organisation phenomena Parent and Weiner (2013) involving feedbacks Allard and Mogilner (2012) and bio-mechanical interactions Giannone et al. (2007); Allard and Mogilner (2012) among different sub-cellular structures. The lack of evidence of any central cellular unit regulating motility and chemotaxis Giannone et al. (2007); Allard and Mogilner (2012); Parent and Weiner (2013) supports this view. Synchronization among the various cell parts is one possibility that was invoked in some works studying cell motility Hermans et al. (2013); Aubry et al. (2015). More recently, synchronization of sub-cellular actin oscillators was shown to be crucial in controlling the motility of ameboid cells Hoeller et al. (2016).

Cell traction forces are essential for cell motility, Trepat et al. (2009), and play an important role for cancer research where tumor cells typically exert stronger tractions than control, non-cancerous, cells Li et al. (2017); Koch et al. (2012). The distribution of cellular forces can be conveniently characterized using multipolar analysis Lauga and Powers (2009); Tanimoto and Sano (2014). In the absence of externally imposed forces, migrating cells are force-free bodies and the monopole term vanishes Tanimoto and Sano (2014). Dipolar forces have been observed in adherent cells both at the scale of the entire cell De et al. (2007); Schwarz and Safran (2013) and at subcellular scales Ghassemi et al. (2012). A dipolar force distribution is however front-back symmetric. Therefore, the dipole term alone is not enough to justify a preferred direction for persistent cell motion, which might be due to the existence of a quadrupolar force distribution Leoni and Sens (2015, 2017).

This theoretical consideration has been shown to hold true by experimental work Tanimoto and Sano (2014) which characterized motility of D. discoideum via traction force microscopy. These authors confirmed that both dipole and quadrupole terms provide robust tools for describing migrating cells – showing that while the axis along which the motion occurs is correlated with the orientation of the dipolar forces, the direction of motion is rather determined by the quadrupolar term along that axis.

The study of synchronization coupled to orientational dynamics (e.g. dipole and quadrupole terms) seems thus relevant for cell motility. So far, to our knowledge, only a few theoretical studies focused on the interplay of orientational order and synchronization. This was done in different contexts: studying the interplay of polar order and synchronization in soft active fluids Leoni and Liverpool (2014); considering swarmalators O’Keeffe et al. (2017), which combine the ability to swarm and to synchronize; investigating collections of self-propelling particles which can synchronize Levis et al. (2017, 2019). Experimentally, a work on bacterial suspensions interpreted collective oscillations as frequency synchronization of bacterial orientational dynamics Chen et al. (2017).

The implications of synchronization for cell motility have been examined in Leoni and Sens (2015), considering the interplay of intracellular and extracellular mechanical interactions for a model cell with dipolar forces Leoni and Sens (2015), already aligned, in contact with a fluid substrate. Such microscopic dipolar forces yield dipolar and quadrupolar terms at the cell scale. In the same spirit, another recent work Leoni and Sens (2017) investigated the role of stochastic adhesion in model cells with oscillatory force-distributions already aligned and synchronized.

Here we pursue the study of synchronization and cell motility considering also orientational order. We adopt an approach which generalizes the use of conserved and broken symmetry variables Chaikin and Lubensky (2000) to characterize the active non-equilibrium dynamics Marchetti et al. (2013) at scales larger than individual microscopic elements. For motile cells, such elements might be acto-myosin assemblies forming micron-sized sarcomeres Wolfenson et al. (2016).

The one particle concentration

[TABLE]

is associated to the probability of finding an active element with position $\mathbf{x}$, orientation $\hat{\bm{u}}$, phase $\phi$ at time $t$ given the microscopic dynamics of active elements with positions $\mathbf{x}_{n}(t)$, orientations $\hat{\bm{u}}_{n}(t)$ and phases $\phi_{n}(t)$, for $n=1,\ldots N$. The large scale behavior of active matter systems, Marchetti et al. (2013), can be characterized using vector and tensor quantities that are related to the vector character (via $\hat{\bm{u}}$) of the microscopic elements. In two dimensions, density, $\rho$, polarization, $\mathbf{p}$, and nematic orientation tensor, $\mathbb{S}$, are written as moments of the concentration as follows

[TABLE]

Here $\rho$ is a conserved quantity while $\mathbf{p}$ and $\mathbb{S}$ are broken symmetry variables: order parameters describing the degree of polar order or of nematic order in the system.

In systems displaying synchronization, however, the phase $\phi$ of each element provides an additional degree of freedom which is not taken into account in the above system of moments. To incorporate the effect of the phase, another series of moments generalizes the above construction as follows

[TABLE]

The first two terms in Eq.(3) were already obtained in Leoni and Liverpool (2014): $\Phi(\mathbf{x},t)$, in the spatially homogeneous limit, is the Kuramoto’s order parameter (Acebrón et al., 2005) and $\bm{\Pi}(\mathbf{x},t)$ describes the combined effect of vectorial symmetry and synchronization. The last moment, $\Sigma(\mathbf{x},t)$, generalizes the nematic orientation tensor $\mathbb{S}$ in presence of synchronization just like $\bm{\Pi}(\mathbf{x},t)$ generalizes $\mathbf{p}(\mathbf{x},t)$.

There is however a crucial difference between $\bm{\Pi}(\mathbf{x},t)$ and $\Sigma(\mathbf{x},t)$: the former can be employed to describe collective dynamics that has no global constraints, like in the case of suspensions of micro-swimmers where $\hat{\bm{u}}$ in Eq.(3) is associated to the orientation (or speed) of the individual swimmer; the latter is suitable for describing global constraints, like in the case of traction forces exerted by migrating cells where $\hat{\bm{u}}$ represents the director of the cell traction vector field which is constrained to be force-free. We note that $\Sigma(\mathbf{x},t)$ could be exploited also in other contexts, like the dynamics of chromatin inside nuclei Saintillan et al. (2018) where specific enzymes, intervening during chromatin remodeling, exert local active dipolar forces Bruinsma R. and et al. (2014).

No global constraints

Using Eq.(1), one can write $\bm{\Pi}(\mathbf{x},t)$ as

[TABLE]

Eq.(4) shows that in the special case of spatially homogeneous polar state, where all the directors point in the same direction, say $\hat{\bm{u}}_{n}\equiv\hat{\bm{U}}$ $\forall n$, then $\mathbf{p}\equiv\hat{\bm{U}}$ and the complex vector order parameter $\bm{\Pi}(\mathbf{x},t)$ is just the product of the two order parameters, one describing polar order and the other synchronization, as $\bm{\Pi}\equiv\mathbf{p}\Phi$. The same holds true for a synchronized state, where $\phi_{n}\equiv\varphi$, $\forall n$, and again $\bm{\Pi}(\mathbf{x},t)$ is the product of $\mathbf{p}$ and $\Phi$. However there are states where $\bm{\Pi}(\mathbf{x},t)$ describes non-trivial configurations. For example, one can construct states where both $\mathbf{p}=0$ and $\Phi=0$ and yet $\bm{\Pi}(\mathbf{x},t)\neq 0$ as it was pointed out in Leoni and Liverpool (2014).

Force-free constraint

A migrating cell is a force-free system, meaning that the forces exerted by the cell add up to zero, as experimentally verified e.g. in Tanimoto and Sano (2014). As it will be shown below, this constraint lead us to consider $\Sigma(\mathbf{x},t)$ for describing cell motility. Using Eq.(1), one can relate $\Sigma(\mathbf{x},t)$ to the microscopic variables as

[TABLE]

In earlier phenomenological studies on active systems which do not include the internal cyclic dynamics, Aditi Simha and Ramaswamy (2002), it was already noted that the force density of force-free active particles has nematic symmetry. Eq.(5) recovers that result in the case where the cyclic elements have all the same phases, $\phi_{n}=\phi_{0}(t)$, $\forall n$. Then, from Eq.(5)

[TABLE]

where $\Phi=e^{i\phi_{0}(t)}$. In particular, we can recover the regime discussed in Aditi Simha and Ramaswamy (2002), and in other following works on active matter Marchetti et al. (2013), by posing $\phi_{0}(t)=const$. This means that there is no cyclic dynamics and $\Phi$ is just a normalization factor in front of the nematic orientation tensor order parameter, $\mathbb{S}$.

The article is organized as follows: in Section .2 we introduce a microscopic model of the cell coupled to the substrate in Section .3, relying on that model, we study the cell migration speed. We use analytical methods to show how the speed relates to the broken symmetry variables Eq.(3) and investigate the speed numerically for simple configurations. Finally, in Section .4 we present the dynamic equations for Eq.(3) discussing the links with motility.

.2 Model cell on a fluid substrate

In this section we provide a concrete realization where, by construction, the vector field associated to the active traction forces exerted by a cell satisfies the force-free constraint: a set of dipolar units distributed in a quasi 2D region of space describing the cell boundary, see Fig.1.

Experimental measurements indicate the existence of sub-cellular contractile units in adherent cells, which can be detected and measured using micro-pillars, see e.g. Ghassemi et al. (2012); Wolfenson et al. (2016). The size of these contractile units is of the order of a few microns and it has been suggested that they are acto-myosin micro-sarcomeres Wolfenson et al. (2016). Similar sub-cellular traction units are seen also in migrating cells Jiang et al. (2015) and can be visualized by means of fluorescence Godeau et al. (2018). The intensity of the fluorescent signal might be a measure of the intensity of the active force exerted by a given cell, useful for comparisons with theoretical models.

.2.1 Traction force unit

The traction force unit, labeled with $\alpha$, is a dipolar element made of two particles, in turn labeled with $n=1,2$. Each particle is driven by an active force $\mathbf{f}^{\alpha}_{n}$. We write $\mathbf{f}^{\alpha}_{n}=(-)^{n+1}\mathbf{F}^{\alpha}$ and $\mathbf{F}^{\alpha}=F^{\alpha}\hat{\bm{u}}^{\alpha}$. For now, we pose $F^{\alpha}=f_{0}\cos\phi_{\alpha}$. As by construction $\hat{\bm{u}}^{\alpha}_{1}=-\hat{\bm{u}}^{\alpha}_{2}$, see Fig.1, the total force of the dipolar unit is zero, $\mathbf{F}^{\alpha}=\mathbf{f}^{\alpha}_{1}+\mathbf{f}^{\alpha}_{2}=f_{0}\cos\phi_{\alpha}(\hat{\bm{u}}^{\alpha}_{1}+\hat{\bm{u}}^{\alpha}_{2})=0$. The same remains true when summing all the forces at the level of the entire cell. The cell is thus force-free.

.2.2 Microscopic force-density

The force density $\bm{\mathcal{F}}(\mathbf{x})$ on the surrounding medium is obtained considering the microscopic density function, given by the Dirac’s delta ‘function’,

[TABLE]

This force distribution is what enters at r.h.s. of the continuum equation describing the medium, see e.g. Aditi Simha and Ramaswamy (2002); Marchetti et al. (2013) and Eq.(8) below.

To gain insight, we expand the delta function Aditi Simha and Ramaswamy (2002), writing $\mathbf{x}^{\alpha}_{n}=\mathbf{c}^{\alpha}+\frac{L^{\alpha}}{2}\hat{\bm{u}}^{\alpha}_{n}$ where $\mathbf{c}^{\alpha}$ is the coordinate of the center of the traction unit $\alpha$. Note that the first term of the expansion generates a contribution of the form $\sum_{\alpha}\cos\phi_{\alpha}\delta(\mathbf{x}-\mathbf{c}^{\alpha})\sum_{n}\hat{\bm{u}}^{\alpha}_{n}=0$ which vanishes due to the force-free condition. This equation can be translated in terms of the complex vector order parameter, $\bm{\Pi}$, compactly as $Re[\bm{\Pi}^{0}]=0$ where the over-script [math] denotes a spatially homogeneous value. Hence, for spatially homogeneous states, $\bm{\Pi}^{0}$ can be associated to the monopolar term of the traction forces.

The next term in the expansion is obtained by noting that $(\hat{\bm{u}}^{\alpha}_{1}\hat{\bm{u}}^{\alpha}_{1}+\hat{\bm{u}}^{\alpha}_{2}\hat{\bm{u}}^{\alpha}_{2})=2\hat{\bm{u}}^{\alpha}_{1}\hat{\bm{u}}^{\alpha}_{1}$. We further simplify our notation by posing $\hat{\bm{u}}^{\alpha}:=\hat{\bm{u}}^{\alpha}_{1}$. As a result we are left with $\bm{\mathcal{F}}(\mathbf{x})\approx f_{0}\sum_{\alpha}L^{\alpha}\cos\phi_{\alpha}\hat{\bm{u}}^{\alpha}(\hat{u}^{\alpha})_{k}\nabla_{k}\delta(\mathbf{x}-\mathbf{c}^{\alpha}).$ If $L^{\alpha}$ performs small oscillations around $l_{0}$ we obtain

[TABLE]

where $Re[...]$ indicates the real part of the term in parenthesis. As anticipated, this generalizes the results obtained in Aditi Simha and Ramaswamy (2002) to the case of particles which possess internal cyclic dynamics. Here $f_{0}l_{0}$ quantifies the magnitude of the microscopic dipole while $\Sigma_{ik}$ incorporates both orientations and phases.

.2.3 Cell coupled with the extra-cellular medium

In order to migrate, a cell needs to exert forces on its surroundings. In vivo, cells are surrounded by an extra-cellular matrix, a polymer gel, that can be reproduced also in vitro Muller et al. (2013); Cukierman et al. (2001). Here we discuss the viscoelastic dynamics of the gel in the long time limit Kruse et al. (2006) where the gel is pictured as a viscous fluid.

Although the fluid model is used here as a simple example of extra-cellular medium, this approach is not unrealistic: a fluid description of an elastic substrate, directly compared with experiments, was already given in Peterson (1996). Furthermore, in elastic substrates the force-distribution determines substrate deformations. Deformations alone would suffice for describing the dynamics of adherent cells but they are not enough for migration. In fact, for migrating cells one has to include also the adhesion dynamics to describe how a cell unbinds and subsequently rebinds at different positions on the substrate. A model of transient, stochastic, adhesion coupled to cell forces and mechanics was proposed in Leoni and Sens (2017). That framework can be generalized to a soft, deformable, substrate Leoni and Sens (2018) but is not considered here.

The equation describing an incompressible viscous fluid of viscosity $\eta$ and velocity $\mathbf{v}$ in the absence of inertia is

[TABLE]

This equation determines the velocity $\mathbf{v}(\mathbf{x})$ generated in the surrounding viscous medium, given the force distribution $\bm{\mathcal{F}}(\mathbf{x})$. In turn, one can express $\mathbf{v}$ as an integral of $\bm{\mathcal{F}}$ using Green’s functions formalism Doi and Edwards (1986). The component $v_{i}$ of the velocity is given by

[TABLE]

(where we use Einstein’s summation convention on repeated indexes for the component $j$). This allows us to connect the theory to the experimental measurements on the velocity of the substrate. For a semi-infinite substrate, with flat surface, $H_{ij}(\mathbf{r})\sim\frac{\delta_{ij}}{2\pi\eta r}$, see Leoni and Sens (2015).

The dynamic equation of a cell in contact with the medium is obtained from force balance ($2^{nd}$ Newton’s law). With the choice of fluid medium, and neglecting inertia, the equation is

[TABLE]

where $n=1\ldots N$ and $\alpha=1,\ldots M$ label different particles representing cellular adhesion sites Leoni and Sens (2015, 2017). $\gamma$ here is a friction coefficient (for a semi-infinite substrate is related to the Stokes’ drag, $\gamma=3\pi\eta a$ where $a$ is the radius of the particle Leoni and Sens (2015)) and $\mathbf{v}(\mathbf{x}^{\alpha}_{n})$ is the flow generated at position $\mathbf{x}^{\alpha}_{n}$ due to all the remaining particles, which can be computed using Eq.(9).

At this point there is still one thing to specify, namely the dynamics of the active forces (or equivalently, the dynamics of the phases $\phi_{\alpha}$ as forces depend on these variables). We will address this in the following.

.2.4 Dynamics of the force generators

The traction force elements are subjected to active forces and coupled with the substrate (here a fluid). The term “active” here means that forces vary in time. For cyclical variations, the time-dependence can be specified using a phase variable. Once this is set, from Eq.(10) we derive the dynamics of other variables such as orientations, amplitudes and phases of oscillations. The dynamic equations for such variables are needed to derive the equations for the moments, Eq.(3), describing the large scale, coarse-grained, dynamics.

Orientational dynamics

By construction, in our model, there is no torque on the traction units. Hence the traction units cannot rotate when isolated but they can rotate thanks to mechanical interactions with other units. The angular speed ${\omega}^{\alpha}$ of such rotations can be computed from $\dot{\hat{\bm{u}}}^{\alpha}$ via $\dot{\mathbf{L}}^{\alpha}=d(L^{\alpha}\hat{\bm{u}}^{\alpha})/dt$, see Fig.1(b), as ${\omega}^{\alpha}=\frac{\hat{\bm{u}}^{\alpha}}{L^{\alpha}}\wedge(\dot{\mathbf{x}}^{\alpha}_{1}-\dot{\mathbf{x}}^{\alpha}_{2})$. Using Eq.(10) we obtain

[TABLE]

Here $\mathcal{O}^{\alpha\beta}_{i}$ follows from Eq.(9). Its detailed expression is reported in the appendix .6. Note that the orientational dynamics contributes to the dynamic equation for the moments $\mathbf{P},\mathbb{S},\bm{\Pi},\Sigma$ (3). However for $\Sigma$, which is one of the relevant quantities, the contribution of $\omega^{\alpha}$ is sub-dominant compared to that of the phase dynamics, as discussed below. Moreover, in the study of the cell speed done below, in .3, for simplicity we consider configurations where the orientations $\hat{\bm{u}}^{\alpha}$ lie along a given direction and we neglect their dynamics.

To derive the equation for the amplitude and phase describing oscillatory dynamics we first need to obtain the equations regulating the deformation and the forces of the traction units.

Deformation dynamics of a traction unit

The deformation dynamics of each dipolar traction unit follows from the definition of $\mathbf{L}^{\alpha}={\mathbf{x}}^{\alpha}_{1}-{\mathbf{x}}^{\alpha}_{2}$, see Fig.1(b), as $\dot{L}^{\alpha}=\hat{\bm{u}}^{\alpha}\cdot(\dot{\mathbf{x}}^{\alpha}_{1}-\dot{\mathbf{x}}^{\alpha}_{2})$. Hence, the dynamic equation for the internal deformation of the element $\alpha$ is obtained from Eq.(10) as

[TABLE]

Here, $\bm{\mathcal{I}}^{\alpha\beta}$ results from Eq.(9). Its detailed expression is reported in the appendix .6. To leading order, $\dot{L}^{\alpha}\approx 2F^{\alpha}/\gamma$ has solutions $F^{\alpha}=-f_{0}\sin(\omega_{0}t+\phi_{\alpha})$ and $L^{\alpha}=2\frac{f_{0}}{\gamma\omega_{0}}\cos(\omega_{0}t+\phi_{\alpha})$. Thus, neglecting mechanical interactions $\mathcal{I}^{\alpha\beta}_{i}$ and in the absence of noise, the oscillating dipoles will maintain their relative phase difference which is controlled solely by the initial conditions. More interestingly in presence of interactions the phases can vary. To describe how the forces, and hence the phases, evolve we use a generic model of self-sustained oscillating forces introduced in Leoni and Liverpool (2012).

Force dynamics of a traction unit

The evolution of the forces, following Leoni and Liverpool (2012), is described by

[TABLE]

where $L^{\alpha}=l_{0}+d_{\alpha}$. As for the other parameters in Eq.(13), $\mathcal{K}>0$ determines the frequency of oscillations, $\mathcal{M}>0$ yields self-sustained oscillations, and $\mathcal{A}$, which can be either positive or negative, determines the non-isochrony of the oscillations. $\mathcal{S}$, associated to the saturation of oscillation amplitude, will be set equal to 1 in the following. As explained in the appendix .6, we can map Eq.(13) and Eq.(12) onto equations for the phase for the amplitude.

Amplitude and phase dynamics

Complex amplitudes are related to $L^{\alpha}$ via $L^{\alpha}\sim l_{0}+[A_{\alpha}e^{i\omega t}+c.c.]/2$. Similarly for the force, we pose $F^{\alpha}\sim i\gamma\omega_{0}[A_{\alpha}e^{i\omega t}-c.c.]/4$. In turn, the complex amplitude is related to real amplitude and phase $R_{\alpha}$, $\phi_{\alpha}$ via $A_{\alpha}=R_{\alpha}e^{i\phi_{\alpha}}$. Using these relations, we obtain an equation for the real amplitude

[TABLE]

with $\lambda=\frac{\mathcal{M}}{\gamma}$; $\beta=\frac{\mathcal{M}}{4\gamma}$ and $T_{ijkl}(x_{\alpha}-x_{\beta})=\nabla_{k}\nabla_{l}H_{ij}(x_{\beta}-x_{\alpha})$.

We consider interactions acting as a small perturbation to the non-interacting dynamics. In this regime, Eq.(14) describes deviations from a fixed point dynamics. The fixed point is the limit cycle ${R}_{\alpha}\sim R_{0}$ with $R_{0}=\sqrt{\lambda/\beta}.$ To study how these perturbations evolve, we write ${R}^{\alpha}\sim R_{0}+\delta{R}^{\alpha}$ and substitute this expression in Eq.(14) keeping terms up to order $\delta$.

The interpretation of the dynamics of $\delta R_{\alpha}$ has been given elsewhere for oscillators moving along one dimensions Leoni and Liverpool (2012). We summarize the main points:

1. due to interactions, the oscillators’ trajectories in the phase space move away from the limit cycle. This is described by terms $\delta R_{\alpha}$;
2. the limit cycle is a stable fixed point: deviations $\delta R_{\alpha}$ relax to the limit cycle and behave as dumped fluctuations. The result is essentially the same here with the main change due to the presence of the directors $\hat{u}^{\alpha}_{k}$ (as motions are no longer in one dimension).

Setting $\dot{\delta R}_{\alpha}=0$ we can eliminate $\delta R_{\alpha}$ in favor of the phase, see appendix .6, obtaining the equation for the phase dynamics

[TABLE]

where $\Delta=-\frac{\omega_{0}}{2}+\frac{\mathcal{K}}{\gamma\omega_{0}}$ and $\chi=\frac{3\mathcal{A}}{4\gamma\omega_{0}}.$ We note that, in particular, the sinusoidal term of Eq.(15) is responsible for synchronization Leoni and Liverpool (2012).

.3 Cell speed

From the mean position of the dipolar units ${\mathbf{C}}=\sum_{n}\sum_{\alpha}{\mathbf{x}}^{\alpha}_{n}/M$ (for $\alpha=1,\ldots M$) we obtain an important readout for cell motility: the cell speed

[TABLE]

which can be computed using Eq.(10). The choice of a quasi-2D model, see Fig.1, simplifies the calculation. In fact, we can write $\mathbf{x}^{\alpha}_{n}=x_{\alpha}\hat{\bm{x}}+(-)^{n+1}\frac{L^{\alpha}}{2}\hat{\bm{u}}^{\alpha}$ and decompose the problem along two independent directions, $\hat{\bm{x}}$ and $\hat{\bm{u}}^{\alpha}$. Performing an expansion valid at distances larger than $L^{\alpha}$ (consistent with our long wavelength description) the tensorial part $H_{ij}[(x_{\alpha}-x_{\beta})\hat{\bm{x}}]$ only depends on the director $\hat{\bm{x}}$ while the other variables are defined in the whole $(x,y)$ plane. We note also that the exact form of $H$ does not play a major role for what concerns the purpose of this work. In fact, due to the quasi-2D nature of the forces, where the centers $x_{\alpha}$ of the force-traction units lie along the x-axis, see Fig.1, different choices of $H$ lead to quantitative but not qualitative differences.

.3.1 Analytical study

To make progress, we need to expand the term $H$ describing interactions. Consistently with the analysis in .2.4 we shall assume that dipolar forces and deformations depend on time as $F^{\beta}=-f_{0}\sin(\omega_{0}t+\phi_{\alpha})$ and $L^{\beta}=l_{0}+R_{0}\cos(\omega_{0}t+\phi_{\alpha})$. The expansion produces terms of the form $F^{\beta}L^{\beta}\sim\sin 2(\omega t+\phi_{\beta})$ which average to zero over a period $T=2\pi/\omega$. However, terms containing $F^{\beta}L^{\alpha}$ give a finite contribution for $\alpha\neq\beta$. Using a mean-field approach Leoni and Liverpool (2010) the speed is given by

[TABLE]

Here the term $U_{ijklu}(x_{\alpha}-x_{\beta})=\nabla_{k}\nabla_{l}\nabla_{u}H_{ij}[(x_{\alpha}-x_{\beta})\hat{\bm{x}}]$ describes the dominant contribution of the interactions between different dipolar units. Along with the interactions, the speed is determined by the term in parenthesis in Eq.(17). The functional $\mathcal{G}_{jklu}[\Sigma,\Phi]$ is defined as

[TABLE]

where we have used the shorthand $\Sigma^{\beta}_{uj}=\Sigma_{uj}(x_{\beta})$ and similarly for $\Phi^{\alpha}$. $\{u\leftrightarrow l\}$ indicate a contribution similar to the previous expression in parenthesis but with the indexes $u$ and $l$ exchanged.

Eq.(17) and Eq.(18) confirm the need of a phase difference along the cell. Since the speed is controlled by $\Sigma$ and $\Phi$, which in turn are functions of two microscopic degrees of freedom $\hat{\bm{u}}$ and $\phi$, it is instructive to examine the impact of these variables separately.

Synchronized and spatially homogeneous states

We begin with examining the case of the phases which are identical everywhere. This means that $\Sigma^{\beta}_{uj}=\Phi^{0}\mathbb{S}^{\beta}_{uj}$ and

[TABLE]

Likewise, terms of the form

[TABLE]

so the associated speed in Eq.(17) is zero. The same conclusion is reached in the (less general) case of spatially homogeneous order where, in addition to $\Phi^{\beta}=\Phi^{\alpha}=\Phi^{0}$, one has also $\Sigma^{\beta}_{uj}=\Sigma^{\alpha}_{uj}=\Sigma^{(0)}_{uj}$.

Non-synchronized states

Let us now consider instead the alternative scenario where the orientations are all constrained (e.g. along $\hat{\bm{x}}$ – see also numerics below) while the phases are free to vary and may synchronize. We pose $\Sigma^{\beta}_{uj}=\mathbb{S}^{0}_{uj}\Phi^{\beta}$ (with e.g. $\mathbb{S}^{0}_{uj}=(\frac{\delta_{uj}}{2}-\hat{x}_{u}\hat{x}_{j})$). In this case

[TABLE]

and a non-zero cell speed in this case is associated to inhomogeneities of the order parameter $\Phi$. In the particular case of just two traction units oriented along the $\hat{\bm{x}}$ axis we recover what expected from previous analysis, Leoni and Sens (2015), for a cell made up of two units $A$ and $B$ where cell speed scales as $\dot{C}\sim\nabla^{3}H(r)\sin(\phi_{A}-\phi_{B}).$

Wave-like perturbations

To gain insight, we consider a wave propagating along $\hat{\bm{x}}$ and consider the effect of such a perturbation on the speed retaining terms up to linear order. Writing $\delta\Sigma^{\alpha}_{kl}=\mathbb{S}^{(0)}_{kl}[e^{i(\Omega t-\mathbf{k}\cdot\hat{\bm{x}}x_{\alpha})}-1]$ that is $\delta\Sigma^{\alpha}_{kl}\approx i\mathbb{S}^{(0)}_{kl}(\Omega t-\mathbf{k}\cdot\hat{\bm{x}}x_{\alpha})$ then $(\delta\Sigma^{\alpha}_{kl})^{*}\approx-\delta\Sigma^{\alpha}_{kl}$. Similarly, we pose $\delta\Phi^{\alpha}\approx i\rho^{(0)}(\Omega t-\mathbf{k}\cdot\hat{\bm{x}}x_{\alpha})$. We obtain as a result the net average speed

[TABLE]

where $\mathcal{G}^{(0)}_{jklu}=[2\mathbb{S}^{(0)}_{kl}\mathbb{S}^{(0)}_{uj}+4\rho^{(0)}\delta_{kl}\mathbb{S}^{(0)}_{uj}]+[l\leftrightarrow u].$ The motion is along $\hat{\bm{x}}$ here but it could be more complex for other choices of $\mathbb{S}^{(0)}_{kl}$. Thus the cell achieves motion by propagating waves along the cell boundary. A wave breaks time reversal symmetry Lauga and Powers (2009) and is consistent with the requirements for swimming at low Reynolds number Purcell (1977).

.3.2 Numerical study

Results analogous to the conclusions of our model were reported in an experimental study Hoeller et al. (2016) where cell motility was related to cellular oscillations and their synchronization, see also Corallino et al. (2018). In Hoeller et al. (2016), actin oscillators distributed along the periphery of ameboid cells were observed. By intervening on biochemical regulators, the authors have been able to induce strong coupling among the oscillators which led to in-phase synchronization along the cell periphery. In this case, as also expected from our analysis, migration is suppressed. For normal untreated cells, migration (and response to chemical gradients) was instead observed.

In Fig.2 we show an example (obtained from our numerics – see appendix .7 for details) reproducing the phenomenology observed in experiments of ref. Hoeller et al. (2016): we consider $N=33$ oscillators along a line and study their collective behavior. In (a) oscillators are synchronized in-phase. In (b) oscillators are synchronized in anti-phase. We have computed the cell speed numerically from Eq.(16), and studied the speed as a function of the non-isochrony parameter $\mathcal{A}$ for different numbers of oscillators, i.e. different numbers of traction units. We find that when particles are synchronized in anti-phase, $\mathcal{A}>0$ in (13), then the cell speed is higher than what observed in case of in-phase synchronization, $\mathcal{A}<0$. Although our analytical expression for the migration speed in (17) indicates that the cell should not move in case of spatially homogeneous synchronization, the small motion observed in the numerics might be due to inhomogeneities and instabilities of the ordered phase. Instead, anti-phase synchronization leads to the propagation of waves Leoni and Liverpool (2012), which are associated to a finite cell speed as in Eq.(20). The reason why cell speed has a bi-phasic behavior as a function of the number of traction units, in Fig.3, can be understood using a scaling argument reported in appendix .7.

.4 Dynamic equations for $\Phi$ and $\Sigma$

In the previous section we have seen that the cell speed, the important readout for cell motility, depends on certain moments, defined in Eq.(3). It becomes then important to discuss if, and under which conditions, these moments can attain finite values needed for motility. The growth of the moments is related to a change of sign of a parameter in their dynamic equations which signals the transition to order. Here we report the dynamic equations for these moments and we find that mechanical interactions, propagated via the substrate, among the traction units can promote such a transition.

The concentration $c$, in Eq.(1), satisfies a Smoluchowski equation of the form of Eq.(25) from which we derive equations for the moments Eq.(2) and Eq.(3), see appendix .8 for the details. We find that the relevant moments are $\Phi$ and $\Sigma$. Density $\rho$ has a trivial equation due to our neglect of spatial gradients. We also find that the dynamics of $\bm{\Pi}$ shows no transition to order and is slave to that of $\Sigma$. All this is expected as the microscopic model is made up of non-polar elements.

The dynamic equation for $\Phi$ is

[TABLE]

Here $[k\leftrightarrow l]$ indicates the same contribution as in the previous parenthesis but exchanging $k$ and $l$. $D_{\phi}$ describes phase diffusion. Again $\Phi^{\alpha}=\Phi(x_{\alpha},t)$ and we defined $\mathcal{Z}=\frac{l^{2}_{0}\omega_{0}\gamma}{8\beta}(1+i\frac{\beta}{\chi})$. Note that in the aligned case, where $\hat{\bm{u}}^{\alpha}=\hat{\bm{x}}$, the equation simplifies recovering that of oscillators in 1D Leoni and Liverpool (2012) although with some differences due to the nature of the oscillators (dipolar here, while in Leoni and Liverpool (2012) the oscillators were monopoles).

The dynamic equation for the complex nematic tensor reads

[TABLE]

where $D$ describes rotational diffusion, $perm_{fgil}=\delta_{fi}\mathbb{S}_{gl}+\delta_{fl}\mathbb{S}_{ig}+\delta_{ig}\mathbb{S}_{fl}+\delta_{lg}\mathbb{S}_{if}$ and again $[k\leftrightarrow l]$ indicate the same contribution as in the previous parenthesis but with indexes $k$ and $l$ exchanged.

Both Eq.(22) and Eq.(21) depend on the nematic orientation tensor $\mathbb{S}$. Its dynamic equation, which is rather cumbersome and therefore reported in the appendix .8.4, depends on $\Phi$ and $\Sigma$. Setting $\partial_{t}\mathbb{S}_{fg}=0$ one can in principle obtain $\mathbb{S}$ as function of $\Phi$ and $\Sigma$. That, in turn, allows one to obtain closed form equations for Eq.(21) and Eq.(22) which would describe the coupled dynamics of $\Phi$ and $\Sigma$. Instead of analyzing this general case, which involves lengthy and approximate calculations, we would like to summarize the generic aspects which are relevant for the application of this theory to the motility of cells.

First, considering spatially homogeneous quantities, the term $\rho^{0}\Phi^{0}$ in Eq.(21) is responsible for the transition to order, signaled by exponential growth of the order parameter $\Phi^{0}$. That is possible when the non-isochrony parameter of Eq.(13) satisfies $\mathcal{A}<0$, i.e. $\chi<0$. However, as we discussed in .3.1, spatially homogeneous values of $\Phi$ are not enough for motility. A finite cell speed, instead, is associated to wave-like perturbations, see Eq.(20). These could arise as instabilities of the ordered state, as was discussed in Leoni and Liverpool (2014). Alternatively, wave like behavior of the order parameter can be obtained in the opposite regime where $\mathcal{A}>0$, i.e. $\chi>0$, and the system’s synchronization is (close to) anti-phase as our numerical study indicates, see Fig.2.

Secondly, homogeneous order in Eq.(22) is possible thanks to the presence of a term $\rho^{0}\Sigma^{0}_{jk}$. The transition to order (signaled by exponential growth of the order parameter $\Sigma^{0}$) occurs at some finite value of the density $\rho^{0}$ and requires $\chi<0$, hence controlled by the sign of the non-isochrony parameter. This is interesting as it shows that the dynamic behavior of this novel complex order parameter $\Sigma$ is to some extent related to that of $\Phi$. Also here a spatially homogeneous value of $\Sigma^{0}$ is not enough for motility; wave-like perturbations, as examined in Eq.(20), might arise as instabilities of the ordered state Leoni and Liverpool (2014) or be associated to other kinds of order (e.g. the analogue of anti-phase for this model).

.5 Discussion

In this article we have developed a theoretical formalism that relies on soft active matter and synchronization to gain insight on the motility of cells. The theory combines vectorial degrees of freedom associated to the direction of the cell forces and a phase variable associated to the time-dependent character of cell forces. The force-free constraint, saying that cellular forces add up to zero, led us to introduce a novel broken symmetry variable - a complex tensor field - to model the distribution of cellular traction force.

We computed the cell speed and found that it also depends on this tensor field as well as on the complex order parameter describing synchronization in Kuramoto model. We found that in-phase synchronization is not enough to promote motility. Our model requires more complex patterns, e.g. the propagation of waves. These might arise from

1. hydrodynamic instabilities of the ordered synchronized state or 2) instabilities of the anti-phase state. We have performed numerical analysis which confirms this view and suggests that both mechanisms are possible, with 2) giving higher speed.

To make further progress, and test our predictions regarding the correlation between cell speed and its relation with $\Sigma$ and $\Phi$, it would be helpful to compute $\Sigma$ directly from experimental traction patterns of migrating cells. Note instead that for the case of adherent cells, which do not migrate, the complex tensor description is a priori not necessary.

Here we have provided a realization of cell with subcellular traction units that implement, locally, the force-free condition. This allows us to simplify the analysis and to derive general analytical expressions. Other types of force-distributions, where the force-free condition is not satisfied locally, but only recovered at the cell scale, are expected to give qualitatively similar results, with quantitative differences in Eq.(21) and Eq.(22).

Finally, here we have modeled the extra-cellular environment as a fluid substrate. A more realistic description requires to model an elastic, adhesive, substrate. Migration on such complex environments can be tackled relying on the formalism we have developed here, by modifying the force balance, Eq.(10), to include the adhesion’s dynamics Leoni and Sens (2017) and by extending that approach to elastic substrates. This task is left for a future work.

Acknowledgements.

The author acknowledges financial support from the ICAM Branch Contributions and Labex Celtisphybio No. ANR-10- LBX-0038 part of the IDEX PSL No. ANR-10-IDEX-0001- 02 PSL and many stimulating scientific discussions with T. B. Liverpool, M. C. Marchetti, and P. Sens.

APPENDIX

In this section we present some details which complement the main text.

.6 Dynamics of the force generations

Orientational dynamics

In Eq.(11)

[TABLE]

where $V_{ijkluv}(x_{\alpha}-x_{\beta})=\nabla_{v}\nabla_{u}\nabla_{k}\nabla_{l}H_{ij}[(x_{\alpha}-x_{\beta})\hat{\bm{x}}].$

Traction unit deformation dynamics

The dynamic equation for the internal deformation of the element $\alpha$ is given by Eq.(12). There, $\bm{\mathcal{I}}^{\alpha\beta}$ follows from force balance, Eq.(10), and to the leading order is given by

[TABLE]

where the tensorial part is $T_{ijkl}(x_{\alpha}-x_{\beta})=\nabla_{k}\nabla_{l}H_{ij}[(x_{\alpha}-x_{\beta})\hat{\bm{x}}]$ as above in the main text.

Amplitude and phase dynamics

The equations for amplitude and phase $\dot{R}_{\alpha}$ and $\dot{\phi}_{\alpha}$, describing the time-dependent oscillatory dynamics are obtained by mapping $L^{\alpha}$ and $F^{\alpha}$ onto amplitudes ${R}_{\alpha}$ and phases $\dot{\phi}_{\alpha}$. This is achieved by introducing complex amplitudes $A^{\alpha}$ related to the deformations and forces via $L^{\alpha}\sim l_{0}+[A_{\alpha}e^{i\omega t}+c.c.]/2$ and $F^{\alpha}\sim i\gamma\omega_{0}[A_{\alpha}e^{i\omega t}-c.c.]/4$, as explained in the main text. Thanks to this change of variable, we obtain dynamic equations for $\dot{A}_{\alpha}$ and $(\dot{A}_{\alpha})^{*}$. Finally, writing $A_{\alpha}=R_{\alpha}e^{i\phi_{\alpha}}$ we derive equations for the amplitudes from $\dot{R}_{\alpha}=\frac{e^{-i\omega t}\dot{A}_{\alpha}+e^{i\omega t}(\dot{A}_{\alpha})^{*}}{2}$ and for the phases from $\dot{\phi}_{\alpha}=\frac{e^{-i\omega t}\dot{A}_{\alpha}-e^{i\omega t}(\dot{A}_{\alpha})^{*}}{2iR_{\alpha}}$.

Setting $\dot{\delta R}_{\alpha}=0$ in Eq.(14) we can eliminate $\delta R_{\alpha}$ in favor of the phase variable obtaining

[TABLE]

Using this expression, we get Eq.(15)

.7 Cell Speed

.7.1 Analytical study

In the equation Eq.(16): to leading order we find $\dot{C}\propto\frac{1}{M}\sum^{M}_{\alpha=1}\sum_{\beta\neq\alpha}\nabla^{3}(\frac{1}{r_{\alpha\beta}})L^{\alpha}\dot{d}^{\beta}.$ The term $\nabla^{3}(\frac{1}{r_{\alpha\beta}})$ scales as $\nabla^{3}(\frac{1}{r_{\alpha\beta}})\sim\frac{1}{r^{3}_{0}(\alpha-\beta)^{3}}$ where $r_{0}$ is the average separation between the centres of the traction units (here chosen to be constant), while the term $L^{\beta}\dot{d}^{\alpha}$ scales as $L^{\alpha}\dot{d}^{\beta}\sim\omega_{0}d^{2}_{0}\sin(\beta-\alpha)$ where $d_{0}$ and $\omega_{0}$ are the oscillation amplitude and frequency. The sum over the index $\alpha$ brings $M$ such contributions which cancel the $1/M$ pre-factor. To study the cell speed as a function of $M$ we can estimate the previous expression as

[TABLE]

For small values of $M,M\sim 1+\epsilon$ the speed $|\dot{C}|$ grows linearly with $\epsilon$. For large values of $M$, $|\dot{C}|$ saturates as $Si(M)\to\pi$ while the remaining terms in parenthesis tend to zero. Taken together these limits provide an explanation for the bi-phasic behavior seen in Fig.3.

.7.2 Numerical study

We integrate Eq.(12), Eq.(13), Eq.(16) using Euler scheme with time-step $dt=0.0031s$. Parameters characterizing the dynamics of individual oscillators are: $\mathcal{K}=1pN/(\mu m\quad s)$, $\mathcal{S}=1/(\mu m)^{2}$, $\mathcal{M}=0.01pN/\mu m$ and $\mathcal{A}=\pm 0.01pN/(\mu m^{3}\quad s)$. The drag coefficient is $\gamma=3\pi\eta a$ as in Leoni and Sens (2015) (ignoring for simplicity cell viscosity) where $a$ is the radius of each adhesion site, $a=0.1\mu m$. The viscosity of the substrate is chosen to be $\eta\sim 6.3\eta_{H_{2}O}$ where $\eta_{H_{2}O}$ is water’s viscosity and interactions are modelled as $H(r)=\frac{1}{2\pi\eta r}$, Leoni and Sens (2015) , with average separations among units $r_{0}=25\mu$m and equilibrium length $l_{0}=15\mu$m.

.8 Dynamic equations for the moments

The concentration $c$, in Eq.(1), satisfies a Smoluchowski equation of the form Marchetti et al. (2013)

[TABLE]

where $D$ and $D_{\phi}$ are respectively the rotational and the phase diffusion constants. Note that terms of the form $\bm{\nabla}\cdot[v_{0}\hat{\bm{u}}c]$, which are typical non-equilibrium terms describing self-propelling particles (with speed $v_{0}\hat{\bm{u}}$ ), are absent here as the individual elements do not self-propel.

The two terms in parenthesis $[...]$ at r.h.s. are respectively rotational and phase currents. Here $\bm{\mathcal{R}}$ is defined as $\bm{\mathcal{R}}:=\hat{\bm{u}}\wedge\frac{\partial}{\partial\hat{\bm{u}}}$, while $D_{\phi}\frac{\partial}{\partial\phi}$ describes phase diffusion. The other terms in parenthesis comprise

[TABLE]

and

[TABLE]

where $\bm{\omega}$ and $\psi$ are given by Eq.(11) and Eq.(15).

For the angular dynamics, we use the expression of Eq.(11) and get

[TABLE]

where $c.c.$ denotes the complex conjugate and $[v\leftrightarrow k]$ indicates the same contribution as the one in the previous parenthesis but with $v$ replaced by $k$ etc. We find that $\Omega_{a}$ is sub-dominant (in terms of powers of the inverse separation which controls the interaction strength) in the equations for the moments $\bm{\Pi}$ and $\Sigma$. In fact the dominant contribution comes from the term $\Psi$ associated to the term describing the two-body phase dynamics.

The phase dynamics is obtained from Eq.(15) as

[TABLE]

Below we present the dynamic equations for the moments, Eq.(2), Eq.(3). For example, the equation for the moment $\Phi$ is obtained from $\dot{\bm{\Phi}}=\int dxe^{i\phi}\partial_{t}c$ and by inserting, in this expression, the r.h.s. of Eq.(25). A similar procedure is followed for the other moments.

.8.1 Equation for $\rho$

Neglecting, as we did, terms containing the spatial gradients then the density equation is simply $\partial_{t}\rho=0.$

.8.2 Equation for $p$

The dynamic equation for $\mathbf{p}$ is

[TABLE]

where $\Gamma_{bklupf}=\delta_{bk}\Delta_{lupf}+\delta_{bl}\Delta_{kupf}+\delta_{bu}\Delta_{lkpf}+\delta_{bp}\Delta_{lukf}+\delta_{bf}\Delta_{lupk}$ and $\Delta_{abcd}=\delta_{ab}\delta_{cd}+\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}$. Hence, the dynamics of $\mathbf{p}$ is slave to that of other moments. Moreover, when the particles are synchronized, the interaction terms vanish. This is expected based on the symmetry of the microscopic elements which are nematic rather than polar. One can see that by considering that the quantity $\bm{\Pi}$, in the synchronized case, reads $\bm{\Pi}=\Phi^{0}\mathbf{p}$. Hence $\bm{\Pi}^{*}=(\Phi^{0})^{*}\mathbf{p}$ and terms of the form $\bm{\Pi}^{*}\Phi^{0}-c.c.$ vanish in that case. A similar consideration holds true also for terms of the form $\bm{\Pi}^{*}\cdot\Sigma$.

.8.3 Equation for $\Pi$

The equation for $\bm{\Pi}$ is obtained from $\dot{{\Pi}}_{f}=\int dx\hat{u}_{f}e^{i\phi}\partial_{t}c$. To get this equation we use the r.h.s. of Eq.(25) neglecting the term $c\bm{\Omega}$ which is subdominant (as powers of the inverse of the distance, which characterizes the coupling strength) compared to $c\Psi$. We obtain

[TABLE]

Note that also here the dynamics of $\bm{\Pi}$ is slave to the dynamics of $\Sigma$ and $\Phi$ and no transition to spatially homogeneous order is obtained. All this is expected as the microscopic force generators are non-polar elements.

.8.4 Equation for $\mathbb{S}$

The dynamic equation for $\mathbb{S}$ is given by

[TABLE]

which shows that the dynamics of $\mathbb{S}$ is slave to that of $\Sigma$ and $\Phi$. Here the term $perm_{pgbklu}=\Sigma^{*}_{bg}\Delta_{pklu}+\Sigma^{*}_{kg}\Delta_{bplu}+\Sigma^{*}_{lg}\Delta_{bkpu}+\Sigma^{*}_{ug}\Delta_{bklp}.$ Again, $\{f\leftrightarrow g\}$ represents the same contribution of the terms in the previous parenthesis with indexes $f$ and $g$ exchanged.