Cell migration plays a crucial role during many key biological processes, from morphogenesis to cancer progression. As a result, the molecular components involved in cell migration, most notably the actin cytoskeleton, have been intensively investigated. Despite the great progress that was made, it is still an open question how do the different cellular components self-organize in a spatial pattern that maintains the robust motile cell shape. Several theoretical models have been proposed to explain the spontaneous emergence of the motile cell shape.

When cells migrate within the blood and lymphatic vessels, they experience fluid flow, which exerts shear forces on the cells. Outstanding examples include lymphocytes in the lymphatic vessels (von Andrian and Mempel, 2003), neutrophils and T Cells (types of immune cells) rolling and migrating in the blood vessels towards a site of inflamation (Luster et al., 2005; Shulman et al., 2009), endothelial cells (Alghanem et al., 2021) and fibroblasts crawling to sites of injury (Kole et al., 2005). Therefore, an understanding of how shear stresses influence cell movement is essential, and is still lacking.

One reason for this is that there appear different responses to shear in different cells, and within the same cell type under different conditions, as shown by the following examples. The exposure of sparsely plated endothelial cells, or a wounded monolayer, to shear flow inhibits their migration against the flow (Zaidel-Bar et al., 2005). The direction of T-lymphocyte cell migration under shear flow depends on the adhesion receptors (Dominguez et al., 2015; Anderson et al., 2019): When VCAM-I (Vascular Adhesive Molecule-I) is used, the cells migrate with the flow, while the cells migrate against the flow when ICAM-I (Intracellular adhesive molecule-I) is used. As the shear rate increases, T-lymphocytes favour migration against the flow when ICAM-1 is present, even in the presence of VCAM-1 (Dominguez et al., 2015). The migration of T Cells with the shear also depends on previous exposure to the flow (Piechocka et al., 2021).

However, one prominent feature that appears consistently in many cell types, is a tendency to migrate up stream against the flow. This behavior was observed in T-lymphocyte cells (Steiner et al., 2010; Dominguez et al., 2015; Anderson et al., 2019), microvascular endothelial cells (Song et al., 2009; Ostrowski et al., 2014; Alghanem et al., 2021), circulating tumor cells (Follain et al., 2018), and in the single-celled amoeba Dictyostelium discoideum (Décavé et al., 2003; Fache et al., 2005; Dalous et al., 2008). The origin of this prevalent migration response to shear flow is not understood at present.

Here, we utilize a recently developed theoretical model to explore the response of adherent cells to shear flow. The coarse-grained theoretical model describes the shape dynamics of a vesicle that contains curved membrane proteins that recruit active protrusive forces from the cytoskeleton. This model was shown to give rise to spontaneous pattern formation on the membrane, resulting in different shapes of vesicle (Figure 1). In the presence of adhesion to an external substrate, a motile phenotype emerges in this model. Here we exert on this vesicle an external force field that emulates the viscous drag force due to the fluid flow. This is an approximate description, that avoids solving the full flow field around the cell, but may provide us with a qualitative understanding of the main physical effects of the shear forces due to the flow. The simplicity of the model makes the calculated results very general, not cell-type-specific. They may therefore shed light on basic physical processes that apply to many cell types, such as the observed tendency of many types of motile cells to migrate upstream.

Our model is based on Monte Carlo (MC) simulations to evolve the shape of a vesicle in time (see Supplementary Material). The vesicle is described by a three-dimensional surface (Figure 2A) of N vertices, each connected with bonds of length l, to form a closed, dynamically triangulated, self-avoiding network, with the topology of a sphere (Fošnarič et al., 2019). The position vector of the ith vertex is r i ⃗ , where i ∈ [1, N]. The vesicle contains mobile curved membrane complexes, which are also sites of force application, representing the protrusive force exerted by actin polymerization. The vesicle is placed on a flat adhesive surface parallel to the x-y plane at z = z _ad.

The total energy of the vesicle is the sum of various contributions (Sadhu et al., 2021): (a) the local bending energy due to the membrane curvature, (b) the energy due to binding between nearest-neighbour membrane protein complexes, (c) the energy due to the active cytoskeleton force, (d) the adhesive energy due to the attractive interaction between the vesicle and the substrate and (e) the energy due to the force experienced by the vesicle due to shear flow.

The bending energy is given by the Helfrich expression (Helfrich, 1973) as W b = κ 2 ∫ A C 1 + C 2 − C 0 2 d A (1) where C ₁ and C ₂ are the two principal curvatures, and κ = 20k _B T is the bending rigidity. We set the spontaneous (isotropic) curvature C 0 = 1 l min − 1 for the nodes that contain the curved protein complexes, while it is set to zero for the bare membrane. The percentage concentration of the vertices with curved protein is denoted as ρ. The interaction energy between nearest-neighbour proteins is expressed as W d = − ∑ i < j w H r 0 − r i j (2)

Where H is the Heaviside step function, r ₀ is the interaction range, r i j = | r i ⃗ − r j ⃗ | is the distance between proteins, and w = 1 k _B T is the interaction energy between neighboring proteins in all the simulations in this paper.

The energy (work) due to the active protrusive force exerted by actin polymerization at the positions of the curved protein complexes δ W F = − F ∑ i n i ̂ ⋅ δ r i ⃗ , (3) Where, n i ̂ is the outward normal to the membrane and index i runs over the positions of all proteins, F is the strength of the active force, and δ r i ⃗ is the MC shift in the position of the node. The total active force is denoted by F ^act (Figure 2A).

The vesicle can adhere to the adhesive surface located at z = z _ad, and this energy contribution is W ad = − ∫ A V z d A , (4) Where V(z) is the interaction potential between the adhesive surface and the vesicle. If the node is close to the surface, z _ad ≤ z _i ≤ z _ad + δz, then the adhesion energy is V(z) = E _ad, while it is zero for all other nodes. We set δz = l _min for this whole paper, which is the minimal permitted bond length (to prevent pathological triangulation). The adhesive surface acts as a rigid barrier that the membrane can not penetrate.

Within this model we do not explicitly describe the fluid surrounding and within the vesicle. The MC calculation does not describe the correct time-scale of shape changes, as the dissipative processes involving the fluid flow are not included. This model can predict the shape changes of the vesicle as it minimizes the energy terms listed above.

This limitation means that when we wish to add the effects of shear forces due to fluid flow, we have to implement some approximate way for exerting these forces. We consider a fluid flow that has a far-field linear profile close to the surface on which the vesicle is adhered (Figure 2A), in the direction v ̂ shear . We do not solve the exact flow field around the vesicle, but we assume that the force exerted on the vertices of the vesicle by the flow is everywhere tangential to the vesicle surface. This assumption is motivated by being in the regime of the low Reynolds number that is applicable for cells, where the non-linear (inertial) terms in the Navier-Stokes equation are neglected (Cantat and Misbah, 1999). In this regime the viscous drag term is tangential to the membrane surface, if the membrane is roughly stationary on the time-scale of changes in the flow.

This tangential force due to the flow is in the direction of the projection of the shear flow direction on the local tangent plane ( n ̂ × ( v ̂ shear × n ̂ ) ) , where n ̂ is the local outwards normal to the surface (Figure 2B). The force on the vertex due to the shear flow is given by F ⃗ shear = F shear n ̂ × v ̂ shear × n ̂ (5) Where the force magnitude due to the shear is assumed to be given by the linear far-field flow velocity at the corresponding distance from the adhesive surface F shear = a z − z ad = a Δ z (6) as shown in Figures 2A,C. Here, we assumed a linear shear profile of the flow speed as function of the distance from the adhesive surface (due to non-slip boundary condition), and a gives the shear rate parameter. It determines how fast the shear force increases with the height from the adhesion surface, so depends on both the imposed flow speed and the viscosity of the fluid.

The force due to the shear flow is applied on the nodes of the vesicle as an external force, which gives the following contribution to the energy (work) of the system due to each MC node move (similar to Eq. (3)) δ W s = − a ∑ i z i − z ad n ̂ × v ̂ shear × n ̂ ⋅ δ r i ⃗ . (7) Note that the shear force is applied to each node along the local tangent, which fluctuates due to local membrane shape undulations (inset of Figure 2C).

The total energy change of the system, per MC mode, is given by δ W = W b + W d + δ W F + W ad + δ W s . (8)

We first verified that our implementation of the effective force due to shear flow produces reasonable results, by calculating the cross-sectional shape of a protein-free vesicle. We found that the vesicle deforms, and tends to lift and detach from the adhesive surface, as the shear flow increases (Supplementary Figure S1), similar to the results of previous experimental and theoretical work (Cantat and Misbah, 1999). These results validate that our implementation of the flow-induced shear forces is qualitatively realistic, although it is not quantitatively accurate as it does not include solving the real flow field around the time-dependent vesicle shape.

Our approximation neglects an additional pressure term that can induce normal forces on the membrane, which is only significant near stagnation points of the flow, such as at the vesicle-surface contact line (Cantat and Misbah, 1999). Note that in three dimensions this local pressure at the stagnation points should have a lesser effect compared to the two-dimensional calculation that is given in (Cantat and Misbah, 1999). Nevertheless, despite neglecting this stagnation pressure effect, our vesicle deforms and lifts in a manner that is qualitatively very similar to the full hydrodynamic solution for the passive adhered vesicle (Cantat and Misbah, 1999) (Supplementary Figure S1). In addition, this effect of stagnation pressure at the cell-surface contact line should play a less important role for the highly spread-out motile vesicle which resembles the motile cell (Figures 2, 3).

Despite the approximate treatment of the shear-flow induced forces on the membrane, our tangential force (Eqs.(5) and (6) is very similar to the full numerical solution of shear Stokes flow over a two-dimensional (Gaver and Kute, 1998) and three-dimensional hemispherical hump (compare Figure 2A with (Wang, 2004)).

Adding passive curved proteins can enhance the adhesion of the vesicle, mitigating the tendency of the vesicle to lift and detach due to the shear flow (Supplementary Figure S1). Note that similar detachments were also observed for cells exposed to strong shear flows (Décavé et al., 2002). Other flow regimes can be implemented in our model, for example, the case of large “slip-length”, where the flow field becomes (Eq. (6)): F ^shear = a (z − z _ad) + b, and the fluid has a non-zero speed on the solid substrate (see, for example, (Zhu and Granick, 2002; Wang et al., 2021)). In Supplementary Figure S2A we show the cross-section shape of the passive vesicle for the extreme case of large slip-length, a = 0, b ≠ 0, which is very different from the case of the linear gradient regime (where a ≠ 0, b = 0).

We next investigated the response of our vesicle that contains active curved proteins to the shear flow.

We have previously found that our minimal-cell model can describe a variety of steady-state shapes of the adhered vesicle (Figure 1) (Sadhu et al., 2021). One such phenotype is a motile, crescent-shaped vesicle, which appears for strong adhesive interaction E _ad = 3k _B T and sufficiently large active protrusive forces (parameters indicated by the blue star in Figure 1). Such a motile vesicle has a direction of polarity, determined by the direction of the total active forces (F ^act), due to the local forces applied by the curved proteins that form the leading edge cluster (Figure 2A). This shape is self-sustaining, with the highly curved leading edge maintained by the active forces, thereby stabilizing the cluster of curved proteins that seek to minimize their bending energy at such high curvature.

In order to demonstrate the effects of shear flow on the vesicle shape, we showed plot in Figure 2C the cross-sectional shape of the vesicle for the case of no shear (blue line), shear in the positive x direction (parallel to the initial polarization of the vesicle, green line), and shear in the negative x direction (opposite to the initial polarization of the vesicle, red line). As shown in Figure 2C, when the shear flow acts opposite to the polarity of the vesicle, the two opposing forces tend to deform the vesicle and sharpen its leading edge. This should therefore further stabilize the leading edge cluster of curved proteins. On the contrary, when the shear is in the direction of the polarity, the leading edge becomes less curved, thereby destabilizing the leading edge cluster, as it cannot maintain the high curvature that minimizes the bending energy of the curved proteins. We therefore expect that the motile vesicle will respond and modify its motility due the shear forces. These deformations of the motile vesicle due to the shear-flow are very similar to those calculated using fully numerical solutions of the Stokes equations for adhered fluid droplets in shear-flow (Dimitrakopoulos and Higdon, 2001).

In Figure 3A we show the effects of different shear directions on the trajectories of the motile vesicle. If shear is absent, the motile vesicle moves persistently along its polarity direction, which meanders over time due to random fluctuations. When shear is parallel to the polarity, we find that the protein aggregate reorients and with it the migration path of the vesicle makes a U-turn, to end up facing the shear flow direction. If the shear is opposite to the polarity, the speed is greatly diminished due to the competition between the active force and the shear force, but the migration direction is highly stable, maintaining the upstream path. Finally, when the shear flow is perpendicular to the polarity, we find again that the protein aggregate rotates and reorients to face the shear.

To explain the origin of these responses of the migration to the shear, we investigate how the shear forces modify the shape of the vesicle (Figure 2C), and how these shape changes affect the leading-edge protein aggregate. The average mean-curvature c avg pro of all the vertices with the curved proteins is shown in Figure 3B. We can see a significant drop in the curvature of the leading edge at the early times for the case when the shear is initially parallel to the polarity direction of the vesicle (green solid line). This is a quantification of the effect shown in Figure 2C As the vesicle reorients to face the shear flow, the average mean-curvature c avg pro increases, and at long times, when the vesicle faces the shear, it is slightly higher in the presence of shear compared to the no-shear case. The proteins aggregate at the leading edge of the vesicle, and prefer the configuration with the higher curvature, which is oriented against the shear flow.

We follow this reorientation process in Figure 3C, and in Figure 3D we plot the local mean curvature of the proteins along the leading edge (c ^pro). The initial reduction of the curvature due to the shear is most significant in the direction of the shear (θ ≈ 0, blue line), as the vesicle is pushed from behind and the front fattens. As time evolves, the proteins rotate and the curvature along the leading edge increases, first in the direction facing the flow (red line, compare negative vs. positive angles). Finally, the protein aggregate orients in the direction facing the shear flow, and the middle of the protein aggregate has the highest curvature.

A very similar dynamics is observed for shear that is perpendicular to the initial vesicle polarization, as shown in Figure 3E. The vesicle experiences a shear force that pushes from one side, which makes the farthest side of the protein aggregate fatter. Therefore, the protein aggregate rotates towards the more highly curved region and becomes motile against the shear. The curvature at the sites of the proteins c ^pro is shown in Figure 3F. Initially the far side (θ ≈ π/2) of the vesicle gets fattened, and has lower curvature than the curvature at the side that faces the shear. This gradient in curvature induces the rotation of the leading edge cluster, as shown in Figure 3G.

The dynamics in our model is driven by minimization of energy and work. In Figure 3H we plot the time evolution of the different energy components: the bending energy W _b , the binding energy W _d , and the adhesion energy W _ad. At the steady-state configurations, when the vesicle is polarized against the flow, we find that the bending energy is increased due to shear, compared to the case when shear is absent. So clearly this energy is not minimized during the reorientation process. The protein-protein binding energy at short times is largest when the leading edge cluster is destabilized by the parallel shear flow (green line), but at long times this energy is essentially unaffected by the presence of shear.

Finally, we find that the adhesion energy W _ad is clearly smaller (more negative) at long times in the presence of shear compared to the no-shear case, as indicated by the two-headed orange arrow. We plot the adhesive area of the vesicle that is in contact with the substrate, for the case with and without shear, in the inset of the rightmost panel of Figure 3H, showing the bigger adhesive area when the shear is present. At a lower shear flow parameter, we found that the results are qualitatively the same, but the reorientation dynamics take longer to occur. On the other hand, large shear can cause a complete reversal of the polarity of the motile vesicle, without the U-turn trajectory, as we show in Supplementary Figure S5 (Supplementary Video S5). Such shear-induced polarity reversal is very similar to the dynamics observed in Dictyostelium discoideum cells in shear flow (Dalous et al., 2008). Curved membrane proteins of lower spontaneous curvature ( C 0 = 0.8 l min − 1 ) give rise to a less stable leading-edge cluster, which therefore causes breakup and reversal of the motile vesicle at lower shear flow (Supplementary Figure S6). Note that at even lower spontaneous curvature we lose the motile phenotype (Sadhu et al., 2021).

In Supplementary Figure S2B–D we show that under conditions where the flow is dominated by large slip-length, such that effectively a = 0, b ≠ 0 (Eq. (6)), the motile vesicle does not perform a U-turn to face the flow. The reason is that due to the constant flow profile there is much weaker deformation of the leading edge of the motile vesicle, compared to the case of linear flow gradient a ≠ 0, b = 0, and the vesicle therefore maintains its persistent direction of migration.

We therefore conclude that our simplified, minimal-cell model can provide a physical mechanism for the stabilization of cell migration that is upstream in the presence of shear flow. The basis for this mechanism is the increased cell spreading due to shear flow (Figure 3H), which was also observed in cells (Chotard-Ghodsnia et al., 2007; Dominguez et al., 2015). Note that cells respond to shear also through signalling that modify the overall cell behavior, which corresponds in our model to changes to the model parameters. Nevertheless, the physical mechanism that we find for migration against the flow is not cell-type-specific and is independent of any complex biochemical signalling. It may therefore explain why this behavior appears in many different cell types, as listed in the introduction.

Next, we explore the response to shear flow for non-motile adhered vesicles in our model. The non-polar, non-motile phenotypes of adhered vesicles in our model have several shapes (Sadhu et al., 2021): At low adhesion or low active force, the vesicle is weakly spread and has a roughly hemispherical shape (parameters indicated by the pink star in Figure 1). For high concentration of the curved membrane proteins, and at sufficiently high adhesion or active force, the vesicle spreads into a round pancake-like shape with a closed, circular leading-edge. The response of both of these shapes to the shear is given in the Supplementary Material (Supplementary Figure S3, S4), where its shown that they roll and slide with the flow.

A more interesting shape arises on surfaces with weaker adhesion, or at high active forces, where the vesicle spreads into a two-arc shape (Figure 4A, parameters indicated by the yellow star in Figure 1). The vesicle is elongated by two leading-edge clusters at opposing ends, with the membrane between them being pulled into a cylindrical shape. Adhered cells often have such elongated shapes, with multiple, competing leading edge lamellipodia, which render them non-polar and non-motile (Pankov et al., 2005; Schaufler et al., 2016; Singh et al., 2020; Dimchev et al., 2021).

In Figure 4B we plot the trajectories of the two-arc vesicle for three different shear conditions with respect to the initial long axis of the vesicle (the long axis is calculated as explained in the SI section S3). When shear is absent, the vesicle is almost completely non-motile, while in the presence of flow the vesicle moves with the shear flow. Interestingly, we can see that initially the vesicle moves faster when the shear is in the same direction as its body axis, compared to when the shear is perpendicular to the vesicle’s long axis. However, this migration along the long axis is unstable, and at long times the vesicle rotates to being perpendicular to the flow (Figure 4C), which is the stable configuration. In Figure 4D we show snapshots of the vesicle at different times as it moves with the shear flow, either parallel or perpendicular to the initial long axis of the vesicle.

As shown in Figure 4B, the vesicle is moving faster when the shear is along the body axis compared to when the shear is perpendicular to the body axis. This indicates that the shear is inducing some polarization of the active forces, which now have a net force that contributes to the active motility along the shear flow. To understand the origin of this shear-induced polarization, we analyzed the two leading-edge clusters at the opposing ends of the vesicle. We denote the arc pulling towards positive x direction (with the flow) and negative x direction (against the flow), arc-1 (red) and arc-2 (blue) respectively (Figures 4E,F). In Figure 4E we show that the net active force due to the two clusters is positive at the early times, indicating that indeed there is a net active force from the leading edges, pulling the vesicle with the flow direction. This shear-induced asymmetry is manifested as a larger active force along x-direction due to arc-1 compared to the negative component from arc-2. However, the sizes (number of active proteins) of the two arcs N _cluster do not show any systematic difference between the two leading edge clusters. Nevertheless, the net force in the flow direction due to arc-1 is stronger than the force due to arc-2 since the efficiency, defined as the net force along the flow direction per protein | F x act | / N cluster , of the proteins in arc-1 is larger (Figure 4E). Even when the size of arc-1 is smaller than arc-2, the proteins in arc-1 can be more efficient and produce a stronger net force along the flow direction, inducing motility of the vesicle in the presence of shear.

To get more insight into this efficiency of the two leading edges, we plotted the positions of the proteins on the x-y plane as shown in Figure 4F, at a time where the sizes of the two clusters is almost identical, yet there is a net force in the direction of arc-1. We fit a circular arc and find the radius of curvature for each leading edge cluster using the gradient-descent method. We find that the radius of curvature is bigger for arc-1 (R = 10.2l _min) compared to arc-2 (R = 7.4l _min), and this flatter shape of arc-1 makes its proteins’ active forces more oriented along the flow, compared to the orientations of the active forces in arc-2. The flatter shape of the leading edge of arc-1 is due to the shear forces pushing membrane along the tubular part that connects the two leading edges (see Supplementary Figure S7), such that membrane area is forced to flow from the region of arc-2 to that of arc-1, allowing the fan-shaped region of arc-1 to grow larger in area.

A very recent experimental study (Park et al., 2023), on Microglia (a type of glial cells), found shape and migration dynamics in the presence of shear that remarkably resemble the dynamics that we obtain in our model for the two-arc vesicle. In these experiments it is observed that the shear flow induces an increase in the size of the lamellipodial leading-edge in the direction of the flow, thereby breaking the symmetry of the cell and inducing its migration with the flow, in a manner that is identical to our model’s prediction (Figure 4F).

In the stable phase, where the vesicle moves with the shear flow that is perpendicular to its body axis, We plotted the cross-sectional area of the vesicle along its middle point (x _avg) at different times (Figure 4G). By following one particular node we illustrate that the membrane is rolling on the surface due to the flow-induced shear forces. In the Supplementary Material we present the dynamics of the two-arc shape at different shear flow strengths.

Weakly polarized cells, such as Chinese hamster ovary (CHO) cells, exhibit weak migration with the shear flow, often maintaining an elogated shape that is perpendicular to the shear direction (Dikeman et al., 2008). These CHO cells tend to spread in a circularly symmetric manner, and indeed our vesicles that spread uniformly tend to slide with the shear flow (Supplementary Material; Supplementary Figure S4), as observed for these cells.

Our “minimal-cell” model, where cell spreading and migration emerges due to curved membrane proteins that recruit the protrusive forces of actin polymerization, is used to explore the effects of shear forces applied to the membrane due to an imposed fluid flow. This model shows that since the self-organization of the curved proteins and active forces are dependent on the membrane shape, the system is strongly affected by these flow-induced shear forces.

We found that the motile crescent-shaped vesicle in our model spontaneously migrates against the shear flow due to the reorganization of curved protein in response to the shear flow. This behavior arises simply from the physics of minimizing the adhesion energy to the substrate. Since our mechanism is based on a very simple model, and a physical mechanism, it may explain why the tendency of cells to migrate against the flow appears in many different cell types that migrate using lamellipodia protrusions (Décavé et al., 2003; Fache et al., 2005; Dalous et al., 2008; Steiner et al., 2010; Ostrowski et al., 2014; Dominguez et al., 2015; Follain et al., 2018; Anderson et al., 2019; Alghanem et al., 2021). Though our model does not include many cellular components, it offers an explanation for the origin of this prevalent migration response to shear flow, which is not understood at present.

For the non-motile vesicles we found that they tend to migrate or roll with the shear flow, which may explain why weakly motile cells tend to move with the flow (Décavé et al., 2003; Fache et al., 2005; Dikeman et al., 2008). Note that cells respond to shear flow due to signalling pathways (Rose et al., 2007; Chistiakov et al., 2017), which modify the overall cell-substrate adhesion and cytoskeleton activity. This layer of biochemical control manifests as modifications to the parameters of the vesicle in our model, beyond the shear-induced shape changes that we investigated.

Our results may explain the observations regarding the adhesion-dependence of the direction of T Cell migration under shear flow (Dominguez et al., 2015): The migration upstream on ICAM-I coated surfaces corresponds to the behavior we expect at high surface adhesion for polarized cells. Indeed T Cells on VCAM-1 are more motile compared to VCAM-1 coated surfaces (Dominguez et al., 2015). The less polar and less motile T Cells on VCAM-1 coated surface (Dominguez et al., 2015) correspond to lower adhesion strength in our model (compare points indicated by the yellow and blue stars in Figure 1), which explains why they move with the fluid flow. Mixing the two adhesion molecules allows the stronger adhesion coating (ICAM-1) to dominate, which explains why the combination of ICAM-1 and VCAM-1 leads to upstream migration.

While it is satisfying that our simple model explains many qualitaive features of the responses of cells to shear flow, we need to remember that in a very complex system such as a living cell different mechanisms may end up producing similar-looking behavior. Furthermore, the physical mechanism that emerges in our model for the stabilization of the upstream migration, namely, the increased flattening of the cell and the resultant increase in cell-substrate adhesion, could possibly arise in other models that do not contain the basic ingredients of our model.

The MC model we used here does not describe the full fluid-flow field surrounding the cell, it only qualitatively captures the main properties of the flow-induced shear forces. Note that the lift-pressure that forms at the stagnation point of the flow near the cell-substrate contact line (Cantat and Misbah, 1999), which is currently absent from our calculation, may play an important role in “peeling” the cell from the substrate and will have to be added to the model. Future studies that include explicitly the fluid dynamics (Noguchi and Gompper, 2004; Mauer et al., 2018) and additional cellular components (Dabagh et al., 2017), may be used to explore the dynamics predicted by our model with better physical realism, at the price of greatly increased complexity and computation time. More complex flows, in channels and in the presence of complex geometric constraints, will be challenging to implement in this model, and can be explored in the future.

Our results may also explain the motility response of cells to external forces that are exerted on them by other means, not due to fluid flow. In (Weber et al., 2012) it was shown that when an adhered cell is pulled by a magnetic beads that is attached to the cell, it tends to polarize in the opposite direction to the applied force. Similarly, when cells that are attached to each other exert a pulling force on each other, they tend to polarize in opposite directions to each other (Weber et al., 2012). This plays an important role during collective cell migration (Mayor and Etienne-Manneville, 2016). These responses may arise from the same behavior that we obtain here, namely, the pulling force tends to stabilize the leading edge of the cell in the direction that is opposite to the direction of the external force (Figure 2).