Gilberto L Thomas

Ismael Fortuna

Gabriel C Perrone

François Graner

Rita M C De Almeida

Rita M C De

Shape-velocity correlation defines polarization in migrating cell simulations

Keywords: single-cell migration, Cell polarization, Modified Fürth equation, CompuCell3D PACS: 05.40.-a, 87.17.Aa, 87.17.Jj

Cell migration plays essential roles in development, wound healing, diseases, and in the maintenance of a complex body. Experiments in collective cell migration generally measure quantities such as cell displacement and velocity. The observed short-time diffusion regime for mean square displacement in single-cell migration experiments on flat surfaces calls into question the definition of cell velocity and the measurement protocol. Theoretical results in stochastic modeling for single-cell migration have shown that this fast diffusive regime is explained by a white noise acting on displacement on the direction perpendicular to the migrating cell polarization axis (not on velocity). The prediction is that only the component of velocity parallel to the polarization axis is a well-defined quantity, with a robust measurement protocol. Here, we ask whether we can find a definition of a migrating-cell polarization that is able to predict the cell's subsequent displacement, based on measurements of its shape. Supported by experimental evidence that cell nucleus lags behind the cell center of mass in a migrating cell, we propose a robust parametrization for cell migration where the distance between cell nucleus and the cell's center of mass defines cell shape polarization. We tested the proposed methods by applying to a simulation model for three-dimensional cells performed in the CompuCell3D environment, previously shown to reproduce biological cells kinematics migrating on a flat surface.

Introduction

Polarization plays a key role in cell migration, both regarding single-cell movement on a flat substrate, or cells collectively moving in a tissue. In fact, polarization is important in every active matter model [START_REF] B Ten Hagen | Brownian motion of a self-propelled particle[END_REF][START_REF] Romanczuk | Active brownian particles -from individual to collective stochastic dynamics[END_REF], since the particle's activity necessarily breaks spatial symmetry. Although there are many different polarization definitions, as the one relative to the basal-to-apical axis in epithelia, in this paper we focus on cell polarization related to the rear-to-front axis associated to eukaryotic cell migration. Here we pay special attention to its correlation to cell velocity. In what follows, we carefully define what we mean by cell polarization and its measurement procedure, as well as how we define and measure cell velocity from a succession of snapshots of cell position obtained from numerical simulations. Our proposition for polarization measurement procedure should serve to infer cell velocity modulus and direction in both simulations and experiments.

Migrating cell polarization has its roots in the biological organization inside the cell and hence can be taken as a possible indicator of the cell movement. In fact, Maiuri and collaborators showed that there is a Universal Coupling between Speed and Polarization (UCSP) [START_REF] Maiuri | Actin flows mediate a universal coupling between cell speed and cell persistence[END_REF], such that both polarization and speed are correlated. UCSP is based on the fact that cell speed and persistence depends on a general advection of polarity ingredients, such as myosins or other molecules, with retrograde flow of actin in respect to the cell's center of mass (for a review, see references [START_REF] Insall | Actin dynamics at the leading edge: From simple machinery to complex networks[END_REF][START_REF] Ridley | Life at the leading edge[END_REF]). In this case, there must be a link between cell polarization and kinematic measurements such as MSD curves or velocity auto correlation functions (VACF). To find this link, it is necessary to propose a robust cell speed measurement and provide a precise definition for cell polarization, that could be used in both single-cell and collective migration as a quantitative predictor for cell movement. Here we aim at defining a migration-related polarization that could serve as a predictor for cell migration.

Velocity definition requires care. Reports on quantitative experiments on single-cell migration on a flat substrate date back to 1920 when Fürth [START_REF] Fürth | Die brownsche bewegung bei berücksichtigung einer persistenz der bewegungsrichtung. mit anvendungen auf die bewegung lebender infusorien[END_REF] showed that the mean square displacement (MSD), represented by |∆ r| 2 , is obtained from the stationary solution of an isotropic Ornstein-Uhlenbeck problem for cell velocity [START_REF] Uhlenbeck | On the theory of the brownian motion[END_REF] now known as Fürth equation:

|∆ r| 2 = 4D [∆t -P (1 -exp (-∆t/P ))] , (1)
where ∆t stands for the time interval used to obtain the squared cell displacement |∆ r| 2 . The symbols represent averages over different times and trajectories, after the stationary state is reached. D stands for the diffusion coefficient of the cell, reflecting the fact that for long ∆t we have |∆ r| 2 ∼ 4D∆t. This equation represents the sum over the displacements measured in two equivalent, perpendicular directions. On the short-time interval limit, Eq. 1 yields a ballistic regime for which |∆ r| 2 ∼ 2D P ∆t 2 . Here, P stands for the time scale associated to the transition from ballistic to diffusive regime, and is known as the persistence time of the model. This equation is obtained from the stationary, exact solution for a stochastic, 2D, isotropic Langevin problem for cell velocity (when the memory on the initial state has been lost), that may be written as

d v dt = - γ m v + η(t), (2)
where γ represents the friction of a particle moving in a fluid, m the particle mass and η(t) is a white, 2Dvector noise, whose components have a standard deviation given by parameter g, defined by η i (t)η j (t ) = gδ ij δ(t -t ). Here, P = m γ and D = g 2mγ . Observe that this equation is isotropic in space. Hence, it should not come as a surprise that a direct consequence of the anisotropy of migrating cells comes in the form of a deviation in the MSD equation observed in both experiments [START_REF] Thomas | Parameterizing cell movement when the instantaneous cell migration velocity is ill-defined[END_REF] and simulations [START_REF] Fortuna | Compucell3d simulations reproduce mesenchymal cell migration on flat substrates[END_REF], which follows instead the modified Fürth equation:

|∆ r| 2 = 2D [∆t -P (1 -exp (-∆t/P ))] + 2DS 1 -S ∆t , (3)
where 0 ≤ S < 1 is a non-dimensional parameter. The modified Fürth equation presents a further diffusive regime for short time intervals, since for short timescales |∆ r| 2 ∼ 2DS 1-S ∆t, indicating an additional, fast diffusive regime. SP gives the time scale for which the systems present a transition from this fast, short time interval diffusive regime to a ballistic-like regime, with P signaling a second regime transition, from the ballistic-like to a slow, long time interval diffusion. By defining natural units for time (P ) and length ( 2DP/(1 -S)), Eq. 3 may be recast in terms of non-dimensional quantities |∆ ρ| 2 = |∆ r| 2 (1-S)/2DP and ∆τ = ∆t/P as follows

|∆ ρ| 2 = ∆τ -(1 -S) [1 -exp(-∆τ)] . (4)
The short-time interval diffusion implies that instantaneous velocity is not a well-defined quantity [START_REF] Thomas | Parameterizing cell movement when the instantaneous cell migration velocity is ill-defined[END_REF], and the Langevin problem, written in the form of a time derivative for particle velocity, cannot yield the modified Fürth equation, Eq. 4. De Almeida et al. [START_REF] Rita | Exact solution for an anisotropic ornstein-uhlenbeck process[END_REF] proposed a set of equations where in a given instantaneous direction (that we shall call polarization direction) the dynamics is equivalent to a one dimensional Langevin problem, while in the perpendicular direction(s) the dynamics is driven by a Gaussian noise in the equation for displacement (not velocity), that is, a translational noise. A dynamic equation for the angle describing the direction of the Langevin dynamics completes the model. In Ref. [START_REF] Rita | Exact solution for an anisotropic ornstein-uhlenbeck process[END_REF] the authors analytically obtain the asymptotic solution and show that the modified Fürth equation is the exact form for the MSD. In this model, the short-time diffusive motion comes from the translational noise in the direction perpendicular to polarization, that is, the deviation from the pure Fürth MSD behavior observed in both simulations and experiments is a direct consequence of the spatial-symmetry break in the kinematics of migrating cells.

The consequence of a diffusive motion at short timescale is that as lim ∆t→0 ∆ r ∆t does not converge [START_REF] Rita | Exact solution for an anisotropic ornstein-uhlenbeck process[END_REF]: an experimentalist may get different values for ∆ r ∆t , depending on the value of ∆t chosen to measure ∆ r, regardless how small ∆t is.

In a recent paper, Fortuna and collaborators [START_REF] Fortuna | Compucell3d simulations reproduce mesenchymal cell migration on flat substrates[END_REF] proposed a simulation model for a three dimensional cell migrating on a flat substrate, using the cellular Potts model (CPM) [START_REF] Hogeweg | Computing and organism: on the interface between informatic and dynamic proceses[END_REF][START_REF] Graner | Simulation of biological cell sorting using a two-dimensional extended potts model[END_REF][START_REF] Glazier | Simulation of the differential adhesion driven rearrangement of biological cells[END_REF] in which, like in an experimental image, each cell is represented by a connected set of pixels. This model has successfully reproduced the migration kinematics observed in experiments [START_REF] Thomas | Parameterizing cell movement when the instantaneous cell migration velocity is ill-defined[END_REF]. The computational CC3D simulation project, along with the instructions to run it, can be downloaded from [START_REF]CC3D simulation project[END_REF]. Here we propose using this model to analyze different geometrical definitions of cell polarization as related to migration. We validate the proposed definitions by their correlation to the cell mobility. Next section defines quantities and their measurement procedures, section III presents the results, and in section IV we discuss and present our conclusions.

Definitions and methods

Simulation model

We use the cellular Potts model implemented in the open source CompuCell3D (CC3D) simulation environment [START_REF] Swat | Multi-scale modeling of tissues using compucell3d[END_REF]. We build a cell as composed of three compartments, each being a connected cluster of sites with an identifying label σ, that is related to the cell (the medium and substrate are treated by the algorithm as special types of cells [START_REF] Fortuna | Compucell3d simulations reproduce mesenchymal cell migration on flat substrates[END_REF]). This feature makes the CPM unique when trying to simulate detailed cell shape contours and fluctuations. A second label, C, is used to identify the cell compartment: nucleus (C = 1), cytoplasm (C = 2), or lamellipodium (C = 3). In this an energy-like function is written as

E = E interf ace + E volume (5)
where

E interf ace = r s( r) J ( r, s) , (6)
where J ( r, s) depends on the site labels σ and C, and is the interface energy per lattice-site surface between neighboring lattice sites at sites at r and s. The sum over s( r), runs over the 4th-neighbor range around r (32 neighbors) to reduce lattice anisotropy [START_REF] Holm | Effects of lattice anisotropy and temperature on domain growth in the two-dimensional potts model[END_REF]. J = 0 is assumed for neighboring lattice sites that belong to the same cell and compartment. In all other cases, we set interface energies as positive (i.e. the Js are ferromagnetic). We choose J values to ensure that the cytoplasm always surrounds the nucleus and the lamellipodium remains attached to the cytoplasm, substrate and medium. The volume term on the right hand side (r.h.s.) of Eq. 5 constrains the volumes of the cell compartments to be close to their reference volumes:

E target volume = 3 C=1 λ C V C -V target C 2 , (7)
where V C is the current volume of the C-th cell compartment, V target C is its target volume, and λ C is the inverse compressibility of the compartment.

Configuration changes are performed as usual in Monte Carlo simulations: after randomly picking a pair of nearest neighbors cell-lattice sites, we tentatively copy labels σ and C from the first to the second lattice site and calculate the change in energy ∆E. If ∆E ≤ 0, we accept the change, and if ∆E > 0, we accept the change with probability exp(-∆E/T m ), where T m is a parameter describing the amplitude of cell-membrane fluctuations. A Monte Carlo Step (MCS) consists of a number of copy attempts equal to the number of sites in the lattice.

Non-conservatives forces require some modification to the usual CPM dynamics. We add a further term when calculating ∆E, to extend the usual CPM effective energy to describe the protrusive forces that Factin polymerization exerts on the leading edge of the cell. After choosing a site r and its neighboring copy target site s, we add the following term to ∆E:

∆E F -actin = λ F -actin [F ( s) -F ( r)] , (8)
where the sum is taken only on the boundary between lattice sites lying in the lamellipodium compartment and the surrounding medium. By changing ∆E and consequently changing the probability that a lamellipodium site is copied over a medium site, Eq. 8 simulates a force to the membrane in the direction of the copy attempt (for more details, see Ref. [START_REF] Fortuna | Compucell3d simulations reproduce mesenchymal cell migration on flat substrates[END_REF]). F ( r) is a field specifying the concentration of F-actin, and λ F -actin > 0 is the force per area per unit F-actin field. The actin concentration field obeys the simplified reaction-diffusion equation, as described in Ref. [START_REF] Fortuna | Compucell3d simulations reproduce mesenchymal cell migration on flat substrates[END_REF], that is,

∂F ( r, t) ∂t = D F ∇ 2 F ( r, t) + k s δ(C( r, t) -3) -k d F ( r, t), (9)
where D F is an effective diffusion constant, k s(d) represents the polymerization (depolymerization) rate of F-actin fibers, and the δ-function says that only lamellipodium creates them. This equation is solved during the simulations by a tool included in the CC3D environment. Details on CC3D PDE solvers can be read in [START_REF] Maciej | Compucell3d reference manual -4[END_REF]. The final form assumed by the F-actin field is illustrated in Fig. 1: the parameters choice ensures that F-actin field is nearly constant on lamellipodium sites and zero otherwise.

The term for ∆E F -actin together with Eq. 5 i) increase the rate at which lamellipodium is copied over cytoplasm; ii) with the lamellipodium volume increase, the copy of cytoplasm sites over lamellipodium is favored (due to lamellipodium target volume term), and iii) the same happening to the rate with medium lattice sites overwrite cytoplasm lattice sites, creating a polarization axis that drives cell migration in the direction of the lamellipodium compartment . This polarization and its consequences simulate the effect of substrate forces and the interaction with extra-cellular matrix through adherent complexes in real cells, in a qualitative and simplified way. For a more detailed model that considers elastic substrates and the explicit role of actomyosin fibers and adhesion complexes, we refer the reader to the works by Zhong and collaborators [START_REF] Zhong | Impact of cell shape on cell migrationbehavior on elastic substrate[END_REF][START_REF] Zhong | How do cells produce and regulate the driving force in the process of migration?[END_REF] and He and collaborators [START_REF] He | Some basic questions on mechanosensing in cell-substrate interaction[END_REF] .

Figure 1 illustrates the crawling dynamics in an in silico cell, and Figure 2 shows two snapshots of them, as created by the CC3D, migrating cells (from Ref. [START_REF] Chen | Single-cell migration chip for chemotaxis-based microfluidic selection of heterogeneous cell populations[END_REF] ), depicting different morphologies that they may assume. Observe that lamellipodium may be engulfing the cytoplasm or concentrated at one side of the cell. In the first case, lamellipodium tends to grow over medium in all directions, producing a short timescale diffusive regime. In the second situation, lamellipodium has a preferred direction to grow over medium sites and creates the polarization axis illustrated in 1 that induces a preferred motion direction, contributing to a non-diffusive motion.

Larger values of λ F -actin favors localizing lamellipodium at one side of the cytoplasm, due to a local-excitation, global-inhibition loop: the more lamellipodium at one side of the cytoplasm, the more lamellipodium grows there (local excitation) and when more lamellipodium is created, the target volume term inhibits the creation of new lamellipodium sites everywhere (global inhibition). Consequently, larger values of λ F -actin favors polarization, decreasing the time interval where it is possible to detect the short time diffusive regime, scaled by SP . The cell compartments are: lamellipodium in green, cytoplasm in blue, and nucleus in yellow. Three steps create a polarization axis that favors migration leftwards: 1 -The F actin field, which effectively non-zero at lamellipodium lattice sites (red curve), favors lamellipodium lattice sites to be copied over medium lattice sites. This increases lamellipodium compartment volume and extends it to the left. 2 -Since there is a volume constraint in the lamellipodium compartment, there is now a preferential copy of cytoplasm lattice sites over lamellipodium ones, which decreases lamellipodium compartment volume and increases cytoplasm compartment volume. 3 -Again there is a volume constraint on the cytoplasm compartment, which favors medium lattice sites to be copied over cytoplasm compartment ones.

Polarization Definitions and Measurements

As we show in this paper, the intuitive idea that front-rear cell polarization should be able to predict cell movement critically depends on how we define this movement. Displacement and instantaneous velocity are two possibilities that are interchangeable when velocity is a well defined vector quantity, which frequently is not the case for migrating cells. Instantaneous velocity is obtained by taking the limit of the displacement per unit time interval ∆t as ∆t → 0. This limit does not exist when the short time interval diffusive term is present (see Eq. 3). In discussing the polarization correlation with cell movement we must then use a mean velocity, where the time interval is kept at a finite value. Because of that, in what follows, we use the terms displacement, velocity and mean velocity meaning different quantities, that cannot be used to characterize the same measurement. Also, we analyze both direction and modulus of mean velocity correlations with polarization. We use the term mean speed for mean velocity modulus, keeping the term velocity for the corresponding vector quantity.

To define a migration-related polarization that could serve as a predictor for cell migration, the choice of the candidate quantities to be further analyzed is based on two biological observations: i) there is a universal coupling between speed and polarization as described by Maiuri and collaborators [START_REF] Maiuri | Actin flows mediate a universal coupling between cell speed and cell persistence[END_REF], that was later applied to explain the retrograde actin flow in migrating cells [START_REF] Andrew C Callan-Jones | Actin flows in cell migration: from locomotion and polarity to trajectories[END_REF]; and ii) the nucleus is positioned at the rear of a migrating cell [START_REF] Gundersen | Nuclear positioning[END_REF]. Both observations are explained by the cytoskeleton arrangement in migration, where acto-myosin fibers and microtubules help in both dragging the cell nucleus and in the trafficking of the necessary proteins, creating the required protein gradients that support the trailing edge and rear contraction, responsible for cell migration.

Here we consider the following possible measures for the cell polarization Π, that explicitly consider the nucleus, to grasp the experimental observation by Gundersen and Woman [START_REF] Gundersen | Nuclear positioning[END_REF] that in migrating cells the Snapshots of simulated cells. A (view from above) and A' (vertical cross section) represent a cell whose spatial symmetry has not been broken. The centers of mass of nucleus (yellow), cytoplasm (blue), and lamellipodium (green) almost coincide (red cross). In A', orange, light gray, and green crosses locate the centers of mass of, respectively, nucleus, cytoplasm, and lamellipodium. B (view from above) and B' (vertical cross section) represent a cell with the broken spatial symmetry. In B, the projection on the substrate plane of the centers of mass do not coincide, and Nucleus lags behind. The centers of mass of these three compartments can be combined to produce different definitions of polarization vectors. In the insets, for comparison, experimental images from [START_REF] Chen | Single-cell migration chip for chemotaxis-based microfluidic selection of heterogeneous cell populations[END_REF].

nucleus lags behind:

Π L-N (t) = r L (t) -r N (t) xy , Π L-CN (t) = r L (t) -r CN (t) xy , Π C-N (t) = r C (t) -r N (t) xy , (10)
Π CN -N (t) = r CN (t) -r N (t) xy , Π C-CN (t) = r C (t) -r CN (t) xy ,
where the indices C, L, and N refer to the different cell compartments (Cytoplasm, Lamellipodium, and Nucleus) and CN refers to the combined compartment comprising Cytoplasm and Nucleus (Fig. 2), and r β (t) for β = L, C, N, CN refers to the center of mass of these compartments. Polarization vectors Π β , for β = L -N , L -CN , C -N , CN -N , and C -CN are defined as the difference vector between the centers of mass of two compartments, projected on the substrate, as indicated by the subindices in the right hand side of Eqs. 10. At each MCS these quantities are recorded. We recall that by knowing C and N , one can directly reconstruct the combined compartment CN comprising both Cytoplasm and Nucleus, and for instance:

Π C-N (t) = V C + V N V C Π CN -N (t) and Π C-CN (t) = V N V C Π CN -N (t), (11)
where V C and V N are the compartment volumes.

From the above equations, we have that Π C-CN (t) = Π C-N (t) -Π CN -N (t). Consequently, Π C-N , Π CN -N and Π C-CN contain the same information. However, in experiments we may easily assess position of the center of mass of the combination of cytoplasm and nucleus by the cell's geometrical center. In what follows we present results considering four definitions for cell polarization, disregarding Π C-CN .

To estimate the performance of the polarization measures as a proxy to cell displacement, we must estimate the correlation between these vectors and the subsequent cell displacement. However, we must first choose a suitable time interval. For too small time intervals (smaller than SP ), displacements may be dominated by random fluctuations in the direction perpendicular to cell polarization [START_REF] Thomas | Parameterizing cell movement when the instantaneous cell migration velocity is ill-defined[END_REF][START_REF] Fortuna | Compucell3d simulations reproduce mesenchymal cell migration on flat substrates[END_REF][START_REF] Rita | Exact solution for an anisotropic ornstein-uhlenbeck process[END_REF]. For too large time intervals, larger that the persistence time P , displacement may also decorrelate with cell polarization. Hence, to evaluate the performance of cell polarization definitions as predictors of cell displacement, we must first choose an adequate time interval. We shall define as ∆τ opt the time interval for obtaining cell displacement that has the highest correlation with cell polarization.

A theoretical estimate of ∆τ theor opt is obtained by finding the time interval that maximizes the slope of the log-log curve |∆ ρ| 2 , given in Eq. 4. The slope α is given by

α = d log( |∆ ρ| 2 ) d log(∆τ) = ∆τ(1 -(1 -S)e -∆τ ) ∆τ -(1 -S)(1 -e -∆τ ) . (12)
The inclination α assumes values in the interval [START_REF] B Ten Hagen | Brownian motion of a self-propelled particle[END_REF][START_REF] Romanczuk | Active brownian particles -from individual to collective stochastic dynamics[END_REF], passing through a maximum for S ≤ ∆τ ≤ 1, that is obtained by numerically solving a transcendental equation. We assume that the value of ∆τ that maximizes α defines a suitable time interval to perform correlation measurements. To test this hypothesis, we must compare with adequate correlation measurements as follows.

To determine the time intervals for which the short-time random fluctuations dominate the dynamics, we monitor mean speed u(τ, δ), defined as the modulus of the mean velocity u(τ, δ) for a finite time interval δ as

u(τ, δ) = ρ(τ + δ) -ρ(τ) δ , ( 13
)
where ρ is the cell position in natural units, and the average of both mean velocity u δ and mean speed u δ are taken over all time points τ and different runs, after the stationary state has been reached. Since for small δ we expect the random, short-time noise to dominate the dynamics, we also expect u δ to diverge as δ → 0. On the other hand, when δ = ∆τ theor opt we expect that the mean speed and polarization correlate. At each MCS, we recorded the center of mass of each compartment and of the combined cytoplasm and nucleus compartment together with the different polarization definitions, as given by Eqs. 10. From there it is possible to obtain the angle between polarization at a given time and the mean velocity in the subsequent time interval δ. We monitor the best correlation between the directions of cell polarization and cell mean velocity by monitoring the velocity direction -polarization correlation index ξ β for each polarization β defined in Eqs. 10, defined as

ξ β = cos θ β sin θ β , (14)
where • represents average over the whole trajectory (after the stationary state has been reached) and over different simulation runs. Index β stands for the labels CN -N , L -CN , and L -N so that θ β represents the angle between cell displacement and the β-th polarization direction (C -N is colinear to CN -N and is not explicitly taken into account here). Here ξ β depends on the time interval δ used to obtain cell mean velocity and its values range from 0, when polarization and the mean velocity u(τ, δ) are always orthogonal, to ∞, when polarization and the mean velocity u(τ, δ) are always parallel.

Besides correlation between the adequate measure for polarization and mean velocity modulus and direction, we also assessed the dispersion around the predicted value for both velocity modulus and direction. In this way, besides predicting the expected value for u δ , we also give the uncertainty we should expect from that prediction. For the modulus, we propose a relative error, ε, given as

ε = σ speed u ∆τ theor opt , (15)
where σ speed is the standard deviation of | u(τ, ∆τ theo opt )| over the trajectory. We monitor the dispersion in direction using the standard deviation of θ β , denoted σ θ β , since the average θ β is close to 0 and hence cannot provide a scale for dispersion. Anyway, angles have a finite natural scale of 2π. Table S1 in Supplementary Materials Online lists the main variables and parameters used in this work.

Results

We first verify whether the polarization definitions presented in Eqs. 10 may be used as predictors for cell velocity by analysing our simulation results. In this paper, we keep fixed the parameters associated with energy and temperature as in Ref. [START_REF] Fortuna | Compucell3d simulations reproduce mesenchymal cell migration on flat substrates[END_REF] and varied parameters related to actin-force intensity (λ F -actin = 150, 175, and 200) and the lamellipodium volume fraction (φ = 0.05, 0.10, 0.20, and 0.30). For each set of parameters, we run 10 simulations with 10 5 MCS. We also considered three different cell sizes, R cell = 10, 15, and 20 measured in (lattice sites) 1/3 , to test for finite size effects. The behavior is qualitatively the same and in the main text we considered R cell = 15 while all results concerning R cell = 10, and 20, are shown in Supplementary Materials Online. Table S2 in Supplementary Material Online shows the parameter sets and fitting constants. For more details in the simulations, we refer the readers to Ref. [START_REF] Fortuna | Compucell3d simulations reproduce mesenchymal cell migration on flat substrates[END_REF].

The simulations start with a spherical cell with no lamellipodium engulfing a cubic nucleus resting over the substrate. The dynamics of the model is such that the cytoplasm sites touching the substrate may transform into lamellipodium sites until the lamellipodium target volume is attained. Due to the energy parameters J between sites belonging to different compartments or medium and substrate, the lamellipodium spreads over the substrate, forming a thin layer that surrounds cytoplasm, as shown in Fig. 2 A (see also Fig. 4 from Ref. [START_REF] Fortuna | Compucell3d simulations reproduce mesenchymal cell migration on flat substrates[END_REF]). For the parameters chosen for this manuscript, the lamellipodium eventually assumes a configuration where it localizes at one side of the cell, spontaneously breaking symmetry and initiating migration. The asymmetric lamellipodium configuration drives migration and its continuous direction changes is reponsible for the observed long time interval diffusion. Usually this behavior has been characterized by calculating the MSD. The transient from the initial configuration is detectable by the relaxation of the MSD curve to a stationary form, when the memory from the initial configuration is lost. Here, we have waited 1000 MCS before recording the data that we use in the analyses. After this simulation time MSD curves do not change any more. For more details in the simulations, we refer the readers to Ref. [START_REF] Fortuna | Compucell3d simulations reproduce mesenchymal cell migration on flat substrates[END_REF].

Figure 3 shows the results for a selected set of parameters, as indicated. In panel A, the mean square displacement |∆ ρ| 2 is plotted versus time interval (∆τ), both given in natural units, together with the fits using the modified Fürth equation (Eq. 4). The deviation from the original Fürth equation (dashed, light gray line) for small time intervals is clearly visible. In panel B, we present the average mean speed u δ as a function of δ. Observe that when δ is too small, mean speed diverges as expected, signaling a diffusive dynamics. In both panels A and B, the dashed vertical lines indicate the value of ∆τ theor opt for which the corresponding |∆ ρ| 2 (in panel A) presents the maximum exponent α opt , as predicted by maximizing Eq. 12. Values for ∆τ theor opt for all simulations are presented in Table S1 in Supplementary Materials Online, where we also present figures for all parameter sets similar to panels A and B (Figs. S1-S3). First it shows that for φ = 0.10 and λ F -actin = 200 the polarization definition whose velocity-polarization correlation index ξ β has its maximum more coincident with the theoretical prediction ∆τ theor opt (indicated by the pink arrows) corresponds to either Π C-N (t) or Π CN -N (t).

Second, the case φ = 0.30 and λ F -actin = 150 requires special care. For this case S = 9.46 × 10 -2 , implying i) a small value for the maximum exponent α of their MSD curves (see Fig. 3 A), ii) a short, almost nonexistent plateau for u δ (see Fig. 3 B) , and iii) a short ballistic-like regime in the MSD curves, also shown in Fig. 3 A. These three observations suggest that the cell never reaches a proper ballistic phase. In fact, the weak intensity of λ F -actin and the large lamellipodium fraction reduces the probability of a symmetry break and the cell preferentially stays in a configuration as shown in Fig. 2 A. The persistence time is hence not large and the time intervals for the fast and slow diffusive regimes are similar. It is then not surprising that all three definitions of polarization fail to yield a velocity-polarization correlation index whose maximum is near the predicted value of ∆τ theor opt . Figs. S1-S6 in Supplementary Materials Online present the same information in Figs. 3 and4 for all parameter sets. (magenta) as functions of S. We have used all fitted Svalues from Ref. [START_REF] Fortuna | Compucell3d simulations reproduce mesenchymal cell migration on flat substrates[END_REF] to numerically obtain ∆τ theor opt that maximizes α in Eq. 12. We fitted these points as a power law. Figure 5 shows the fit as a magenta solid line, with dots representing the numerical solutions for those values of S corresponding to each simulation of the whole set, including R cell = 10 and 20. In purple we show the equivalent for the values of α opt . The numerically obtained ∆τ theor opt are well fitted by the power law ∆τ theor opt = 1.932(±8.6 × 10 -3 ) S 0.455(±2.0×10 -3 ) with R 2 > 0.999, while the points associated to α opt are fitted by α opt = 2.022(±2.39 × 10 -3 ) -1.115(±4.56 × 10 -3 )S 0.397(±3.90×10 -3 ) , with R 2 > 0.999. Observe that, as S approaches zero, ∆τ theor opt approaches zero and α opt approaches 2: when the short-time diffusive regime disappears, the optimal time interval goes to zero, as it should, such that the MSD curve presents a maximum exponent equal to 2. Also, as ∆τ theor opt increases, α opt decreases: when the fast diffusive regime is long (S is large), the MSD curve has a short ballistic-like regime, such that the maximum value of α opt does not reach values near 2. In fact, for S > 0.3 we have α opt < 1.3. The consequence of the behavior shown in Figs. 3 to 5 is that by fitting the MSD curve from simulations or experiments using the Modified Fürth equation, it is possible to obtain S and from there to obtain the time interval ∆τ theor opt for which the correlation between polarization direction and displacement direction is optimized. Based on the results shown in Fig. 4, we select Π CN -N to proceed with our analyses.

Figure 6 A shows the expectation value for mean speed u(∆τ theor opt ) and Fig. 6 B shows its relative error ε versus the average polarization modulus | Π CN -N | for 1980 points of 10 typical trajectories (19800 points), after the stationary state is reached. Each set of points is obtained using the same λ F -actin are linked with lines. In these lines, each point corresponds to a different lamellipodium fraction φ (0.05, 0.1, 0.2, and 0.3) increasing from right to left. This figure shows that u(∆τ theor opt ) increases with average polarization modulus, while the relative error in speed decreases: when | Π CN -N | is large, the cell migrates with higher speed and with less dispersion.

Figure 6 C shows the effect that increasing φ decreases the overall average of speed: as the cells have more lamellipodium, the symmetry break that allows the lamellipodium to efficiently drive cell migration is less frequent.

Figure 6 D shows all measured points for two different parameter sets: (φ = 0.3, λ F -actin = 150) that shows a very small polarization and (φ = 0.05, λ F -actin = 200), where the cells are highly polarized most of the time. The cloud of points for each case illustrates that increasing | Π CN -N | correlates with increasing u(∆τ theor opt ) and decreasing dispersion for more polarized cells. Figure 6 D also presents a picture of a cell in a typical configuration from the simulations for the two sets of parameters, evincing the different lamellipodium configurations for each case. The clouds presented in Fig. 6 D bring further information. For the parameter sets presenting a well defined ballistic-like regime, most of the time the cells present a well developed polarization and, hence, most of the points will present larger values of polarization. Cells will sporadically present small polarization when they are changing the direction of their lamellipodia. On the other hand, in case the cells mostly present small values of polarization, most of the time they are not in a ballistic-like regime. For φ = 0.30 and λ F -actin = 150, most points present small values of polarization. Observe that this happens in spite of a large value for φ, when too much lamellipodium engulfing the cytoplasm makes it more difficult to build a high polarization. Accordingly, this parameter set also shows a small value for the maximum exponent of the MSD curve, a short plateau at the u δ versus δ plot, and a very small velocity-polarization correlation index. In Supplementary Materials Online Fig. S7 presents the equivalent of Fig. 6 for for R cell = 10 and 20 and Figs. S8,S9, and S10 present the clouds of points, similar to Fig. 6 B, for all sets of parameters. Also, Figs. S11, S12, and S13 show the histograms for u(∆τ theor opt ) Π β for all 4 polarization definitions and all different simulation sets, further detailing the information summarized in Fig. 6. Now we turn to the direction of cell velocity as compared to the polarization Π CN -N . Figure 7 A plots the standard deviation σ θ CN -N for the angle θ CN -N between the mean velocity calculated using ∆τ theor opt and the cell polarization at the begining of ∆τ theor opt as a function of the average polarization modulus Π CN -N for each parameter set with R cell = 15. As expected, σ θ CN -N decreases as Π CN -N increases: more polarized cells show less deviation of their mean velocity from the polarization axis. The inset in this figure present the clouds for two sets of simulations: for both sets the mean value for θ CN -N is zero (since the deviations from the polarization axis are symmetrical), but the spread of the deviation angles are reduced for highly polarized cells.

Figure 7 B shows that the average polarization modulus decreases with φ, reinforcing the conclusions drawn from Fig. 6: a large fraction of lamellipodium may not support a polarized front that would maintain a large polarization. To illustrate this point, we have produced two animations, SM1-15 0.05 200.mp4 and SM2-15 0.30 150.mp4, for the the two parameter sets considered in the inset of panel A, available in Supplementary Materials Online. Accordingly, the dispersion in direction increases with φ, as shown in Fig. 7 C. Finally, Fig. 7 D shows that the average mean speed u ,CN -N δ measured parallel to Π CN -N converges as δ → 0 and, hence, this component is well defined. In Supplementary Materials Online Figs. S14 and S15 present the equivalent of Fig. 7 for for R cell = 10 and 20 and Figs. S16, S17, and S18 present the clouds of points of the inset in Fig. 7 A for all sets of parameters.

Together, Figs. 3 to 7 show that it is possible to define a polarization vector for the cells, and that it may serve as a proxy for cell displacement, provided the cells are at least minimally polarized. They also show that the most adequate polarization definitions are those that do not consider the lamellipodium center of mass. Regarding β = CN -N or C -N , for all parameter sets there is a clear correlation between the size of amplitude of polarization and both mean speed (calculated using δ = ∆τ theor opt ) and polarization direction. We remark, however, that when the cells are in a poorly migrating phenotyped, as parameterized by φ, persistent motion is rare and the polarization is small. In these cases, the definition of a polarization as a proxy for cell displacement is meaningless.

The correlation of polarization with cell displacement depends on cell displacement being measured for a time interval equal to δ = ∆τ theor opt . In turn, the definition of ∆τ theor opt depends on assuming that cells kinematics may be characterized by a MSD curve given by the modified Fürth equation, Eq. 3 [START_REF] Fortuna | Compucell3d simulations reproduce mesenchymal cell migration on flat substrates[END_REF]. While Fürth equation represents the stationary solution of an isotropic Ornstein-Uhlenbeck problem for the cell velocity, the modified Fürth equation is the stationary solution for an anisotropic Ornstein-Uhlenbeck problem [START_REF] Rita | Exact solution for an anisotropic ornstein-uhlenbeck process[END_REF]. In this model, the cell's instantaneous vector velocity is an ill-defined quantity but the speed in the direction of the polarization is well-defined and follows a one dimensional Langevin problem. As shown by the divergence of speed as δ → 0 in Fig. 3 B, the present simulations are able to reach the ill-defined velocity regime. Since here we have defined a polarization vector, we may now verify whether the cell's velocity in the direction of polarization converges or not as the time interval used to calculate the displacement decreases. The speed in the direction of the polarization is defined as

u ,β (τ, δ) = Π β (τ) • ρ(τ + δ) -ρ(τ) δ , β = CN -N, L -CN, L -N, C -N , (16)
and its average, u ,β δ for a given value of δ, is shown in Fig. 7 D for β = CN -N . We may conclude that the speed in the direction of the polarization is a well defined quantity. Figs. S19, S20, and S21 in Supplementary Materials Online show equivalent plots for all sets of parameters.

Discussion and Conclusions

From the results we presented above, we reach the following conclusions. First, it is possible to define a polarization vector on the basis of its correlation with cell displacement. Although we have considered simulations, we propose cell polarization definitions that are readily transposed to experiments, as the distance between the centers of nucleus and cytoplasm. For our simulations, polarization definitions that do not consider the position or shape of lamellipodia showed a better performance as cell displacement predictors. In migrating cells, lamellipodium position determines the localization of the nucleus behind the cell's geometric center [START_REF] Gundersen | Nuclear positioning[END_REF]. Consequently both nucleus and lamellipodium positions carry information on the direction of the movement. Nuclear positions, however, are more persistent and show smaller fluctuations.

Second, the assessment of the performance for polarization definition as a proxy for cell displacement direction depends on the time interval used to measure displacement. The determination of the ideal time interval depends on the fitting of MSD curve using the modified Fürth equation [START_REF] Fortuna | Compucell3d simulations reproduce mesenchymal cell migration on flat substrates[END_REF][START_REF] Thomas | Parameterizing cell movement when the instantaneous cell migration velocity is ill-defined[END_REF]. The optimal time interval for obtaining cell displacement are shown in Fig. 5 for all values of S, provided the simulated or experimental MSD curve is well fitted by the modified Fürth equation. A detailed procedure for fitting experimental MSD curve is provided in [START_REF] Thomas | Parameterizing cell movement when the instantaneous cell migration velocity is ill-defined[END_REF] and its associated supplementary materials, where the fitting was performed for 12 different experimental setups from 5 different laboratories. A second, faster but less precise method is 1) Plot the mean square displacement curve (in a log-log plot) and then 2) find the time interval that yields the maximum exponent in that curve. That would give an estimate for the optimal time interval in laboratory units to correlate polarization and mean speed. Either way, to obtain the polarization, it is necessary to follow the cell geometrical center of mass (without lamellipodium) and the cell nucleus.

Third, the performance of the polarization vector is quantified by Figs. 6 and 7. The polarization vector is a good predictor for cell displacement provided it is not small. For cells that spend most of the time polarized, and hence have a migrating phenotype, our definitions apply well to the simulations. That must be verified in experiments.

Fourth, the simulations verified that velocity is well defined when measured parallel to the polarization direction, that is, the average mean speed in the polarization direction, u ,CN -N δ , converges to a defined value as δ → 0, in contrast to the average mean speed, u δ , that diverges in the same limit, as shown in Fig. 3B. This behavior is in agreement with the theoretical predictions of the Anisotropic Ornstein-Uhlenbeck model [START_REF] Rita | Exact solution for an anisotropic ornstein-uhlenbeck process[END_REF]. In the direction perpendicular to the polarization axis, the velocity is ill-defined, diverging when the time interval used for estimating displacement approaches zero. This must be verified in experiments.

The polarization definition and measurement that we propose can be applied to migration experiments of single cells on flat substrates. Migration is an important feature in characterizing cell phenotypes as, for example, in cancer where tumor malignancy correlates with cell migration capabilities, or where the tumor cell diversity may indicate different disease outcomes. In these and many other instances, a robust and reproducible protocol for characterizing migrating cells is a first, necessary step. A second step is to obtain robust distributions of quantities (properly measured), such that phenotypic classes may be proposed, there is the need of obtaining not only the value of quantities as S, D, P , or the mean speed ∆τ theor opt , but also their distributions in a given population. This is a work in progress.

The extension of these definitions and measurement procedures to collectively migrating cells may also contribute to the field. This is another work in progress and will be presented elsewhere.

Figure 1 :

 1 Figure 1: Crawling dynamics in an in silico cell.The cell compartments are: lamellipodium in green, cytoplasm in blue, and nucleus in yellow. Three steps create a polarization axis that favors migration leftwards: 1 -The F actin field, which effectively non-zero at lamellipodium lattice sites (red curve), favors lamellipodium lattice sites to be copied over medium lattice sites. This increases lamellipodium compartment volume and extends it to the left. 2 -Since there is a volume constraint in the lamellipodium compartment, there is now a preferential copy of cytoplasm lattice sites over lamellipodium ones, which decreases lamellipodium compartment volume and increases cytoplasm compartment volume. 3 -Again there is a volume constraint on the cytoplasm compartment, which favors medium lattice sites to be copied over cytoplasm compartment ones.

Figure 2 :

 2 Figure 2: Shape polarization definitions.

Figure 3 :

 3 Figure 3: MSD and average mean speed curves reveal departure from classic Fürth equation. (A) |∆ ρ| 2 versus time interval (∆τ) in comparison with Fürth's equation (grey dashed line). (B) u δ versus δ. The dots indicate the predicted value for the optimal time interval and maximum exponent for the MSD curves, corresponding to the values of S obtained for the simulations. The curves in (A) and (B) stand for different parameter sets; each curve is averaged over 10 simulation samples for R cell = 15, φ = 0.10, for λ F -actin = 150 (black squares), 175 (red circles), 200 (green triangles). The dashed vertical lines indicate corresponding values for ∆τ theor opt in both panels. Standard errors in (A) and (B) are of the order of 10 -2 or less of the respective values of |∆ ρ| 2 and u δ and are smaller than the symbols size.

Figure 4

 4 Figure4presents ξ β , defined in Eq. 14, as a function of δ, the time interval used to calculate cell displacement for three definitions of polarization and for two sets of simulations parameters. Π C-N is not shown because its direction is identical to Π CN -N direction. What this figure shows us two important points, as follows.First it shows that for φ = 0.10 and λ F -actin = 200 the polarization definition whose velocity-polarization correlation index ξ β has its maximum more coincident with the theoretical prediction ∆τ theor

Figure 4 :

 4 Figure 4: Π CN -N is strongest correlator with mean velocity when calculated using δ = ∆τ theor opt . ξ β as a function of the time interval δ used to obtain the displacement direction, considering the three polarization definitions: CN -N (green diamonds), L -CN (blue circles), and L -N (red triangles), for two sets of parameters with R cell = 15: φ = 0.10, λ F -actin = 200 (filled symbols, ∆τ theor opt = 0.048) and φ = 0.30, λ F -actin = 150 (open symbols, ∆τ theor opt = 0.672). The pink arrows indicate the ∆τ theor opt for both parameter sets. Standard errors are smaller than the symbols size. Similar panels considering different values of φ for R cell = 10, 15, and 20 are presented in Supplementary Materials Online.

Figure 5

 5 Figure 5 shows α opt (purple) and ∆τ theor opt

Figure 5 :

 5 Figure 5: Dots represent ∆τ theor opt and αopt as functions of S, as numerically obtained from the modified Fürth equation, by maximizing Eq. 12 and using the values of S obtained from the simulated MSD curves. Solid lines represent the fitting functions in the legend.

Figure 6 :

 6 Figure 6: Performance of Π CN -N as a predictor for cell speed. (A) u ∆τ theor opt versus Π CN -N in units of cell radius for the CN-N polarization definition Eq.(10). (B) Relative error for the speed ε, as defined in Eq.(15), versus Π CN -N . The pink arrows in (A) and (B) indicate the direction of increasing φ. (C) The average speed u ∆τ theor opt versus φ. Each curve in (A,B,C) has 4 points, one for each value of φ. All curves are for all parameter sets with R cell = 15. The legend in (A) applies to (B) and (C) as well. (D) Individual cell speed u ∆τ theor opt versus its polarization modulus Π CN -N measured for all points of the stationary trajectory and 10 different simulations, considering two parameter sets (representative simulated cells are shown). For R cell = 10 and R cell = 20 see Fig. S7 in Supplementary Materials Online.

Figure 7 :

 7 Figure 7: Performance of Π CN -N as a predictor for cell displacement direction. (A) The angular standard deviation σ θ CN-N versus average polarization modulus ΠCN-N(τ) in units of cell radius for all parameter sets with R cell =15. The displacements for the calculating θCN-N considered δ = ∆τ theor opt . The inset presents individual measurements considering all points of the stationary trajectory, for 10 simulation runs. (B) Average polarization modulus Π CN -N as a function of φ. (C) σ θ CN-N versus φ. (D) Plot of the average mean speed u ,CN -N δ measured parallel to Π CN -N as a function of δ. Equivalent panels for R cell = 10 and 20 are shown in Figs. S14 and S15 in Supplementary Material Online.

Acknowledgemnts

We thank Dr. Hélène Delanoë-Ayari, of Université Claude Bernard -Lyon 1, and Dr Andrew Callan-Jones, from Université Paris Diderot, for carefully reading and providing suggestions for the manuscript. This work has been performed under the CAPES-COFECUB agreement 88887/190134. We acknowledge partial support from Brazilian Agency Conselho Nacional de Desenvolvimento Científico e Tecnológico -CNPq.

Bibliography