## 1 Introduction

Cells change their shape and crawl in response to external stimuli by extruding and remodelling cell membrane protrusions, typically in the form of filopodia (finger-like protrusions) or lamellipodia (sheet-like protrusions). These cell shape changes are mainly due to the reorganisation of the filamentous actin (F-actin), a complex and highly dynamic composite contained in the outer rim of eukaryotic cells. Indeed, these filaments are formed from monomeric G-actin and exhibit a fast growing end (barbed end) and a slow growing end (pointed end), according to the ratio of the association rate to the dissociation rate of globular G-actin monomers. Actin polymers have a very high gelation efficiency and form viscoelastic gels at volume fraction below 0.1, thus controlling the viscoelasticity of the cytoplasm. Changes in F-actin levels and F-actin organisation are then intimately associated with the cell biological function and its specific morphologies.

However, the mechanisms driving the formation of cellular protrusions are still controversial. Cell protrusions can be generated by the formation of F-actin bundles which push forward the cytoplasmic membrane in specific places. An alternative view considers that the initiating process is the local polymerisation of an actin gel in a free volume located between the plasma membrane and the cell actin cortex 〚1〛.

While it is usually considered that cell membrane fluctuates without any privileged spatial direction in space nor particular time coherence, it is thus of particular interest to study cell protrusive activity as a self-organised dynamical process. Indeed, different results support the self-organised character of these cell shape changes. Ehrengruber et al. 〚2〛 show that neutrophils undergo periodic cytoskeletal rearrangements that lead to cycles of shape changes, ultimately resulting in cell translocation. In neutrophils, sinusoidal oscillations of F-actin, with periods of 8–10 s, were reported 〚3〛. Other studies also reported other types of self-organised behaviour, like the existence of rotating waves around the cell body 〚4, 5〛. This suggests isovolumetric extension and retraction, i.e. cyclic alteration of cell body size related to veil-like lamellipods changes.

Different models have been proposed to explain the periodic oscillatory movement of the cell membrane and the role of actin polymerisation/depolymerisation cycles in the recurring lamellipod protrusion and retraction dynamic. They usually consider the mechanical interactions of the actin cytoskeleton with: (i) the cytoplasm, through viscoelastic stresses generated by the actomyosin complex; (ii) the cell membrane, through tension forces due to the local membrane curvature. Oster and Perelson 〚6〛 proposed a model of lamellipodial motion based on the physical chemistry of actomyosin gels. The model consists of a sheet of cytogel attached to the substratum by elastic tethers. Permeation at the leading membrane raises the internal calcium concentration, activating the solation factors and the actomyosin contractile machinery. Peskin et al. 〚7〛 formulated a theory to account for the force generated by the polymerisation process itself when the filaments are rigid. They proposed that the addition of sub-units to the end of growing filaments rectified the Brownian motion of any diffusing object in the front of the filament, and showed that this ratcheting of diffusive motion could generate sufficient force to account for a number of motile phenomena. Molginer and Oster 〚8〛 generalise this model when the thermal fluctuations of either the actin fibres or membrane are significant. This extension provides a mechanical explanation for the protrusion of lamellipodia. It is indeed known that temperature can induced sol–gel transition in cross-linked actin networks.

Other models have been proposed, in which biochemical processes are described through a reaction-diffusion model which takes into account the interactions of two states of F-actin, free and linked to nucleation proteins. Oscillatory actin waves are generated within a fixed area modelling the cell periphery 〚9〛. Conversely, in most of the mechanical models, the periodic oscillating behaviour of the cell membrane protrusion is generated by the coupling between intracellular pressure (hydrostatic...) and F-actin network contractility, as in the cytomechanical model of Alt and Tranquillo 〚5〛.

In this paper, we consider this latter model as a theoretical basis for further theoretical investigations of cell–cell interactions. In particular, we propose different extensions of the original model which take into account cell–cell interactions and the influence of extracellular factors through modulation of the membrane protrusivity and cell membrane elasticity. We further analyse how the coupling of the self-oscillating state of the cell with an extracellular gradient could initiate cell migration. As an experimental counterpart of this cell simulated behaviour, we consider fibroblasts morphological changes as they can be observed on videomiscroscopy sequences. Starting from a previous characterisation of the spontaneous pulsating behaviour of isolated fibroblast 〚10〛, we discussed to which extent modifications of the experimentally observed fibroblast shapes can be interpreted as perturbations or bifurcations of the spontaneous self-oscillating cell protrusive behaviour. In a second part of the paper, we further analyse the response of this cytomechanical model to an extracellular gradient of chemoattractant and show that subsequent cell migration can be initiated by only considering local modifications of the cell cortex–membrane properties.

## 2 Experimental observations of autonomous cell pulsating behaviour

Characterisation of cell membrane deformation dynamics emerges from the contour maps which gives a spatio-temporal picture of the cell morphology. These maps correspond to the reconstruction of the cell boundary (or contour) in polar coordinates where the membrane radial exension L_{exp} measured from the cell estimated barycentre is plotted against the angular position θ at each time, for distance values greater than ${\overline{L}}_{\mathrm{exp}},$ where ${\overline{L}}_{\mathrm{exp}}$ is the mean length of the membrane extension. This allows to exhibit the protrusive zones of the cell (Fig. 1c).

Observations of sub-confluent L929 fibroblasts cultured in vitro on two-dimensional rigid substratum showed that isolated cells spontaneously developed cyclic and pulsating movements of extension/retraction as exemplified in Fig. 1a and 1b, where co-ordinated periodic changes of cell morphology develop between two quasi-perpendicular directions of protrusion. The four protrusive directions involved are clearly exhibited by the associated contour map in Fig. 1c which also reveals their persistence with time.

When cells are no longer spatially isolated, this spontaneous cyclic behaviour interacts with the pulsating protrusive activity of neighbouring cells. We then observe modifications or disappearance of the directions of protrusions, as if cells try to avoid interactions with the neighbouring cells through some kind of physical contact–inhibition mechanism 〚11〛.

In order to characterise interacting cells protrusive behaviour, we define a so called influence domain for each cell within a given group of cells (Fig. 2). This domain is schematically represented by a circle centred on the cell nucleus barycentre, with a radius defined by the maximum recorded amplitude of the cell protrusions in the video-microscopy sequence during the considered observation duration. The intersection of two (or more) influence domains thus represents potential interaction zones, which can be superimposed to the time-varying contour maps determined for each cell.

This characterisation provides a global picture of the contact inhibition phenomenon, which has been reported for various types of cells. We observed that cell protrusions preferentially appear outside the darkest interaction zones (Fig. 2), supporting the idea that contact inhibition is triggered in most cases by physical contacts between cells, which then retract and try to develop their protrusions in other empty spaces. However, since cell protrusive activity is a spatio-temporal coordinated process, it remains to precise how simultaneous competing lamellipodia can integrate such external constraints and trigger morphological cell bifurcation between a free running shape to a constrained one.

## 3 A cytomechanical model for cell deformations

### 3.1 The model

We used as a modelling basis the cytomechanical model proposed by Alt and Tranquillo 〚5〛 to describe the spontaneous dynamics of cell deformations, with specific application to the analysis of keratinocytes spatio-temporal deformations. Indeed, the model involves a minimum number of variables, while incorporating salient features of biomechanical cell deformations.

The model considers that cell protrusions dynamics are due to the biophysical properties (visco-elasticity, contractility,...) of the cortical polymer layer or network of actin and myosin filaments which underlies the cell membrane and surrounds the cell body and nucleus. Moreover, this more or less dense network is able (i) to disassemble at locations where it becomes condensed, (ii) to reassemble in cell protrusions, like lamellipodia. Cell protrusions then mainly results from a mechanical non-equilibrium state of the pressure and forces acting on the cell periphery. Subsequent F-actin polymerisation/depolymerisation in a moving protrusion is a consequence of the so-induced intracellular space variation. The networks in the cortex and lamellipodia are assumed to have the same biophysical properties, except the boundary condition imposed by the additional cortex–membrane surface tension. The first part of the model depicts the dynamics of the actin in the cell cortex which is mainly responsible for cell shape reorganisation through its ability (i) to polymerise into F-actin and depolymerise into G-actin, (ii) to interact with myosin to generate the contractile activity in the cell, (iii) to move through the cytoplasm via convection. The local amount of F-actin also determines the intensity of the resistive stress applied on the membrane as the result of actin cytoskeleton-membrane proteins crosslinks. This resistive stress, together with the cortex–membrane curvature-induced stress, balances the intracellular hydrostatic pressure which acts as the ‘driving force’ in the model.

In order to transform the two-dimensional original problem into a simpler one-dimensional version, we moreover assume as in Alt and Tranquillo 〚5〛 that the F-actin density as well as its convective tangential velocity is constant in the radial direction. The first assumption is justified when the width of the cell cortex $L\left(\theta ,\text{\u2006}t\right),$ remains small in front of the radius of the cell body R_{0}. The second implies that the cortical network can slip with respect to the membrane, which is also partly justified if considering that membrane proteins involved in the network-membrane connections are mobile within the membrane.

The three variables considered in the cytomechanical model are thus: (i) the F-actin concentration in the cortex $a\left(\theta ,\text{\u2006}t\right),$ (ii) the F-actin tangential velocity $v\text{\u2006}\left(\theta ,\text{\u2006}t\right)$ (taken at the periphery of the cell body delimited by a unit circle of radius R_{0}=1), (iii) the width of the cell cortex annulus $L\text{\u2006}\left(\theta ,\text{\u2006}t\right)$ – or the membrane position taken from the surface of the cell. The spatio-temporal evolution of these three variables is given by the following set of partial differential equations.

The first one stands for the conservation of F-actin in the cortex. Actin polymerisation and depolymerisation are considered to occur with the same rate η. Moreover, the state of polymerisation or depolymerisation depends on the local value of F-actin concentration relatively to the stationary concentration a_{*}.

$$\frac{\partial}{\partial \mathrm{t}}\text{\u2006}\left(L\text{\u2006}a\right)+\frac{\partial}{\partial \theta}\text{\u2006}\left(L\text{\u2006}a\text{\u2006}v\right)=\eta \text{\u2006}L\text{\u2006}\left({a}_{*}-\mathrm{a}\right)$$ |

The second equation describes the balance of forces applied on the cortex–membrane in the radial direction. The model takes into account a viscous friction that represents the cell adhesion on the substratum, modulated by the coefficient Φ_{1}, the intracellular hydrostatic pressure ${\u03d0}_{1},$ the resistive elastic stress due to the cortical network membrane attachment with the elasticity coefficient γ_{1} and a curvature-dependent stress due to the surface tension of the cortex–membrane complex, characterised by the coefficient τ_{1}.

$${\Phi}_{1}\text{\u2006}a\text{\u2006}\frac{\partial \mathrm{L}}{\partial \mathrm{t}}={\u03d0}_{1}-{\gamma}_{1}\text{\u2006}L\text{\u2006}\mathrm{a}+\frac{\partial}{\partial \theta}\text{\u2006}\left({\tau}_{1}\text{\u2006}a\text{\u2006}\frac{\partial \mathrm{L}}{\partial \theta}\right)$$ |

The third equation concerns the balance of forces in the tangential direction. It includes the frictional force on the network moving in the cytosol, with magnitude controlled by the drag coefficient Φ_{0} (related to the attachment of the network with the membrane and the other types of filaments), a viscous stress with viscosity coefficient μ_{0} and the membrane curvature induced stress with coefficient τ_{0}.

$${\Phi}_{0}\text{\u2006}a\text{\u2006}\mathrm{v}=\frac{\partial}{\partial \theta}\text{\u2006}\left[{\mu}_{0}\text{\u2006}a\text{\u2006}\frac{\partial \mathrm{v}}{\partial \theta}+{\sigma}_{0}\left(a\right)-\frac{\partial}{\partial \theta}\text{\u2006}\left({\tau}_{0}\text{\u2006}a\text{\u2006}\frac{\partial \mathrm{L}}{\partial \theta}\right)\right]$$ |

The contractile activity of the actomyosin network is modelled by the nonlinear function ${\sigma}_{0}\left(a\right).$ Two mechanical states can be distinguished according to the value of the network F-actin concentration $a\left(\theta ,\text{\u2006}t\right)$ with respect to a saturation concentration a_{sat}. If $a\left(\theta ,\text{\u2006}t\right)$ is below the saturation concentration value, the contractile stress increases according to a parabolic law. Above the saturation threshold, the stress decreases exponentially as a consequence of the network swelling. The nonlinear function proposed by Alt and Tranquillo (1995) 〚5〛 is the following:

$${\sigma}_{0}\left(a\right)={\psi}_{0}\text{\u2006}{a}^{2}\text{\u2006}{e}^{-\mathrm{a}/{\mathrm{a}}_{\mathrm{sat}}}$$ |

_{0}controls the magnitude of the contractile stress.

The partial differential equations system is nondimensionalised by setting the following dimensionless variables:

$$\begin{array}{ccc}\multicolumn{1}{c}{\tilde{t}=\mathrm{t}\text{\u2006}\eta}& \multicolumn{1}{c}{\tilde{a}=\frac{a}{{a}_{0}}}& \multicolumn{1}{c}{\tilde{L}=\frac{L}{{R}_{0}}}\end{array}$$ |

The normalised parameters are then:

$$\begin{array}{ccc}\multicolumn{1}{c}{\tilde{\u03d0}=\frac{{\u03d0}_{1}}{{\Phi}_{1}\text{\u2006}{a}_{0}\text{\u2006}{R}_{0}\text{\u2006}\eta}}& \multicolumn{1}{c}{\tilde{\gamma}=\frac{{\gamma}_{1}}{{\Phi}_{1}\text{\u2006}\eta}}& \multicolumn{1}{c}{\tilde{\mu}=\frac{{\mu}_{0}}{{\Phi}_{0}}}\end{array}$$ |

$$\begin{array}{cccc}\multicolumn{1}{c}{\tilde{\psi}=\frac{{\psi}_{0}\text{\u2006}{a}_{0}}{{\Phi}_{0}\text{\u2006}{R}_{0}\text{\u2006}\eta}}& \multicolumn{1}{c}{{\tilde{a}}_{\mathrm{sat}}=\frac{{a}_{\mathrm{sat}}}{{a}_{0}}}& \multicolumn{1}{c}{\tilde{\tau}=\frac{{\tau}_{0}}{{\Phi}_{0}\text{\u2006}{\eta}^{2}}}& \multicolumn{1}{c}{\tilde{d}=\frac{{\tau}_{1}\text{\u2006}{\Phi}_{0}}{{\tau}_{0}\text{\u2006}{\Phi}_{1}}}\end{array}$$ |

Dropping the tilde for notational simplicity, we obtain the following system of dimensionless differential equations:

$$\frac{\partial}{\partial \mathrm{t}}\text{\u2006}\left(L\text{\u2006}a\right)=-\frac{\partial}{\partial \theta}\text{\u2006}\left(L\text{\u2006}a\text{\u2006}v\right)+\mathrm{L}\text{\u2006}\left(1-\mathrm{a}\right)$$ |

$$a\text{\u2006}\frac{\partial \mathrm{L}}{\partial \mathrm{t}}=\u03d0-\gamma \text{\u2006}L\text{\u2006}\mathrm{a}+\mathrm{d}\text{\u2006}\frac{\partial}{\partial \theta}\text{\u2006}\left(\tau \text{\u2006}a\text{\u2006}\frac{\partial \mathrm{L}}{\partial \theta}\right)$$ |

$$a\text{\u2006}\mathrm{v}=\frac{\partial}{\partial \theta}\text{\u2006}\text{\u2006}\left[\mu \text{\u2006}a\text{\u2006}\frac{\partial \mathrm{v}}{\partial \theta}+\sigma \text{\u2006}\left(a\right)-\text{\u2006}\frac{\partial}{\partial \theta}\text{\u2006}\left(\tau \text{\u2006}a\text{\u2006}\frac{\partial \mathrm{L}}{\partial \theta}\right)\right]$$ |

The linear stability analysis of this system of equations around the homogeneous steady state ${L}_{\mathrm{s}}=\u03d0/\gamma \text{\u2006},$ a_{s}=1 and v_{s}=0 is performed in a standard way by analysing perturbations of the form $\mathrm{exp}\left(\lambda \text{\u2006}\mathrm{t}+\mathrm{i}\text{\u2006}m\text{\u2006}\theta \right).$ The evolution with time of perturbations with spatial mode m is given by the roots of the associated dispersion equation, i.e. the pair of eigenvalues

$$\lambda \text{\u2006}\left(m\right)=\frac{1}{2}\text{\u2006}\left[-1+\frac{{m}^{2}\text{\u2006}\rho}{1+\mu \text{\u2006}{m}^{2}}-{\mathrm{m}}^{2}\text{\u2006}d\text{\u2006}\tau \pm \text{\u2006}\sqrt{{\left(2\text{\u2006}\gamma -1+\frac{{m}^{2}\text{\u2006}\rho}{1+\mu \text{\u2006}{m}^{2}}+{\mathrm{m}}^{2}\text{\u2006}d\text{\u2006}\tau \right)}^{2}-4\text{\u2006}\gamma \text{\u2006}\left(\gamma +{\mathrm{m}}^{2}\text{\u2006}d\text{\u2006}\tau \right)-4\text{\u2006}\frac{{m}^{4}\text{\u2006}\u03d0\text{\u2006}\tau}{1+\mu \text{\u2006}{m}^{2}}}\right]$$ |

Looking for oscillatory solutions of the model leads to the analysis of Hopf bifurcations and complex conjugated eigenvalues. The real part of the eigenvalues is independent of the coefficient $\u03d0$ and γ. The key parameters for unstable modes determination are thus the viscosity coefficient μ, the surface-tension coefficient τ and the drag coefficient d. As it can be intuitively understood, high values of μ and τ tends to strengthen the cell cortex–membrane complex and thus tends to reduce the number of unstable modes. At the opposite, high value of σ, which represents the variation amplitude of the actomyosin network contractility, can favour the destabilisation of a larger number of spatial modes.

### 3.2 How to take into account cell–cell interactions and cell migration with the cytomechanical model

The model has been successfully used to describe the spontaneous dynamics of various cell lines, including leucocytes 〚12〛 and keratinocytes 〚13〛. Our aim is now to investigate if such a model is a valuable representation for analysing the various aspects of cell interactions through a continuous and integrated description of the cell dynamics. We will especially focus on cell–cell interactions as well as cell migration induced by a gradient of extracellular factor. Indeed, in both situations cell detects and integrates external stimuli, and exhibits a variety of responses, such as contact inhibition and/or the appearence of new dynamical states like migration. Complex pathways of signal transduction are involved in the elaboration of the cell response and these processes are clearly not explicitely described in our present cytomechanical model. However, we can at least simulate the corresponding effects of this signalisation, as reported experimentally, and then examine their consequences on the cell simulated dynamics. Then, adequation of the model simulations and observed cell dynamics response will provide an additional step for the model validation. For example, the observed polarisation of a cell submitted to a gradient of extracellular factor suggests that the membrane morphological instability leading to the polarisation of the membrane is the result of the decrease of the membrane tension 〚14〛. In our model, this effect can be modelled by considering a dependence of the tension coefficient with regard to the local concentration of extracellular factor on the membrane. One further can assume that the inclusion of molecules in the membrane leads to a decrease of the tension, as suggested by various studies on membrane properties 〚15, 16〛.

## 4 Modelling cell–cell interactions as perturbations of the spontaneous oscillatory protrusions dynamics

### 4.1 The model

Following the experimental results briefly reported in the second paragraph, we consider the cytomechanical model presented in the previous section as a basis for theoretically analysing physical cell–cell contact inhibitions as possible perturbations of the cell spontaneous periodic protrusive behaviours. The basic assumption is that such a mechanical external signal can be transduced, through different signalling pathways, into structural modifications of the cell cortex, thus resulting in a globally re-adapted co-ordination of cell protrusions. At the cell level, such response would lead to modifications of cell morphologies and/or changes of the characteristic parameters of the cell protrusive behaviour (periodicity, protrusion amplitude...).

We investigated here two different mechanical responses. The first one assumes a local decrease of the protrusive force (coefficient $\u03d0$) in the cell–cell interaction area. The second one considers that stress receptor activation can lead to a local re-organisation of the cortical network–membrane attachment, leading to a decrease of the curvature-induced surface tension. We thus consider the following step functions to describe the angular variation of the respective model parameters in any given interaction window $\left[{\theta}_{1},\text{\u2006}{\theta}_{2}\right]:$

$$\u03d0\left(\theta \right)=\{\begin{array}{ccc}\multicolumn{1}{c}{\u03d0}& \multicolumn{1}{c}{\mathrm{if}}& \multicolumn{1}{c}{\theta <{\theta}_{1}\mathrm{and}\text{\u2006}\theta >{\theta}_{2}}\\ \multicolumn{1}{c}{{\u03d0}_{\mathrm{int}}}& \multicolumn{1}{c}{\mathrm{if}}& \multicolumn{1}{c}{{\theta}_{1}\le \theta \le {\theta}_{2}}\end{array}\tau \left(\theta \right)=\{\begin{array}{ccc}\multicolumn{1}{c}{\tau}& \multicolumn{1}{c}{\mathrm{if}}& \multicolumn{1}{c}{\theta <{\theta}_{1}\mathrm{and}\text{\u2006}\theta >{\theta}_{2}}\\ \multicolumn{1}{c}{{\tau}_{\mathrm{int}}}& \multicolumn{1}{c}{\mathrm{if}}& \multicolumn{1}{c}{{\theta}_{1}\le \theta \le {\theta}_{2}}\end{array}$$ |

The model thus becomes:

$$\frac{\partial}{\partial \mathrm{t}}\text{\u2006}\left(L\text{\u2006}a\right)=\text{\u2006}-\text{\u2006}\frac{\partial}{\partial \theta}\text{\u2006}\left(L\text{\u2006}a\text{\u2006}v\right)+\mathrm{L}\text{\u2006}\left(1-\mathrm{a}\right)$$ |

$$a\text{\u2006}\frac{\partial \mathrm{L}}{\partial \mathrm{t}}=\u03d0\text{\u2006}\left(\theta \right)-\gamma \text{\u2006}L\text{\u2006}\mathrm{a}+\mathrm{d}\text{\u2006}\frac{\partial}{\partial \theta}\text{\u2006}\left(\tau \text{\u2006}\left(\theta \right)a\text{\u2006}\frac{\partial \mathrm{L}}{\partial \theta}\right)$$ |

$$a\text{\u2006}\mathrm{v}=\frac{\partial}{\partial \theta}\text{\u2006}\left[\mu \text{\u2006}a\text{\u2006}\frac{\partial \mathrm{v}}{\partial \theta}+\sigma \text{\u2006}\left(a\right)-\frac{\partial}{\partial \theta}\text{\u2006}\left(\tau \text{\u2006}\left(\theta \right)\text{\u2006}a\text{\u2006}\frac{\partial \mathrm{L}}{\partial \theta}\right)\right]$$ |

$$\frac{\partial \tau \text{\u2006}\left(\theta \right)}{\partial \theta}=\left(\tau -{\tau}_{\mathrm{int}}\right)\text{\u2006}\left[\delta \text{\u2006}\left(\theta -{\theta}_{2}\right)-\delta \text{\u2006}\left(\theta -{\theta}_{1}\right)\right]$$ |

### 4.2 Simulation of mechanically perturbed cell oscillatory morphologies

As initial cell morphology and dynamical state, we consider the spontaneous pulsating movements of extension/retraction of the cell membrane between two perpendicular directions observed for isolated L929 fibroblasts 〚10〛. This behaviour has been simulated with the following set of model parameters:$\u03d0=\gamma =0.5,$ μ=2, τ=1.5, δ=0.001 a_{sat}=2 and ψ=5.5. The corresponding spatio-temporal map $L\text{\u2006}\left(\theta ,\text{\u2006}t\right)$ of the simulated membrane deformations is plotted in Fig. 3a. This asymptotic periodic steady-state is taken as the initial condition for the simulation of cell–cell mechanical interactions. Above described modifications of the model parameters are then applied to simulate the cell response to the external mechanical stimulus and the progressively emerging dynamical state is recorded. Except in Fig. 3a, where there is no interaction, the interaction window $\left[{\theta}_{1},\text{\u2006}{\theta}_{2}\right]$ (75 degrees for each simulation) is indicated at the bottom of the map.

In the simulation presented in Fig. 3b, we observed the modification of the protrusive dynamics starting from the interaction zone (t=0 to t=12), then influencing the overall dynamic of protrusion (t=12 to t=30). In this case, both the hydrostatic pressure and the surface tension have been decreased, with:

$$\u03d0\left(\theta \right)=\{\begin{array}{ccc}\multicolumn{1}{c}{0.5}& \multicolumn{1}{c}{\mathrm{if}}& \multicolumn{1}{c}{\theta <135\text{\u2006}\mathrm{and}\text{\u2006}\theta >210}\\ \multicolumn{1}{c}{0.4}& \multicolumn{1}{c}{\mathrm{if}}& \multicolumn{1}{c}{135\le \theta \le 210}\end{array}\tau \left(\theta \right)=\{\begin{array}{ccc}\multicolumn{1}{c}{1.5}& \multicolumn{1}{c}{\mathrm{if}}& \multicolumn{1}{c}{\theta <135\text{\u2006}\mathrm{and}\text{\u2006}\theta >210}\\ \multicolumn{1}{c}{0.75}& \multicolumn{1}{c}{\mathrm{if}}& \multicolumn{1}{c}{135\le \theta \le 210}\end{array}$$ |

In the second case (Fig. 4), only the alteration of the protrusive pressure $\u03d0$ has been considered.

We are now studying the effect on the protrusion dynamics of the interaction zone orientation respectively to the cell position. The simulated behaviour reported in Fig. 4a corresponds to:

$$\u03d0\left(\theta \right)=\{\begin{array}{ccc}\multicolumn{1}{c}{0.5}& \multicolumn{1}{c}{\mathrm{if}}& \multicolumn{1}{c}{\theta <165\text{\u2006}\mathrm{et}\text{\u2006}\theta >240}\\ \multicolumn{1}{c}{0.4}& \multicolumn{1}{c}{\mathrm{if}}& \multicolumn{1}{c}{165\le \theta \le 240}\end{array}$$ |

The simulation in Fig. 4b is realised in a similar way, except that the interaction zone has been shifted counter clockwise with a 30° angle, i.e.:

$$\u03d0\left(\theta \right)=\{\begin{array}{ccc}\multicolumn{1}{c}{0.5}& \multicolumn{1}{c}{\mathrm{if}}& \multicolumn{1}{c}{\theta <135\text{\u2006}\mathrm{et}\text{\u2006}\theta >210}\\ \multicolumn{1}{c}{0.4}& \multicolumn{1}{c}{\mathrm{if}}& \multicolumn{1}{c}{135\le \theta \le 210}\end{array}$$ |

In the first case (Fig. 4a), the cell protrusion dynamics is modified mainly in the interaction zone but the cell keeps a pulsating state. However, when the spatial location of the external stimulus has been shifted (Fig. 4b), we observe a global change of the cell dynamical state with the appearance of a rotating wave of deformation around the cell body.

## 5 From cell deformations to cell chemotaxis

### 5.1 A model of cellular migration induced by an extracellular factor

Cell chemotaxis can be decomposed into three basic dynamical steps, namely chemo-attractant receptor binding, intracellular signal transduction and induction of biased translocation. Facing the explosion of molecular information, mainly provided through the development of specialised probes, theoretical models can help to identify basic understanding of coupled cell dynamical features including spontaneous oscillatory cell shape changes, random migration and directional response to extracellular gradients. In this section, we analysed the possible induction of cell migration as a phenotypic response to the extracellular stimulus provided by a chemoattractant factor. We assume that this diffusible factor is released from a ponctual source located in the neighbourhood of the cell.

We do not explicitly consider the series of complex biochemical and biophysical responses elicited after chemoattractant binding to specific surface receptors (internalisation, activated G-protein transduction cascade, cytoskeletal association...) On the basis of the here considered cytomechanical model, our aim is rather to correlate temporally and spatially lamellipodia formation with cortical F-actin network dynamical and physical properties, associated with modified cell cortex–membrane surface tension changes. We thus assume that the emerging effect of the signalling pathways triggered by chemoattractant molecule binding to the cell membrane is a modification of the cortex–membrane surface tension, which can be related to structural modification of the large protein complexes and cluster of chemotactic receptors associated with the cell membrane 〚17〛. More precisely, we consider that the associated parameter τ depends linearly on the local concentration C of the extracellular factor at the membrane:

$$\tau \text{\u2006}\left(C\right)=\tau -\Lambda \text{\u2006}C$$ |

Λ is a coefficient which characterises the sensitivity to the extracellular factor and thus provides a control parameter for different down-regulation scheme of receptors sensitivity (decrease in the number of receptors with increasing chemoattractant concentration, integration of signal through the coupling of different transducers 〚17〛). On the other hand, surface tension still remains proportional to the cortical network F-actin concentration.

The corresponding cytomechanical model is thus defined by the following new system of equations:

$$\frac{\partial \left(L\text{\u2006}a\right)}{\partial \mathrm{t}}=-\text{\u2006}\frac{\partial \left(L\text{\u2006}a\text{\u2006}v\right)}{\partial \theta}+\mathrm{L}\text{\u2006}\left(1-\mathrm{a}\right)$$ |

$$a\text{\u2006}\frac{\partial \mathrm{L}}{\partial \mathrm{t}}=\u03d0-\gamma \text{\u2006}L\text{\u2006}\mathrm{a}+\mathrm{d}\text{\u2006}\frac{\partial}{\partial \theta}\text{\u2006}\left(a\text{\u2006}\tau \left(C\right)\text{\u2006}\frac{\partial \mathrm{L}}{\partial \theta}\right)$$ |

$$a\text{\u2006}\mathrm{v}=\frac{\partial}{\partial \theta}\text{\u2006}\left[\mu \text{\u2006}a\text{\u2006}\frac{\partial \mathrm{v}}{\partial \theta}+\sigma \left(a\right)-\text{\u2006}\frac{\partial}{\partial \theta}\text{\u2006}\left(a\text{\u2006}\tau \left(C\right)\text{\u2006}\frac{\partial \mathrm{L}}{\partial \theta}\right)\right]$$ |

When the extracellular factor concentration C increases, the membrane stress decreases: this induces both a morphological instability and a cell polarisation toward the extracellular gradient (see next section).

However, it is clear that a temporary cell adhesion to the substratum is a necessary condition for the cell to exert the force required for forward cell body translocation. Cell adhesion is not described explicitly in the model but taken into account implicitly through the definition of a stress threshold above which the breakdown of cell–substrate links will occur. Cell–substratum adhesion/de-adhesion dynamics is thus computed as follows: cell polarisation leads to an asymmetry of the cell morphology relatively to the perpendicular axis of the translocation direction. At the same time, we assume that at the leading lamellipodia, the cell establishes new focal adhesions with the substratum. This morphological asymmetry is associated with the displacement of the cell geometric barycentre from its initial position. When this displacement exceeds a given threshold, interpreted as the admissible limit for the resulting stress applied on the cell cytoskeleton, the rear part of the cell is allowed to move forward up to the point where its centre of mass coincides with the cell geometric barycentre: then, reorganisation of the cell morphology occurs before the next translocation step. Therefore, as long as an extracellular gradient can be detected by the cell, the cell moves toward the source according to this three-phase mechanism, i.e. cell polarisation, displacement of the centre of mass, cell remodelling and reorganisation of intracellular stresses. Such modelling of cell crawling motion, which is characteristic of mesenchymal cells migration, summarises the five-step process described by Sheetz et al. 〚18〛.

### 5.2 Simulation of cell chemotaxis

We assume that the cell is initially non-stimulated, i.e. circular and non polar with a uniform cortical distribution of F-actin, together with a uniform distribution of receptors on the plasma membrane. Even if the chemoattractant is a diffusible factor, we assume for simplicity that the time scale of the diffusion process is large enough compared to the cell migration velocity, such that a quasi steady-state distribution $C\text{\u2006}\left(r\right)$ of the chemoattractant in the radial direction r can be considered. We take:

$$C\text{\u2006}\left(r\right)=\frac{Q}{2\text{\u2006}\sqrt{\pi \text{\u2006}D\text{\u2006}{t}_{0}}}\text{\u2006}\mathrm{exp}\text{\u2006}\left(-\frac{{r}^{2}}{4\text{\u2006}D\text{\u2006}{t}_{0}}\right)$$ |

_{0}of the 1D-diffusion equation

$$\frac{\partial \mathrm{C}}{\partial \mathrm{t}}=\mathrm{D}\text{\u2006}\frac{{\partial}^{2}C}{\partial {\mathrm{r}}^{2}}$$ |

Simulations have been realised with the following set of normalized model parameters: Q=45, D=4.0 and t_{0}=1 for the chemoattractant, $\u03d0=\gamma =0.5,$ μ=4, τ=1.5, δ=0.001, a_{max}=10, a_{sat}=2 and ψ=5.5 for the cell cytomechanical parameters. The initial cell state and morphology correspond to the spatially homogeneous steady state given by $\mathrm{L}=\u03d0/\gamma ,$ a=1 and v=0. The cell is initially located at a normalised distance d_{0}=4.5 from the punctual source of chemoattractant and the threshold for the displacement of the geometric barycentre is set to 0.1.

Three simulations have been performed, according to possible experimental situations. Fig. 5 shows the cell migration towards the chemoattractant source. The progressive destabilisation of the cell membrane, with concomitant appearance of leading protrusions and cell polarisation, are observed, in pace with a slow translocation in the direction of the source. When the cell reaches the source, it tends to collapse as the membrane still responds to the chemoattractant stimulus.

To avoid this effect, which is related to our steady-state approximation of the chemoattractant gradient, two new simulations have been performed. The first one (Fig. 6) simulates the displacement of the source. When the cell is about to reach the source, the source is moved instantaneously back to the initial position of the cell. The cell then reorganises its shape and polarises in the opposite direction, toward the new source location. The second case (Fig. 7) corresponds to an extinction of the source. Then, as expected, the cell stops moving and recovers its original homogeneous shape.

In order to simulate possible experimental control of cell migration through variations of the extracellular gradients, we simulated cell chemotaxis for different stiffness of the $C\left(r\right)$ spatial distribution. Modulation of the chemoattractant gradient is obtained by changing the parameters D or t_{0} in the analytical expression of $C\left(r\right).$ In the first simulation (Fig. 5), the migration velocity, evaluated for a cell with diameter of $15\text{\u2006}\mu \mathrm{m},$ is $0.65\text{\u2006}\mu \mathrm{m}$ per normalised-time unit when t_{0}=1. For increasing values of the parameter t_{0}, i.e for a smaller gradient of chemoattractant (Fig. 8), a nonlinear variation of the migration velocity is observed, initially increasing from 0.37 to $0.77\text{\u2006}\mu \mathrm{m}$ per normalised-time unit, then decreasing to $0.48\text{\u2006}\mu \mathrm{m}$ per time unit.

Another way of controlling the cell migration velocity is to change the admissible threshold of cell deformation, in relation with the stress threshold above which the link between the membrane receptors, like integrins, and the substratum are broken. The results are presented in Fig. 9, where the cell migration velocity has been computed for thresholds ranging from 0.06 to 0.12. As expected, stronger cell–substratum attachment is correlated with slower cell translocations. Globally, the cell velocity decreases in a nonlinear way from 0.67 to $0.56\text{\u2006}\mu \mathrm{m}$ per normalised-time unit with increasing values of the threshold.

## 6 Discussion

A large body of experiments has shown that cells change their shape in response to various environmental signals, let them be chemical or mechanical. Moreover, experimental data obtained with different cell-types have shown that cell shape changes can occur spontaneously, i.e. in absence of external stimuli. Thus, consideration of spontaneous cell shape changes appears as a pre-requisite for further investigations of directional cell migration induced by the dynamics of membrane protrusions. In this paper, we consider a ‘cortical layer’ model of the cell, i.e. an idealised circular cell with a contractive outer annulus mainly composed of F-actin. Dynamic changes of the cell radius result from a mechanical non-equilibrium state between an intracellular hydrostatic pressure, a curvature dependant surface tension stress, which is proportional to the network density of F-actin, and F-actin network contractility and swelling. This model, initially proposed by Alt and Tranquillo (1995) 〚5〛, is able to generate various patterns of spontaneous cycles of membrane extension and retraction, corresponding to oscillatory non-steady states that exist beyond some bifurcation values of the model parameters.

We successfully used this modelling approach to reproduce the self-sustained oscillatory protrusive activity experimentally observed with isolated L929 fibroblast cells 〚10〛. We here extend this model and further analyse cell–cell mechanical interactions and chemotactically induced cell migration by considering spatially inhomogeneous variation of the surface tension of the cortical actin/cell membrane complex. Our working hypothesis is that cell membrane protrusion is associated with a decrease in cell surface tension. This is in agreement with recent experimental data of Raucher and Sheetz (2000) 〚19〛 and with scanning acoustic microscopy measurements 〚20〛 which reveal that lamellipodia extensions take place in cellular region of relatively low stiffness. It has also been reported experimentally for phospholipids bilayers 〚16〛 that membrane molecules inclusion decreases the membrane stiffness.

Our model thus provides a continuous description of the successive stages of the cell morphological changes, from spontaneous oscillatory deformations up to polarisation and then migration. It furthermore provides realistic simulations of cell migration, which keeps some random aspects of the cell trajectory and in which the direction of locomotion is determined by the result of competing lamellipodia occurring simultaneously around the cell and cell progressive polarisation. We especially focus on chemotaxis, which is simulated by considering that the binding of chemoattractant molecules to membrane receptors triggers a transduction pathway which finally results in a decrease of the cortical surface tension. On the other hand, cell adhesion is necessary to insure further development of the cellular traction forces pulling the cell body; this is one of the five-step scenario described by Sheetz et al. (1999) 〚18〛 for cell migration: membrane extension, attachment to the substratum, cell contraction, release of the attachment at the trailing edge and recycling of the receptors. We do not explicitly consider all these steps but rather proposed a simplified view of this mechanism as an all or none process. If the polarisation of the cell exceeds some threshold corresponding to some value of intracellular stress, then local breakdown of cell–substrate will occur. Even if simple, this approach can still take into account data on cell adhesion forces as well as observations of integrins tracks left behind a moving cell as a testimony of the focal adhesions breakdown.

Another extension of this theoretical approach is the model ability to define optimum conditions for cell migration. Controlling cell migration is indeed crucial in many events such as wound healing, where fibroblasts or myofibroblasts migration is critical, or as tumour invasion, where a decrease of the migration speed would help to reduce cancer cell infiltration in normal tissue.

This article illustrates the influence of first the slope of the chemoattractant gradient, second the strength of the cell–substratum attachment. Another way to control cell migration is a direct modification of the cell characteristic parameters, including the cytomechanical parameters. One can thus define a ‘model-driven acquisition’ approach, where biological hypotheses will be translated into modifications of the cytomechanical model parameters values. Associated model simulations will identify the corresponding cell responses (membrane extension rate, number and frequency of membrane extensions, cell velocity...) and thus the range of experimental conditions within which such theoretical cell responses could be potentially observed. For example, the cell response and adaptation to an extracellular factors can be investigated in connection with the results of Raucher and Sheetz 〚19〛, who reported a greater than expected increase of lamellipodial extension rate in PDGF-activated cells. This was interpreted as the result of a localised decrease of the apparent membrane tension, compared to the more global effects of detergent addition, like desoxycholate, which still decrease in a similar way the membrane tension but initiate cell membrane extensions at many new sites. The cytomechanical model here presented can take into account such global shape changes, while distinguishing between different parameters which determine the apparent membrane tension, including the bilayer tension and the membrane-cytoskeleton adhesion modulated by the F-actin concentrations.

Beyond its ability to describe the local but globally coordinated nature of lamelipodial activity, the considered cytomechanical model could thus provide a common basis for understanding morphological changes observed with other cell lines or with cells submitted to different experimental conditions in the framework of nonlinear dynamical systems. For example, it has been reported that slight periodic variations of the pH can induce periodic modification of the cell shape 〚21〛. More generally, a given cellular dynamical state corresponds to a given set of cytomechanical parameters (membrane elasticity, cytoplasm viscosity...) that could differentiate the various cell lines 〚11〛. The different self-organised behaviours experimentally observed on isolated fibroblasts, as well as on other cells 〚4, 12, 13〛 favour this hypothesis. But further data are needed to assess that, in addition to the molecular characterisation based on the use of specific fluorescent probes, quantification of spontaneous membrane deformation could be used as a specific cell line signature.

## Acknowledgements

A.S. is supported by a fellowship from the ‘Région Rhône-Alpes’. This work is granted by the program ‘Métrologie du vivant’ from the ‘Région Rhône-Alpes’. We thank Dr Xavier Ronot and Anne Doisy for providing videomicroscopy sequences of cultured fibroblasts and for helpful discussions.

## Version abrégée

Les changements de forme cellulaire sont dus à la réorganisation de l'actine filamenteuse (actine-F), un composé complexe et très dynamique, contenu à la périphérie du corps cellulaire dans les cellules eukaryotes. Cependant, les mécanismes à l'origine de la formation des protrusions cellulaires sont toujours controversés. Les observations expérimentales des déformations cellulaires spontanées sur des substrats bi-dimensionnels montrent que les cellules manifestent un comportement auto-organisé, avec l'apparition de schémas de déformation récurrents. Un modèle cytomécanique, basé sur les propriétés viscoélastiques du réseau d'actine-F et sur les propriétés mécaniques de la membrane cellulaire, a permis de simuler la plupart des comportements dynamiques spontanés observés expérimentalement.

Nous nous sommes alors appuyés sur ce modèle pour étudier les interactions cellule–cellule et la migration cellulaire par chimiotaxie. L'observation expérimentale de fibroblastes L929 isolés montre que ces cellules ont spontanément un mouvement cyclique pulsant d'extension/rétraction de la membrane, selon deux directions perpendiculaires. Cette dynamique spontanée est altérée dans le cas des cellules en interaction : les cellules perdent leur configuration symétrique et semblent adopter de nouvelles morphologies limitant les contacts cellule–cellule. À partir de ces observations, nous proposons de prendre en compte ces interactions par une modification locale des paramètres du modèle cytomécanique : la pression hydrostatique protrusive et la tension de surface induite par la courbure membranaire. La localisation de la zone d'interaction par rapport à la géométrie symétrique initiale de la cellule a également été analysée. On observe que la dynamique de déformation est altérée dans les zones d'interactions, tandis qu'une transition dynamique entre une déformation membranaire oscillante stationnaire et une onde de déformation en rotation a été générée.

Le modèle a ensuite été étendu à l'étude de la migration cellulaire induite par un gradient de facteur extracellulaire. Les particules du facteur extracellulaire interagissent avec la membrane et le cortex cellulaire, en modifiant leurs propriétés mécaniques. Ces interactions induisent une instabilité morphologique, qui conduit à la migration de la cellule. Celle-ci remonte alors le gradient extracellulaire selon un processus impliquant deux étapes principales : (i) la polarisation de la cellule dans la direction du gradient, puis (ii) le déplacement de son centre de masse lorsque qu'une polarisation suffisante est atteinte. Ce seuil représente la force nécessaire à la cellule pour contrebalancer sa force d'adhésion au substrat. Trois simulations, associées à des situations expérimentales concrètes, sont proposées : migration de la cellule vers une source fixe, migration de la cellule avec déplacement de la source et migration de la cellule avec extinction de la source. Dans les trois cas, le comportement de la cellule simulée est en accord avec l'expérience. Deux paramètres de contrôle de la vitesse de migration sont ensuite analysés. Une augmentation du coefficient de diffusion du facteur extracellulaire conduit à un ralentissement de la migration. Une augmentation du seuil caractérisant la force d'adhésion de la cellule à son substrat conduit également à un tel ralentissement.