Laure Laforgue

Arnold Fertin

Yves Usson

Claude Verdier
email: claude.verdier@univ-grenoble-alpes.fr

Val Érie

M Laurent

Efficient deformation mechanisms enable invasive cancer cells to migrate faster in 3D collagen networks

Cancer cell migration is a widely studied topic but has been very often limited to two dimensional motion on various substrates. Indeed, less is known about cancer cell migration in 3D fibrous-extracellular matrix (ECM) including variations of the microenvironment. Here we used 3D time lapse imaging on a confocal microscope and a phase correlation method to follow fiber deformations, as well as cell morphology and live actin distribution during the migration of cancer cells. Different collagen concentrations together with three bladder cancer cell lines were used to investigate the role of the metastatic potential on 3D cell migration characteristics. We found that grade-3 cells (T24 and J82) are characterized by a great diversity of shapes in comparison with grade-2 cells (RT112). Moreover, grade-3 cells with the highest metastatic potential (J82) showed the highest values of migration speeds and diffusivities at low collagen concentration and the greatest sensitivity to collagen concentration. Our results also suggested that the small shape fluctuations of J82 cells are the signature of larger migration velocities. Moreover, the displacement fields generated by J82 cells showed significantly higher fiber displacements as compared to T24 and RT112 cells, regardless of collagen concentration. The analysis of cell movements enhanced the fact that bladder cancer cells were able to exhibit different phenotypes (mesenchymal, amoeboid). Furthermore, the analysis of spatio-temporal migration mechanisms showed that cancer cells are able to push or pull on collagen fibers, therefore producing efficient local collagen deformations in the vicinity of cells. Our results also revealed that dense actin regions are correlated with the largest displacement fields, and this correlation is enhanced for the most invasive J82 cancer cells. Therefore this work opens up new routes to understand cancer cell migration in soft biological networks.

Introduction

Cell migration is fundamental for many biological processes such as the immune response, wound repair, tissue homeostasis and also for pathologies such as cancer metastasis, inflammation or autoimmunity. The migration mechanisms have been widely studied but were initially limited to 2D motion on various substrates simplifying cell-substrate interactions [START_REF] Condeelis | Life at the leading edge : the formation of cell protusions[END_REF][START_REF] Lauffenburger | Cell migration: a physically integrated molecular process[END_REF][START_REF] Sheetz | Cell migration as a five-step cycle[END_REF] . Studies over the last decades have allowed to identify the different factors regulating 2D cell movement: substrate rigidity [START_REF] Lo | Cell movement is guided by the rigidity of the substrate[END_REF][START_REF] Ambrosi | Traction patterns of tumor cells[END_REF] , substrate anisotropy [START_REF] Aubry | Mechanical link between durotaxis, cell polarity and anisotropy during cell migration[END_REF] or the effect of extracellular protein density [START_REF] Palecek | Integrin-ligand binding properties govern cell migration speed through cell-substratum adhesiveness[END_REF] . These effects can be measured thanks to the substrate deformation (usually functionnalized gels) and subsequent forces obtained by Traction Force Microscopy (TFM) [START_REF] Dembo | Stresses at the cell-to-substrate interface during locomotion of fibroblasts[END_REF][START_REF] Butler | Traction fields, moments, and strain energy that cells exert on their surroundings[END_REF][START_REF] Peschetola | Time-dependent traction force microscopy for cancer cells as a measure of invasiveness[END_REF] . Following these works, two main migration modes have been identified in the literature, the mesenchymal and amoeboid types of motion [START_REF] Friedl | Plasticity of cell migration: a multiscale tuning model[END_REF] .

In 3D networks, efforts have been made thanks to the use of confocal microscopy, in particular fluorescent reflection techniques [START_REF] Friedl | Amoeboid leukocyte crawling through extracellular matrix: lessons from the dictyostelium paradigm of cell movement[END_REF][START_REF] Iordan | Breakdown of cell-collagen networks through collagen remodeling[END_REF] allowed to vizualise biological fibers -such as collagen -mimicking the physiological environment of the cells. Other studies showed that the two modes of migration (mesenchymal or amoeboid) observed on 2D substrates are also found during migration in 3D matrices [START_REF] Friedl | Tumour-cell invasion and migration: diversity and escape mechanisms[END_REF][START_REF] Torka | Rock signaling mediates the adoption of different modes of migration and invasion in human mammary epithelial tumor cells[END_REF][START_REF] Mierke | Physical view on migration modes[END_REF] . The proteolysis-independent amoeboid mode generally describes a way used by migrating cells constantly developping protrusions and retracting them by changing their shape, in order to move forward. Interestingly, cells do not necessary use adhesion molecules such as integrins to migrate in 3D [START_REF] Lämmermann | Rapid leukocyte migration by integrin-independent flowing and squeezing[END_REF] . Regarding the mesenchymal mode, cells usually exhibit a highly polarized shape with the development of protrusions, the formation and pulling on adhesions, cell contraction, release of bonds at the rear and recycling of proteins [START_REF] Friedl | Plasticity of cell migration: a multiscale tuning model[END_REF] . In this mode, the molecular mechanisms involved are based on actin polymerization and acto-myosin contractility, this being modulated by the level of attachment [START_REF] Lämmermann | Mechanical modes of 'amoeboid' cell migration[END_REF] , leading to traction force generation. In contrast to cells migrating on 2D substrates, large F-actin protrusions like lamellipodia are less present in 3D situations [START_REF] Friedl | Plasticity of cell migration: a multiscale tuning model[END_REF][START_REF] Mierke | Physical view on migration modes[END_REF][START_REF] Petrie | At the leading edge of three-dimensional cell migration[END_REF] . Moreover, cells cultured in 3D matrices show sparse stress fiber structures, more complex F-actin organisation at the front of the cell [START_REF] Olson | The actin cytoskeleton in cancer cell motility[END_REF] and the acto-myosin cytoskeleton is mainly localised in the sub-membranous actin cortex. Interestingly, some in vivo studies revealed the role of cytoskeleton regulators and adhesion receptors in pathological cell motility. Alterations in regulation pathways of acto-myosin contraction or actin dynamics contribute to the ability of tumor cells to switch between a variety of strategies [START_REF] Olson | The actin cytoskeleton in cancer cell motility[END_REF][START_REF] Clark | Modes of cancer cell invasion and the role of the microenvironment[END_REF] . Thereby, migrating tumor cells can adopt diverse morphologies and actin organization depending on cancer type [START_REF] Pandya | Modes of invasion during tumour dissemination[END_REF] . When adopting preferentially a migration mode or switching between different modes, cancer cells can adapt to their surrounding ECM, interact and modify the microenvironment for optimum motility [START_REF] Lautscham | Migration in confined 3d environments is determined by a combination of adhesiveness, nuclear volume, contractility, and cell stiffness[END_REF] . Moreover, the temporal fluctuation of cancer cell shapes can facilitate their migration [START_REF] Eddy | Morphodynamics facilitate cancer cells to navigate 3d extracellular matrix[END_REF] . Taken together, these results highlight the complexity of 3D migration and the need to bring new insights about the spatio-temporal cell behavior [START_REF] Wu | Statistical analysis of cell migration in 3d using the anisotropic persistent random walk model[END_REF] , in particular cell-induced 3D environment deformations. This is why recent studies focused on the ability of cells to deform and remodel 3D matrices using natural biomaterials [START_REF] Baker | Deconstructing the third dimension: how 3d culture microenvironments alter cellular cues[END_REF] . To measure traction forces exerted by cells at a given time, cells were first encapsulated in PEG containing fluorescent microbeads [START_REF] Legant | Measurement of mechanical tractions exerted by cells in three-dimensional matrices[END_REF] . In collagen gels containing fluorescent beads, the strain energy exerted by invasive or non invasive cells on their surroundings showed complex and anisotropic patterns [START_REF] Koch | 3d traction forces in cancer cell invasion[END_REF] . Using confocal reflection microscopy, local displacement fields [START_REF] Steinwachs | Three-dimensional force microscopy of cells in biopolymer networks[END_REF] measured in collagen networks surrounding breast carcinoma cells showed a surprising constant cell contractility regardless of the collagen matrix stiffness. Time-lapse imaging of migrating cells was also used to elucidate the molecular mechanisms underlying force transmission in 3D fibrin matrices [START_REF] Owen | A cytoskeletal clutch mediates cellular force transmission in a soft, three-dimensional extracellular matrix[END_REF] . Additionally, a few in vitro studies investigated the correlation between metastatic potential and migration of cancer cells in synthetic or physiological matrices. It was found that aggressive tumorigenic cells in 3D synthetic environments migrate differently as compared to normal primary cells [START_REF] Schwartz | A quantitative comparison of human ht-1080 fibrosarcoma cells and primary human dermal fibroblasts identifies a 3d migration mechanism with properties unique to the transformed phenotype[END_REF] . Finally, breast cancer biomarkers in association with morphological characteristics appear to be the dominant factor influencing cell motility in compliant collagen matrices [START_REF] Baker | Deconstructing the third dimension: how 3d culture microenvironments alter cellular cues[END_REF] .

Here, we addressed the role of cancer cell invasiveness on 3D motility, using three bladder cancer cell lines with different metastatic potentials [START_REF] Abidine | Mechanosensitivity of cancer cells in contact with soft substrates using afm[END_REF] in various collagen matrices, together with a recently developped method for calculating collagen displacements [START_REF] Fertin | Displacement fields using correlation methods as a tool to investigate cell migration in 3D collagen gels[END_REF] . We determined migration speeds at the lowest collagen concentration and found a significant difference, as well as a collagen concentration-dependent persistence. Then, using cell morphology and sphericity index, we showed that grade-3 cells were able to migrate using long protrusions and elongated morphology whereas grade-2 cells move with a round or slightly elongated shape showing short extensions. Then we used the Phase Correlation method [START_REF] Fertin | Displacement fields using correlation methods as a tool to investigate cell migration in 3D collagen gels[END_REF] to calculate the collagen displacement fields resulting from forces exerted by cells, together with the detailed analysis of actin. This analysis showed that collagen remodelling was closely connected to cell shape changes, associated with several actin-rich regions, playing a fundamental role in 3D migration. Furthermore, our results emphasized the various mechanisms (mesenchymal and amoeboid) used by cancer cells to migrate through these soft gels.

Results

Effects of metastatic potential and collagen concentration on cell motility

To determine how the metastatic potential of cancer cells is correlated with migration, we investigated the migration and effective speed of bladder cancer cells of increasing metastatic potential (RT112 < T24 < J82) in gels of different concentrations (0.95 mg/mL, 1.8 mg/mL, 4.5 mg/mL). Supplemental movie S1 gives an idea of the motion of a T24 cell migrating in the low concentration collagen (0.95 mg/mL). Cell migration speeds (i.e. average of instantaneous speeds) and effective speeds (total distance traveled divided by total migration time) at various collagen concentrations are shown in Figure 1A and 1C. The results reveal that the migration and effective speeds decrease with collagen concentration regardless of the metastatic potential. A strong decrease in migration speed is observed for J82 cells (65%), larger than in the case of T24 (50%) and RT112 cells (53%), as collagen concentration increases from 0.95 to 4.5 mg/mL, and similarly from 0.95 to 1.8 mg/mL (respective decrease of 56%, 35% and 38%). Regarding mean migration speed values, J82 cancer cells move with velocities in the range [7.3 -20.9 µm/h], as compared to [7.5 -14.9 µm/h] for T24 cells and [4.3 -9.1 µm/h] for RT112 cells. A strong decrease is also observed in terms of effective speed of J82 cells as collagen concentration increases from 0.95 to 4.5 mg/mL (80%), larger than in the case of T24 (63%) and RT112 cells (57%).

Interestingly, these results also reveal that grade-2 cancer cells (RT112) migrate slower than J82 cells at the lowest collagen concentration. Differences in cell speed are also visible at high collagen content: the migration and effective speeds of grade-2 cancer cells are statistically lower than T24 cells. On the other hand, no significant differences were measured between T24 and J82 cell speeds whereas these two cell lines do not have the same cancer stage. Note that J82 cells display a broad range of speed values at low collagen concentration, and this will be discussed later.

To further characterize the migration behavior of cancer cells at different collagen concentrations, we studied the persistence of cell migration using the MSD power-law exponent α and diffusivity D (see Materials and Methods section). Our data reveal that α and D are decreasing functions of the collagen concentration [START_REF] Laforgue | 3d cancer cell migration in collagen matrices[END_REF] . For the most invasive J82 cells (see Figure 1B), α ranges from super-diffusive values (mean α ∼ 1.39) at 0.95 mg/mL to random ones (mean α ∼ 1.09) at 1.8 mg/mL, and sub-diffusive (mean α ∼ 0.7) at 4.5 mg/mL. Similar variations in power-law exponents were observed for T24 and RT112 cells but were not significant. Moreover, the diffusivity decreases by a factor of about 10 for J82 cells of about 4 for T24 and RT112 cells between 0.95 mg/mL and 4.5 mg/mL (see Table 1). Finally, at low collagen concentration, we found that RT112 migration behavior is random (mean α ∼ 1.1) as compared to the super-diffusive motions of T24 and J82 cells (respective means α ∼ 1.37 and 1.39) (Figure 1B). A similar remark can be made in terms of mean directionality d(t) (see Materials and Methods) : RT112 cells show an important decrease of the directionality in time as compared to T24 and J82 cells (Figure 1D). Furthermore, these RT112 cell trajectories seem to reach limiting values around 0.2 whereas the T24 and J82 cell trajectories tend to 0.3-0.4. At this lowest collagen concentration, we also observed an effective speed of RT112 cells significantly lower than the J82 cell effective speed (Figure 1C) confirming a very different behavior of J82 cells.

Collagen concentration

Cancer cell morphology in 3D collagen matrices

To investigate the migration strategies used by the three cancer cell lines, we analysed the 3D cell morphology parameters (i.e. cell sphericity index, major radius) and introduced a new concept based on contour fluctuations. Figures 2A and2B show means and SEM of sphericity indexes as well as major radii for the three cell lines at various concentrations. At low concentration, grade-3 cells (T24 and J82 cells) were found to have similar small values of sphericity index, significantly smaller than grade-2 cells (0.426 for J82 cells, 0.429 for T24 cells and 0.467 for RT112 cells). At intermediate (1.8 mg/mL) and high concentrations (4.5 mg/mL), T24 cells showed a higher sphericity (0.495 and 0.514 respectively) as compared to RT112 cells (0.407 and 0.464) and J82 cells (0.402 and 0.432). J82 cells showed a major radius significantly higher (24.4 µm, 23,9 µm and 18.1 µm at increasing concentrations) than T24 and RT112 cells for all concentrations. To summarise, the general morphology of cancer cells migrating in the high concentration gels is different from the one found at low concentration:

• J82 cells display a significantly smaller major radius compared to low collagen concentration, with similar sphericity index (around 0.43).

• T24 cells also show a significant smaller major radius compared to low collagen concentration but with a rounder shape (sphericity index of 0.514 significantly higher than 0.429).

• RT112 cells display a small major radius, with the same sphericity (index around 0.46).

To investigate further the diversity in cell phenotypes (i.e. morphologies) at low collagen concentration, the contour fluctuations vs sphericity index were plotted for each time (Figure 2C). High values of the contour fluctuations are the signature of many cell protrusions (see insets in Figure 2C) whereas low values (< 20%) correspond to migrating cells exhibiting few protrusions. Consistent with these observations, the graph shows that the lower the sphericity index, the higher the shape irregularity. The same trend was found for medium and high collagen concentrations (data not shown). We observed that cancer cells exhibit a range of sphericity indexes between 0.3 and 0.6, corresponding to morphologies going from highly polarized cells with many extending protrusions (contour fluctuations > 40%) to very round cell shapes with few short extensions such as small membrane blebs (contour fluctuations < 20%). However, a certain fraction of cancer cells with low sphericity index (index < 0.4) exhibits low values of contour fluctuations (< 40%). In particular, metastatic T24 and J82 cells display a larger range of sphericity indices (0.3 to 0.6) with important shape differences (contour fluctuations from 5 to 60 %) whereas RT112 cells are restricted to a range of high sphericity indices (0.35 to 0.58) and smaller contour fluctuations (10 to 55 %). To further characterize the effect of invasiveness on cell morphology, the relationship between contour fluctuations and major radius was plotted in Figure 2D. This is an interesting criterion to determine cell invasiveness at low collagen concentration: J82 cells clearly exhibit the lowest fluctuation curve of all cell lines for a given size (Figure 2D) revealing that large fluctuations are 3/17 inefficient to enhance cell migration velocity (case of RT112 and T24 cells). This is also emphasized in Figure 2E where the migration speed is plotted as a function of contour fuctuations. Straight lines are used as a fit and reveal that the correlation (or slope) is larger for J82 cells.

To summarize, these parameters reveal the round morphology or slightly elongated shape of RT112 cells with a small major radius, an elongated morphology for T24 and J82 cells and the ability of grade-3 cells to show a great diversity of shapes at low collagen concentration.

Dynamic cell-induced 3D displacements

To study how cancer cells deform their environment, we determined the 3D displacements of collagen fibers induced by the migration of RT112, T24 and J82 cells. Figure 3A-C show the superposition of two successive reflectance confocal images obtained at 10 min time interval for each cell. The corresponding displacement fields around the cells are presented from two viewing angles (1rst angle in Figure 3D-F, 2nd angle in Figure 3G-I). The angle distributions of displacement vectors -relative to the direction of migration -are also shown in Figure 3J-L.

RT112 cells show a quite round morphology (Figure 3A and Figure 3D) at both times with no clear cell shape polarization. F-actin is located at the cell periphery in the form of aggregates. A short extension is visible at the leading edge. The cell moves without important shape changes. Despite the round cell shape, the displacement field spans a relatively large zone (see Figure 3D and Figure 3G) with a maximal value of 2.6 µm. Large fiber movements (> 1.4 µm) are oriented preferentially between 10 • and 65 • at the front and at the rear (see Figure 3J).

In Figure 3B, the images of the T24 cell correspond to an elongated morphology with a protrusion at the leading edge and a contraction of the rear (see white arrow on the overlay and cell shapes on the left panel). The cell is polarized as seen from its actin distribution. Its motion corresponds to a highly asymmetric 3D displacement field with a maximal value of 3.6 µm (see Figure 3E and Figure 3H). Collagen displacements are seen far from the cell (as far as 50 µm from the cell surface). The largest displacements (> 1.2 µm) are along the cell major axis and important displacements localized on the lateral side are also visible. Moreover, collagen fibers move in opposite directions between the leading edge and the rear edge: displacement vector angles (see materials and methods) lie in the range [10 • -100 • ] at the front as compared to [60 • -180 • ] at the rear (see Figure 3E and Figure 3K).

In the case of the J82 cell, the image overlay shows an elongated cell shape at both times (see Figure 3C). These images together with the 3D-displacement field reveal a long cylindrical protrusion (see white arrow on the left panels of Figure 3C and see Figure 3F and Figure 3I). Note that the protrusion displays several actin-rich regions at the tip, and also blebs at the surface. The magnitude of the fiber displacements reaches important values (up to 7.1 µm) with vectors showing [0-80 • ] orientation at the front and [0-180 • ] at the rear (see Figure 3L). The largest displacements (magnitude > 5 µm) are observed at the front along the migration axis and also at the rear with a mean orientation around 150 • (see Figure 3F and Figure 3I). At the rear part, some displacements are also visible at long distances with a mean orientation around 50 • .

Cell dynamics and collagen displacements can be summarized as follows:

-RT112 cells stay quite round with short extensions, while inducing extended fiber displacement fields with fibers moving in the direction of migration.

-The T24 cell is highly polarized with a protrusion inducing fiber movements towards the cell w.r.t. to the migration direction and a rear retraction generating fiber displacements directed away from the cell.

-The J82 cell motion induces collagen fiber movements in the direction of migration at the leading edge. Important and localized displacements are also visible at the rear.

These results are completed by a study of the collagen displacement dynamic in association with the cell behavior (see supplemental movies S2, S3 and S4). Movies show that T24 and J82 cells migrate with long or large protrusions. These cells show various displacement patterns and are able to push and pull on the collagen fibers. In contrast, the RT112 cell exhibits a slightly elongated shape with short extensions and only pulls on fibers located around protrusion/cell surface.

Correlation between fiber displacements and the F-actin distribution

In order to explore the role of metastatic potential on collagen deformations, we calculated the largest fiber displacements induced by the migration of bladder cancer cells for all times as well as the mean value of this maximal displacement for each cell line. This analysis is summarized in Table 2. At low collagen concentration (0.95 mg/mL), J82 cells -with the highest metastatic potential -generate significantly higher fiber displacements (mean value of 4.7 µm) as compared to T24 (2.7 µm) and RT112 cells (2.6 µm). Moreover, the magnitude of the largest displacements decreases significantly when collagen concentration increases for the three cell lines (see Table 2). Note that J82 cells show the highest value of displacement regardless of the collagen concentration.

Furthermore, to relate the largest fiber-displacements to the presence of high actin-density in the cytoskeleton, we compared the magnitudes and directions of displacement vectors in different regions. This analysis is shown in Figure 4 for two migration stages of a T24 cell at low collagen concentration (0.95mg/mL) (image #1 and image # 5 of supplemental movie S3). For the first stage, the cell shows a highly polarized morphology at times t and t+∆t (Figure 4A). We notice actin-rich regions at leading and rear edges for both times and actin speckles appear along the cell body (Figure 4A). The cell-induced displacement field during ∆t (= 10 min) exhibits large displacements located in regions 1 and 4 (∼ 2 µm in region 1 at the front, ∼ 2.8 µm in region 4 at the rear, Figure 4B). A study of the vector orientation with respect to the direction of migration showed a mean angle on the order of 42 • at the leading edge (region 1) whereas the movement of fibers at the rear (region 4) is opposed to migration (mean angle 141 • , Figure 4E). This is consistent with the displacement field (Figure 4B) showing vectors directed away from the actin zone at the front and back (Figure 4A-B). This pattern corresponds to the progression of the leading extension associated with the cell rear-contraction (see Figure 4A). Note that lateral sides (regions 2 and 3) correspond to high-displacement vectors (Figure 4B) but these are very localized and parallel to the cell surface. The mean angle in regions 2 and 3 is about 82 • and 63 • respectively. As anticipated, we note that the amplitude of the maximal displacements decreases as a function of cell-to-fiber distance for the actin-rich regions (regions 1 and 4 in Figure 4F), but this decrease is not observed for the other side areas (regions 2 and 3).

During the second stage (Figure 4C), the cell is polarized and exhibits two cylindrical protrusions. The actin-rich regions (aggregates/filaments, see dotted arrow in Figure 4C) are mainly localized in front of the nucleus and small actin spots are also visible at the end of protrusions. Surprisingly, in the nucleus periphery, we observe high displacements oriented towards the nucleus (regions 2 and 3 in Figure 4C-D). Indeed, displacement angles take values between 49 • and 65 • with respect to the direction of migration (Figure 4G). We also pay attention to the large movements of collagen fibers around the cylindrical protrusion and analyze their orientation (regions 1 in Figure 4G). In this region, the mean angle is 79 • . Thus, very localized and perpendicular fiber movements can be seen in the vicinity of protrusions. As before, close to the large actin-rich regions (2 and 3), we observe large displacements of collagen fibers with their amplitude decreasing with cell-fiber distance (Figure 4D). Taken together, these results clearly indicate that high fiber displacements are correlated with the presence of actin-rich regions within the cell. The displacement vectors are either directed away from or toward these large actin-rich areas. In the vicinity of these regions, displacement fields can be large and decrease as a function of cell-fiber distance.

This correlation between actin density and displacement fields can also be found in the collagen displacement movies (see supplemental movies S2, S3 and S4). It is enhanced for T24 and J82 in comparison with RT112 cells.

Finally, to prove the reliability of this previous arbitrary choice, we carried out a systematic analysis to obtain correlations between actin intensites and displacement magnitudes (or norms) for all cases. Figure 5 explains the procedure used. For example, a J82 cell is followed in time. Here, two time steps are shown (Figure 5A, t=20 min and Figure 5B, t=30 min) and Intensity Weighted Distances (IWDs, as defined in Mat. & Methods section) are visualized. IWDs are splitted into two sub-ranges corresponding to low and high values (respectively [0-25%] and [75-100%] ranges). This cell is elongated at t=20 min, then retracts at t=30 min. Therefore we expect strong displacements to be correlated with intense actin spots at time t=20 min. To check this, the time-evolution of the corresponding displacement norms (low and high domains) is described in Figure 5C-D (violin plot of displacement norm distributions and median values). The cell behaves differently in time. The corresponding actin spots can also be observed in Figure 5E. A selection of the more interesting times can be done using the criterion Norm high /Norm low >1. Doing this allows to present a cloud of data points (Figure 5F) for vector norm vs. actin intensity, for example here at t=20 min. The behavior, as shown by the linear regression (in red) gives a positive slope. All cases for the 11 cells studied and a total of 49 times gave positive slopes. These slopes are compared for the three cell types (J82, T24 and RT112) in Figure 5G and show differences, with a larger significant slope in the case of the most invasive J82 cells.

5/17

Discussion

The mechanisms of cell migration have been widely studied but were very often limited to 2D motion on various substrates simplifying cell-substrate interactions, cell morphology and obviously cytoskeleton dynamics. Recently, there has been considerable interest to investigate the cell motion in 3D matrices that better mimick the physiological environment. These results highlight the biophysical role of the matrix microarchitecture on cell behavior [START_REF] Carey | Biophysical control of invasive tumor cell behavior by extracellular matrix microarchitecture[END_REF] , the contribution of cell contractility, the role of fiber stiffness in 3D migration [START_REF] Doyle | Local 3d matrix microenvironment regulates cell migration through spatiotemporal dynamics of contractility-dependent adhesions[END_REF] , the role of actin bundles and integrins in force transmission to the surrounding environment [START_REF] Owen | A cytoskeletal clutch mediates cellular force transmission in a soft, three-dimensional extracellular matrix[END_REF] . Moreover, other studies emphasize a cell plasticity mechanism (mesenchymal-to-amoeboid transition) in cell migration favoring tumor dissemination and metastasis [START_REF] Pandya | Modes of invasion during tumour dissemination[END_REF][START_REF] Wolf | Compensation mechanism in tumor cell migration: mesenchymal-amoeboid transition after blocking of pericellular proteolysis[END_REF] . Cancer cells can use different modes of migration depending, among others, on the physical properties of the ECM [START_REF] Petrie | At the leading edge of three-dimensional cell migration[END_REF][START_REF] Parekh | Sensing and modulation of invadopodia across a wide range of rigidities[END_REF] . However, much of our knowledge about the migration in 3D environments comes from studies involving cells from different cancer types. Moreover, these studies do not combine simultaneously spatial and temporal analysis of cell behavior and cell-induced environment deformations.

In this study, we investigated how the malignancy of bladder cancer cells is related to their migration strategy in 3D collagen matrices mimicking physiological environments. We used three bladder cancer cell lines with different degrees of metastatic potential: RT112, T24 and J82 cells. The choice of these cells comes from previous studies, where various properties associated with their malignancy, such as adhesion and rheology, were reported [START_REF] Abidine | Mechanosensitivity of cancer cells in contact with soft substrates using afm[END_REF][START_REF] Laurent | Atomic force microscopy reveals a role for endothelial cell icam-1 expression in bladder cancer cell adherence[END_REF][START_REF] Abidine | Local mechanical properties of bladder cancer cells measured by afm as a signature of metastatic potential[END_REF] . RT112 cancer cells are moderately differentiated (grade 2), whereas T24 and J82 cancer cells are poorly differentiated (grade 3) and the J82 cell line has higher metastatic potential (see materials and methods). In addition, these cancer cells were used before as models for molecular classification and exhibit different mTOR-associated gene expression which is relevant for chemotherapeutic strategy [START_REF] Hau | Differential mtor pathway profiles in bladder cancer cell line subtypes to predict sensitivity to mtor inhibition[END_REF] . They express α5 and β 1 integrin subunits which are known to play a role in cell invasion by enhancing transmission and generation of contractile forces [START_REF] Knowles | Clt1 targets bladder cancer through integrin alpha5beta1 and clic3[END_REF][START_REF] Nista | beta 1-integrin expression and function in human bladder cancer cells: modulation by tnf alpha[END_REF] .

To avoid collagen degradation during cell migration, we limited experiments to a maximum time of 100 minutes, which restricts the MMP activation, therefore enzymatic degradation of matrices. To test if the cancer cell motility response to increasing matrix rigidity is malignancy dependent, we measured migration speeds at different collagen concentrations, 0.95, 1.8 and 4.5 mg/mL, corresponding to respective Young moduli E∼ 20, 250 and 1400 Pa [START_REF] Iordan | Breakdown of cell-collagen networks through collagen remodeling[END_REF] . We found values of migration velocities in agreement with recent data obtained in 3D rat tail collagen gels [START_REF] Doyle | Local 3d matrix microenvironment regulates cell migration through spatiotemporal dynamics of contractility-dependent adhesions[END_REF][START_REF] Wolf | Compensation mechanism in tumor cell migration: mesenchymal-amoeboid transition after blocking of pericellular proteolysis[END_REF][START_REF] Wu | Three-dimensional cell migration does not follow a random walk[END_REF] . It should be noted that Fraley et al. [START_REF] Fraley | Three-dimensional matrix fiber alignment modulates cell migration and mt1-mmp utility by spatially and temporally directing protrusion[END_REF] obtained higher values of speeds as they conduced cell motility studies up to 16.5h and they noticed that speed values decrease after MMP inhibitor treatment. As already noticed by Wolf et al. [START_REF] Wolf | Compensation mechanism in tumor cell migration: mesenchymal-amoeboid transition after blocking of pericellular proteolysis[END_REF] , we also observed that increasing the collagen concentration induced a decrease in the migration speed, regardless of the metastatic potential of bladder cancer cells (see Figure 1A). This result suggests that an increase in stiffness from 10 Pa to 1000 Pa has a significant effect on cell migration speed regardless of bladder cancer cell lines because of the limited space available due to smaller matrix pore size. Nevertheless the impact is more important for J82 cells: the relative decrease of J82 mean speed (65%) is more important than T24 (50%) and RT112 cells (53%) when collagen concentration increased from 0.95 to 4.5 mg/mL. A recent study demonstrated that the migration velocity of partially-or fully-transformed breast cells decreases when collagen density increases from 2 mg/mL to 4 mg/mL [START_REF] Baker | Deconstructing the third dimension: how 3d culture microenvironments alter cellular cues[END_REF] , this effect being more important for partially transformed cells. Our findings do not show the same trend as the migration speed decreases similarly for the three cell lines (∼ 20% from 1.5 to 4.5mg/mL). Possible explanations for this discrepancy could be related to cell type. Indeed, Baker et al. [START_REF] Baker | Deconstructing the third dimension: how 3d culture microenvironments alter cellular cues[END_REF] examined breast cancer cell progression series of lines established from a non cancerous human line whereas we used bladder cells lines established from human bladder cell carcinomas.

In addition, we also identified that collagen concentration plays a significant role for the movement persistence of J82 cells: the analysis of the mean square displacement (MSD) showed that the power law exponent (Figure 1B) significantly decreased when increasing collagen concentration for grade-3 J82 cells whereas no significant reduction was measured for the other cell lines. Moreover, this analysis highlighted a super-diffusive motion at low concentrations (α ∼ 1.39) and a sub-diffusive one at larger concentrations (α ∼ 0.7). These results were confirmed by higher values of diffusivities and mean directionalities for grade-3 cells (T24 and J82 cells) in comparison of RT112 cells (see Table 1 and Figure 1D). We concluded that grade-3 cell cells move more efficiently as compared to RT112 cells. Moreover, J82 cells exhibit the most effective migration. All these findings go along with morphological changes: an increase in collagen concentration (0.95 to 4.5 mg/mL) induces a decrease of the major radius of J82 cells (Figure 2B) and a reduction of contour fluctuations (30.7 ± 1.0 % to 24.0 ± 1.0 %) which is the signature of the number of extensions. In summary, the impact of collagen stiffness increase is more important for malignant J82 cells in terms of speed reduction, movement persistence and diffusivity. These results also suggest that matrix rigidity sensing is of utmost importance for malignant tumor cell motion, and that the response of grade-3 cells is characterized by changes of the major radius and contour fluctuations allowing cells to move more easily in less dense matrices.

Since J82 and T24 cancer cells migrate faster than RT112 cells in 0.95 mg/mL collagen, and the same for T24 vs. RT112 at 4.5 mg/mL, we compared the morphological parameters of the bladder cancer cell lines. At 0.95 mg/mL, grade-3 cells (T24 and J82) displayed mean sphericity indexes significantly lower than grade-2 cells (RT112). Moreover, by looking at the contour fluctuations as a function of sphericity index, we observed that grade-3 cells exhibited a larger range of sphericity indexes and contour fluctuations than grade-2 cells. Therefore grade-3 cells show a great diversity of shapes (from highly elongated to very 6/17 round ones) with variable numbers and types of protrusions (Figure 2C). This diversity of shapes and protrusions is coherent with the recent study of Eddy et al. [START_REF] Eddy | Morphodynamics facilitate cancer cells to navigate 3d extracellular matrix[END_REF] . However, J82 cells have a larger size as compared to other cells but exhibit the lowest shape fluctuations for a given major radius. These results suggest that a combination of low sphericity index and ability of cells to display a large range of contour fluctuations is optimal for migration in low collagen concentration. But the most interesting criteria to recognize the migration of a given cancer cell is by comparing its contour fluctuations at a given size (see Figure 2D) or look at the migration speed as a function of contour fluctuations (Figure 2E). These two graphs definitely exhibit different modes of migration used by invasive cells.

The less invasive cells (RT112 cells) showed a round or slightly elongated shape with short extensions or fine extensions or membrane blebs (see inserts in Figure 2C and movie S2). The analysis of the displacement field indicates that RT112 cells can move without important shape changes (movie S2) and no clear shape polarization. RT112 cells show sparse fiber movements directed in the migration direction along the cell surface (Figure 3D and images #2, #3 and #5 movie S2). Furthermore, we observed that high-density F-actin spots are visible at the cell periphery but not distributed isotropically (see Figure 3A and movie S2). This could signify that asymmetric actin accumulations at the cell periphery are needed to achieve directional migration and are consistent with previous studies about actomyosin contractility in amoeboid mode [START_REF] Chi | Rear actomyosin contractility-driven directional cell migration in three-dimensional matrices: a mechanochemical coupling mechanism[END_REF] . We also noticed that RT112 cells showing a relatively wide and short protrusion can generate fiber displacements around the edge of protrusions (image #1 movie S2). Notably, the cell produces scattered pulling movements on fibers. These results are consistent with an amoeboid cell behavior whose migration is independent of cell matrix adhesions [START_REF] Petrie | Activating the nuclear piston mechanism of 3d migration in tumor cells[END_REF] .

Previous studies reported that highly polarized cancer cells in 3D matrices induce large displacements near the protrusive tip [START_REF] Legant | Measurement of mechanical tractions exerted by cells in three-dimensional matrices[END_REF] and that forces are used almost exclusively to pull collagen fibers towards the cell [START_REF] Steinwachs | Three-dimensional force microscopy of cells in biopolymer networks[END_REF][START_REF] Hall | Toward single cell traction microscopy within 3d collagen matrices[END_REF][START_REF] Cóndor | Breast cancer cells adapt contractile forces to overcome steric hindrance[END_REF] . The results presented here showed that the migration of an elongated T24 cell can exhibit a series of various displacement patterns (see supplemental movie S3). First, protrusions push fibers at the front whereas the cell retracts at the rear. This creates relaxation of fibers away from the cell (Figure 3E and Figure 3H and first image of movie S3). In the next step, elongation of protrusions and adhesion/anchoring to collagen fibers results in a pulling force in the polarization direction as well as on the lateral sides (image #3 of movie S3). To our knowledge, these are the first time-dependent measurements of displacement fields showing that bladder cancer cells use a mesenchymal mode of migration thanks to protrusions, attachment, contraction and detachment, as in 2D [START_REF] Friedl | Tumour-cell invasion and migration: diversity and escape mechanisms[END_REF] . Obviously, since our study used 10 minute time lags, we observed the different mechanims at the same time. The formation of a 3D-lamellipodium at the leading edge generating a pushing force in the direction of migration and the relaxation of the collagen fibers at the rear (after the rear detachment) are visible in image #1 of movie S3. Interestingly, a clear difference appears between 2D and 3D migration: during the contractile phase, several actin fibers are localized in front of the nucleus, generating a pulling force. The amplitude of these contractile displacements decreases with cell-surface distance (Figure 4D). We also noticed the same decrease of the pushing displacement at the front while fibers relax at the rear, close to the large actin-rich structures (Figure 4F).

The migration data showed that invasive J82 cells are elongated and migrate using long cylindrical protrusions similar to a 3D-lamellipodium with high concentration of actin at the tip, with rapid membrane blebbing along the sides (Figure 3F, Figure 3I and supplemental movie S4) [START_REF] Petrie | Generation of compartmentalized pressure by a nuclear piston governs cell motility in a 3d matrix[END_REF] . J82 cells induce large displacement fields (maximum value of 7.2 µm) and the membrane is very dynamic (see Table 2 and supplemental movie S4). These results agree with recent studies describing mesenchymal-like phenotype of T24 and J82 cells [START_REF] Deb | Phosphoproteomic profiling identifies aberrant activation of integrin signaling in aggressive non-type bladder carcinoma[END_REF] . Deb et al. showed that the integrin pathway was significantly enriched and some major proteins regulating actin-polymerisation activity, cell motility and invasion are hyperphosphorylated for these two cell lines. Displacement fields for invasive J82 cells seem to be different from T24 cell: the cell pushes at the front whereas collagen fibers are pulled at the back of the nucleus to move the cell forward, while protrusions are growing (Figure 3F, Figure 3I and images #1, #2 and #5 in movie S4). Then, cell movement produces relaxation of those fibers located near the end of protrusions (where actin-rich areas are visible, images #2, #3 and #4 in movie S4). Note that the migration of the J82 cell gives rise to large displacements scattered in all directions (movie S4). These fiber displacements are the largest in comparison to T24 and RT112 cells, regardless of collagen content. This suggests that J82 cells generate stronger forces.

This was finally assessed by a systematic analysis of the correlations between actin intensites (as obtained with the Intensity Weighted Distance, the so-called IWD) and corresponding displacement magnitudes. As shown in Figure 5A-B, cells develop a variety of mechanisms. After selecting particular times when cells are polarized or exhibit a loss of symmetry (Figure 5C-D), we study the correlation between displacement norm and actin intensity at these times, and showed positive correlation (Figure 5F). Therefore larger displacements are always found in the regions close to the actin-rich cell domains. In addition, this correlation is significantly stronger (Figure 5G) for the more invasive J82 cells.

Conclusions

Cell migration in 3D matrices is a complex mechanism: the spatial deformation of the environment is not uniform and follows steps depending of cell phenotype. It is crucial to understand this spatio-temporal process in the case of cancer cells. Using 3D time-lapse imaging on a confocal microscope and a phase correlation method, we followed fibers deformations, as well 7/17 as cell morphology and live-actin distribution during cell movement. Our findings agree with a recent study describing the mesenchymal-like phenotype of T24 and J82 cells [START_REF] Deb | Phosphoproteomic profiling identifies aberrant activation of integrin signaling in aggressive non-type bladder carcinoma[END_REF] . In addition, our study revealed that grade-3 cells (T24 and J82 cells) were characterized by a great diversity of shapes in comparison to grade-2 cells (RT112 cells). Furthermore, we showed that the small shape fluctuations of highly invasive cells (J82) enable their migration in fibrous matrices with low collagen content. Our results also revealed that J82 cells produce the largest fiber displacements whatever the matrix rigidity. The detailed analysis of cell movements enhanced the fact that bladder cancer cells were also able to exhibit different phenotypes (mesenchymal and amoeboid) according to their metastatic potential. The different spatio-temporal mechanisms were emphasized and we showed that cells can push and/or pull collagen fibers in order to deform collagen efficiently. The actin-rich regions were found to be correlated with the largest displacement fields and the correlation is enhanced in the case of grade-3 cells.

Materials and Methods

Cell culture

Three epithelial bladder cancer cell lines, representing increasing malignancy states of progression, were used, RT112, T24 and J82 (ATCC, Rockville, MD). RT112 cells show an intermediate grade-2 and a stage Ta, i.e. these cells are moderately differentiated and have a small metastatic potential. T24 and J82 cells show a higher grade (grade-3) compared to RT112 cells and exhibit higher stages (T2-T3 and T3 respectively) i.e. they are poorly differentiated. J82 cells have the highest metastatic potential [START_REF] Abidine | Mechanosensitivity of cancer cells in contact with soft substrates using afm[END_REF][START_REF] Sundar Rajan | Unraveling the receptor-ligand interactions between bladder cancer cells and the endothelium using afm[END_REF] .

Cancer cells were cultured in RPMI 1640 complete medium (Gibco, Saint Aubin, France) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin. The cell lines were stably transfected with the LifeAct GFP plasmid to stain F-actin [START_REF] Riedl | Lifeact: a versatile marker to visualize f-actin[END_REF] . Cultures were grown at 37 • C in 5% CO 2 atmosphere. Cells were let to grow until 50% confluency was reached. Then they were detached from culture dishes using trypsin-EDTA (0.25%), were re-suspended at 2. 10 5 cells/mL concentration in culture medium, and finally seeded in collagen gels.

Collagen gel sample preparation

Gels were prepared from a rat tail collagen type I solution (Corning, New-York, USA) in Nunc LabTek II chambered coverglass 8 wells (Gibco, Saint Aubin, France). In order for the collagen gel to strongly adhere to the bottom glass surface, LabTeks were first treated with sodium hydroxide (NaOH 0.1N) for 1 min, and allowed to dry. Then 3-Aminopropyl-trimethoxysilane (APTMS) was added for 10 min. After this treatment, LabTeks were cleaned 3 times with distilled water and dried. Finally, glutaraldehyde (0.5% in PBS) was added for 30 min, and after the solution was removed, LabTeks were rinced 3 times with distilled water and allowed to dry before adding the collagen solution. All the following steps were then carried out in ice (around 4 • C) to prevent polymerization of the collagen before adding it into the LabTeks. Solutions containing collagen (0.95, 1.8 and 4.5 mg/mL), RPMI medium supplemented with glutamine, FBS, antibiotics and fibronectin (10 µg/mL) were prepared. After adding NaOH (0.1 M) to reach a pH of 7.4, freshly harvested cancer cells were included in the collagen mixture to obtain a final concentration of 2.10 6 cells/mL. The solution was poured into the LabTeks and transfered to an incubator for gel polymerization at 37 • C and 5% CO 2 for 30 min.

Live cell and collagen fibers imaging

The collagen fibers and the cell actin cytoskeleton were visualized simultaneously using the reflection and fluorescence imaging modes on a confocal microscope (Zeiss LSM710 model, Germany, equipped with a 40X water immersion objective, NA = 1.1). For this purpose, the 488 nm argon-laser wavelength was used to visualize the actin cytoskeleton of GFP-transfected cancer cells RT112, T24 and J82 [START_REF] Sundar Rajan | Unraveling the receptor-ligand interactions between bladder cancer cells and the endothelium using afm[END_REF] , in combination with a 633 nm HeNe Laser for the reflectance imaging of collagen fibers. The microscope set-up used a 80/20 beam splitter that directed 20% of the reflected light by the sample to the detector [START_REF] Iordan | Breakdown of cell-collagen networks through collagen remodeling[END_REF] . A specific chamber (37 • C, 5% CO 2 ) was mounted on the microscope stage in order to maintain cells in physiological conditions. Only isolated cells were investigated.

To construct a 3D image of migrating cells in a collagen matrix, we acquired Z-stacks of fluorescent/reflection images. The number of optical slices was different depending on cell acquisitions but the distance between two optical slices remained equal to 0.77 µm for all acquisitions. Image sizes were 512 x 512 x 128 pixels, with corresponding voxels 0.20 µm x 0.20 µm x 0.77 µm (or 0.32 µm x 0.32 µm x 0.77 µm for zooming purposes). To study the migration of cancer cells, the total acquisition time was around two hours and stacks were acquired every 10 minutes.

3D migration characterization

To quantify 3D migration [START_REF] Wu | Statistical analysis of cell migration in 3d using the anisotropic persistent random walk model[END_REF] , the velocity of cancer cells was calculated by tracking their center of geometry, using a home-made plugin under ImageJ. We used the green channel to locate the cell, and a threshold to determine its contour. The voxels inside the contour were kept, the volume was calculated and the center of geometry was determined precisely. This is the best way to obtain accurate data, as compared to other tracking methods such as bounding box center or ellipse fitting. Its position vector 8/17 was written as r(t), using the initial cell position as the origin. From the position vector, the instantaneous migration velocities were obtained and averaged in time [START_REF] Laforgue | 3d cancer cell migration in collagen matrices[END_REF] . Cells were classified as migrating cells if the center of geometry moved with a velocity higher than 3 µm/h 54 . The effective speed was calculated as the ratio of the distance traveled (from the initial to final positions of the cells) over total migration time. To evaluate the ability of the cancer cell to move in a complex Extra-Cellular Matrix (ECM), we calculated the mean square displacement (MSD) from the cell trajectories as follows:

MSD(t) = |r(t + τ) -r(τ)| 2 τ (1)
where τ is the lag time, and the brackets represent the sum of possible changes as τ varies from 0 to the longest time T. The MSD is usually characterized using a power-law relationship of the type MSD(t) = D(t/t 0 ) α where the exponent α represents the persistence of migration [START_REF] Mierke | The role of vinculin in the regulation of the mechanical properties of cells[END_REF] , D is the diffusivity and t 0 is the lag time. α = 1 corresponds to a ramdomly migrating cell while a value of 2 indicates a ballistically migrating cell. For the intermediate cases, 0 < α < 1 corresponds to a sub-diffusive motion whereas 1 < α < 2 corresponds to a super-diffusive migration. The MSD power-law exponent was computed from the cell trajectories recorded during 2 hours. When the MSD was not linear at large time lags, we computed the power-law exponent on a limited linear domain (fit on the first 60 minutes, to avoid the role of MMPs that could affect cell motion). Directionality d(t) was calculated from the ratio of the euclidean distance e d (t) between the starting and current positions divided by the real distance r d (t) travelled by the cell at time t : d(t)= e d (t) r d (t) .

Analysis of 3D-cell morphology

To investigate the morphology of migrating cancer cells in 3D-collagen matrices, we analysed their shape using a set of 3D-geometrical parameters. To capture this morphology with precision, we acquired confocal Z-stacks with 0.77 µm thickness slices. Then we applied the following procedure for all cells : a 3D gaussian blur (σ x = σ y = 1 in x, y and σ z = 2.6 in z), a median filter and thresholding using the OTSU algorithm. The cell major radius (R), obtained by fitting the cell contour by an ellipsoid, was defined as the length of the major axis of the ellipsoid. The sphericity ψ (ranging from 0 to 1) was defined as the ratio of the surface area of a sphere (with the same volume V p as the cell) to the cell surface A p .

ψ = π 1/3 (6V p ) 2/3 A p (2)
ψ equals 1 for a sphere and is close to 0 for a very elongated cell. Finally, we calculated the variation coefficient ν of the distance from the cell surface to its center, i.e. the ratio of the standard deviation (SD) to the mean distance (d m ): ν = SD d m . This ratio of contour fluctuations (in %) was evaluated in order to indicate the cell shape irregularity (ranging from 5% to about 60 % in this study).

3D displacement of collagen fibers

The 3D displacement of collagen fibers induced by the migrating cells was determined using a home-made phase correlation algorithm [START_REF] Fertin | Displacement fields using correlation methods as a tool to investigate cell migration in 3D collagen gels[END_REF] applied to the reflection image stacks obtained in time. The basic idea is to correlate only the positions of these fibers that are bright enough. Then we obtain incremental 3D displacements between times t and t + ∆t. The phase-only correlation (POC) was used and reads as:

r(x) = F -1 I * 1 (ω) * I 2 (ω) I * 1 (ω) * I 2 (ω) (3)
where i 1 (x) is the image intensity at time t, I 1 (ω) the Fourier transform of i 1 , I 1 (ω) = F (i 1 (x)), F -1 is the inverse Fourier transform, ω the angular frequency, x the position vector, and I * 1 (ω) denotes the complex conjugate of I 1 (ω). The same applies for i 2 and I 2 at time t + ∆t. The most noticeable property of the POC function is that it exhibits a unique sharp peak when the signal i 2 is a shifted version of i 1 . Using this property, we obtain the required displacement u in voxels, by locating the position of the peak of r(x):

u = arg max x (r(x)). (4)
The discrete displacement was calculated using two passes on a 512 x 512 x 128 pixels image. We started with a large window (64x64x16) located around the point of interest, and obtained a first estimate u u u of the displacement using phase correlation as described above, where the displacement is expressed in voxels. Then a smaller window (16x16x16) was used to estimate the extra displacement δ u u u using sub-pixel resolution as explained previously [START_REF] Foroosh | Extension of phase correlation to subpixel registration[END_REF] . The final displacement u u u sub is the sum of the two displacements u u u sub = u u u + δ u u u. Window parameters have been optimized and correspond to the ones given above. Note that de-noising was used by removing data with a small correlation coefficient, and other irrelevant displacement vectors showing very large differences with their neighbors [START_REF] Fertin | Displacement fields using correlation methods as a tool to investigate cell migration in 3D collagen gels[END_REF] .

9/17

Correlation between actin intensity and displacements

In order to capture precisely the relationship between actin intensities (or spots) with the displacements, we developed a systematic analysis to relate the displacement at position x with the closest actin region (within the cell at positions x c ). Different methods were used but the most satisfactory one consists in calculating the weighted actin density at the required displacement position x, i.e. the so-called IWD (Intensity Weighted Distance) :

IW D(x) = ∑ c I(x c )w(x -x c )/ ∑ w(x -x c ) (5)
where I(x c ) is the actin intensity at any voxel within the cell (the sum is on the cell voxels), and w(x -x c ) = 1/ xx c α is a weighting function based on the Euclidian distance between x and x c . α = 8 was found to be well adapted, indeed the weight should be stronger for the closest areas. This gives the actin weight at the point of application of the displacement.

To go further, we split the IWD regions (colored ones) into two major regions corresponding to lower [0-25%] and upper IWD quartiles [75%-100%] for simplicity. These two regions are named «high» (triangles) and «low» (circles). The corresponding magnitudes (norms) displacements have to be correlated to the values of IWD. They are plotted in time.

Statistics

Statistical analyses were performed using ANOVA. p < 0.05 was considered to be significantly different (*), ** corresponded to p < 0.01, *** corresponded to p < 0.001. All error bars were calculated using the standard error of the mean (SEM). When normality tests were not verified, statistical analyses were calculated using the Kruskal-Wallis test (morphological parameters and effective speed). p values < 0.05 were considered statistically significant.

Figure 1 .Figure 2 .Figure 3 .Figure 4 .

 1234 Figure1. Effect of cell metastatic potential and collagen concentration on (A) migration speed, (B) power law exponent, (C) cell effective speed, (D) Mean directionality (at 0.95 mg/mL). Three cancer cell lines of increasing invasiveness (RT112 < T24 < J82) and collagen gel concentration (0.95 mg/mL, 1.8 mg/mL and 4.5 mg/mL) were used. Power-law exponents α in (B) were obtained from MSD(t) = D (t/t 0 ) α . Box-whisker plots indicate median value, 25% and 75% quartiles and whiskers extend to the 5% and 95% percentiles (A, B and C). Statistical significance from another collagen concentration * corresponds to p < 0.05, ** corresponds to p < 0.01, *** corresponds to p < 0.001. Significant difference (+) from RT112 cells (+ p < 0.05, ++ p < 0.01). Cell number: (A) and (C) For 0.95 mg/mL, N= 9, 10, 7 respectively for J82, T24, RT112 cells. For 1.8 mg/mL, N= 8, 9, 7. For 4.5 mg/mL, N= 10, 13, 9. (B) For 0.95 mg/mL, N= 8, 9, 7 respectively for J82, T24, RT112 cells. For 1.8 mg/mL, N= 7, 8, 7. For 4.5 mg/mL, N= 8, 9, 9. (D) Graph representing the mean directionality vs. time for the three cancer cell lines at 0.95 mg/mL.

Figure 5 .

 5 Figure 5. Correlation between actin intensities and displacement magnitudes (J82 cell). A)-B) Values of IWDs are shown at two different times, 20 min, and 30 min. Symbols represent low and high IWDs, respectively [0-25%] and [75-100%] ranges. C) Violin plots of the displacement norm distributions corresponding to low (blue) and high (green) IWD values in time. D) Median of the norms (low and high IWDs) vs. time, with significant difference from Kruskal-Wallis test. E) Corresponding cell shapes. F) Displacement norms vs. weighted actin intensities (IWD) for t=20 min. The red line is the regression curve. G) Slopes of regression curves at times for which Norm high /Norm low >1, for the three cell types. Significance of Kruskal-Wallis test, followed by Dunn's test with Bonferroni adjustment.

Table 1 .

 1 Mean values of the diffusivity (in µm 2 ) ± standard error of the mean.

		Cell type	D (µm 2 )
	0.95 mg/mL	J82	11.4 ± 3.8
		T24	8.3 ± 3.4
		RT112	2.8 ± 0.7
	1.8 mg/mL	J82	3.4 ± 1.0
		T24	5.2 ± 1.7
		RT112	1.1 ± 0.2
	4.5 mg/mL	J82	1.3 ± 0.3
		T24	1.9 ± 0.7
		RT112	0.7 ± 0.1

Table 2 .

 2 Mean values of the largest fiber displacements (in µm) ± standard error of the mean. Statistical significance p (in comparison to J82 cells) for each collagen concentration is also indicated: * p < 0.05, ** p < 0.01, and *** p < 0.001.

	4/17

Acknowledgements

We thank the "Agence Nationale de la Recherche" (grant # 12-BS09-020-01, TRANSMIG). The authors are thankful to A. Duperray for fruitful discussions and for providing cancer cell lines. C.V. and V.M.L. are members of the LabeX Tec 21 (Investissements d'Avenir: grant agreement No. ANR-11-LABX-0030)

Author contributions statement

C.V. and V.M.L. designed research, wrote the manuscript. A.F., Y.U., C.V. and V.M.L. analysed the data. L.L. performed the experiments and computed migration speed, MSD, mean displacement. All authors reviewed the manuscript.

12/17


Samuel Valiquette

Frédéric Mortier

Jean Peyhardi

Gwladys Toulemonde

Choice of mixture Poisson models based on Extreme value theory

Keywords: Poisson mixture, Discrete extreme value theory, Peak-over-threshold

Count data are omnipresent in many applied fields, often with overdispersion due to an excess of zeroes or extreme values. With mixtures of Poisson distributions representing an elegant and appealing modelling strategy, we focus here on the challenging problem of identifying a suitable mixing distribution and study how extreme value theory can be used. We propose an original strategy to select the most appropriate candidate among three categories: Fréchet, Gumbel and pseudo-Gumbel. Such an approach is presented with the aid of a decision tree and evaluated with numerical simulations.

Introduction

Count data are classically observed in many applied fields such as in actuarial science when evaluating risk and the pricing of insurance contracts [e.g., [START_REF] Bartoszewicz | Modelling the claim count with Poisson regression and negative binomial regression[END_REF], in genetics to model the number of genes involved in phenotype variability [e.g., [START_REF] Anders | Differential expression analysis for sequence count data[END_REF] or in ecology to model species abundance [e.g., [START_REF] Wenger | Estimating species occurrence, abundance, and detection probability using zero-inflated distributions[END_REF]. While Poisson models and regression are well established choices for these type of data, they are not suitable for overdispersed data, which typically occur with an excess of zeroes or extreme large values. To overcome such limitations the use of Poisson mixture models has been proposed. This assumes the Poisson's intensity is no longer an unknown fixed value, but a positive random variable. Mixture approach induces overdispersed distributions with more zeroes and high values compared to the classical Poisson model [START_REF] Shaked | On mixtures from exponential families[END_REF]. A variety of mixture distributions has been already proposed [START_REF] Karlis | Mixed Poisson distributions[END_REF]. Classical examples include the gamma distribution [START_REF] Greenwood | An inquiry into the nature of frequency distributions representative of multiple happenings with particular reference to the occurrence of multiple attacks of disease or of repeated accidents[END_REF], the lognormal [START_REF] Bulmer | On fitting the Poisson lognormal distribution to species-abundance data[END_REF] or the Bernoulli [START_REF] Lambert | Zero-inflated Poisson regression, with an application to defects in manufacturing[END_REF]. From a general point of view, any distribution with a non negative support, finite or not, can be a potential candidate as a mixing distributions. Some tests exist to verify whether data are overdispersed [START_REF] Yang | A score test for overdispersion in Poisson regression based on the generalized Poisson-2 model[END_REF] or if they come from a mixed Poisson distribution [START_REF] Carriere | Nonparametric tests for mixed Poisson distributions[END_REF]. These tests justify the use of Poisson mixtures, but do not make claims on what type of mixing distribution should be selected. To our knowledge, there are no studies that propose a solution to this problem. This paper aims to propose a new strategy to select an appropriate and efficient mixing distribution family.

Usually one may choose to fit some predetermined Poisson mixtures and keep the best model based on some criteria. However the following example shows the difficulties in the mixing distribution choice. A sample of size n = 500 has been simulated from a Poisson-Beta type II with parameters equal to a = 1 and b = 2.2. This choice ensures finite expectation and variance for the mixed Poisson distribution and are equal to 0.83 and 8.47 respectively. Such a distribution has been used to model accident proneness [START_REF] Holla | On a discrete compound distribution[END_REF]. First and as expected, the use of a simple Poisson model fails to properly fit data (see Figure 1). In particular, it does not capture well the high frequency of zeroes and large values once its parameter is estimated ( λ = 1.02). Alternative latent distributions on λ have been used for this simulation, either the gamma or lognormal as classically proposed in ecology [START_REF] Tikhonov | Joint species distribution modelling with the R-package Hmsc[END_REF] or the inverse-gamma applied in actuarial science for liability insurance claims [START_REF] Tzougas | EM estimation for the Poisson-Inverse-Gamma regression model with varying dispersion: an application to insurance ratemaking[END_REF]. Inference is performed in a Bayesian framework (details in Section 2.2). While the Poissongamma, i.e. negative binomial, or the Poisson-lognormal are popular choices, the Poissoninverse-gamma is privileged for this example. Indeed, the posterior model probabilities of the negative binomial and Poisson-lognormal are respectively 0% and 38.8% compared to 61.2% for the Poisson-inverse-gamma. The latter mixture distinguishes itself even if the three models behave similarly with regards to Figure 1.

The reason that the inverse-gamma is privileged in this example is due to its tail behavior which is similar to the beta type II distribution. We establish in this paper that an adequate choice for the distribution on λ can be made by analysing such a property. This paper is organized as follow. Section 2 presents a classification of various Poisson mixtures based on their tail behavior using extreme value theory. Using these new classifications, we construct in Section 3 a strategy to choose a family of distributions for λ. This strategy is presented in the form of a decision tree where each step leads to an adequate category of Poisson mixtures. Simulations for each section are also presented to attest the relevance and usefulness of such an approach.

Poisson mixture tail behavior

In this section, we present the fundamental results in extreme value theory and the restrictions when it comes to discrete distributions. Following this, we present three categories of Poisson mixtures that characterise different tail behaviors and conclude with simulations to assess the interest of selecting such a category.

Theoretical foundations and Poisson mixtures categories

The tail behavior of a distribution can be studied using extreme value theory. Such a statistical approach analyses how the maximum of a distribution F stabilizes. The theory says that F belongs to a max domain of attraction if it exists two normalizing sequences a n > 0 and b n such that F n (a n x -b n ) converges to a non-degenerate distribution when n → ∞ [START_REF] Resnick | Extreme Values, Regular Variation and Point Processes[END_REF]. There are three possible domains of attraction named Weibull, Gumbel and Fréchet. These domains describe the asymptotic tail behavior of F and correspond respectively to finite, exponential and heavy tailed distributions. In the sequel, they will be denoted by D -, D 0 and D + , and we will write this property by F ∈ D where D is one of the three domains. In the following, we assume that the mixing distribution λ is supported on R + , then only D 0 and D + are considered. Usual continuous distributions belong to a domain, this is not the case for discrete random variables. Indeed, a necessary condition for a discrete distribution F to be in a domain of attraction is the property of long tailed [START_REF] Anderson | Extreme value theory for a class of discrete distributions with applications to some stochastic processes[END_REF] defined by:

lim n→∞ 1 -F (n + 1) 1 -F (n) = 1. (1)
In particular well known discrete distributions do not satisfy this property such as among others Poisson, negative binomial or yet geometric distributions. Even though the latter two are Poisson mixtures with a gamma distributed mixing parameter, which belongs to D 0 , the domain of attraction is not carried over by the mixture distribution. However [START_REF] Anderson | Extreme value theory for a class of discrete distributions with applications to some stochastic processes[END_REF] and [START_REF] Shimura | Discretization of distributions in the maximum domain of attraction[END_REF] showed that if a discrete distribution verifies

lim n→∞ 1 -F (n + 1) 1 -F (n) = L ∈ (0, 1) (2)
then F is, in a sense, 'close' to the Gumbel domain. More precisely, [START_REF] Shimura | Discretization of distributions in the maximum domain of attraction[END_REF] showed that property [START_REF] Stephenson | evd: Extreme value distributions[END_REF] implies that F is the discretization of a unique continuous distribution belonging to D 0 . On the other hand, [START_REF] Anderson | Extreme value theory for a class of discrete distributions with applications to some stochastic processes[END_REF] showed that there is a sequence b n such that F n (x + b n ) has infimum and supremum limits bounded by two different Gumbel distributions, implying that F is not far from this domain.

Because Poisson mixture distributions are discrete distributions, they are constrained to the long tailed property in order to have a domain. Otherwise, they may be close to the Gumbel domain. But Poisson mixtures are uniquely identifiable by the distribution on λ [START_REF] Feller | On a General Class of "Contagious" Distributions[END_REF], this means that its tail behavior is dependent of the latter. Therefore, we need to understand conditions on λ that allow Poisson mixture distribution to inherit a domain or to be close to the Gumbel one. [START_REF] Perline | Mixed Poisson distributions tail equivalent to their mixing distributions[END_REF] established some conditions for such preservation. From this point forward, we note by F and f the cumulative distribution function (cdf) and the probability density function (pdf) for λ and F M the cdf of the resulting Poisson mixture. Firstly, if F ∈ D + and is such that lim

x→∞ xf (x) 1-F (x) = α for some α > 0 (1 st Von Mises condition), then F M ∈ D + . Secondly, if F ∈ D 0 , lim x→∞ d dx 1-F (x) f (x) = 0 (3 rd Von Mises condition) and f (x) 1-F (x) = o(x -δ ) as x → ∞ for some δ ≥ 1 2 , then F M ∈ D 0 .
These results clarify some conditions for which domain of attraction of the mixing distribution is propagated to the associated Poisson mixture distribution. Very naturally, we denote these situations by two categories of Poisson mixtures: Fréchet and Gumbel. A broad set of distributions satisfies the 1 st Von Mises condition, for example the Fréchet, folded-Cauchy, beta type II, the inverse-gamma or the gamma/beta type II mixture [START_REF] Irwin | The generalized waring distribution applied to accident theory[END_REF]. Unfortunately, examples are scarce for the Gumbel domain. Indeed the extra condition on the hazard function is quite restrictive. Some examples are the lognormal, the Benktander type I and II [START_REF] Kleiber | Statistical Size Distributions in Economics and Actuarial Sciences[END_REF] or the Weibull distribution, with further restrictions on the parameters for the latter two cases. The Gumbel category does not encompass cases like the Poisson-gamma; while the mixing distribution belongs to D 0 , it does not satisfy the additional condition on the hazard function. In order to categorise such Poisson mixtures, we study distributions on λ that behave like the gamma.

Definition 1 [START_REF] Willmot | Asymptotic tail behaviour of Poisson mixtures with applications[END_REF]). A density f is a gamma type if

lim x→∞ f (x) C(x)x α e -βx = 1 (3)
where C(x) is a locally bounded function on R + and slowly varying, i.e. lim t→∞ C(tx) [START_REF] Bingham | Regular Variation[END_REF]), α ∈ R and β > 0.

C(t) = 1 for every x ∈ R + (see
Using gamma type distributions in the Poisson mixture context allows us to extend the categorisation to cases where F ∈ D 0 , but F M ∈ D 0 . First, we prove that such a mixing distribution belongs to D 0 (see Proposition 1). Finally Theorem 1 establishes that F M ∈ D 0 and quantifies the closeness to the Gumbel domain.

Proposition 1. If F is a gamma type distribution, then F ∈ D 0 .
Proof. Let F be the survival function of a gamma type distribution. A sufficient condition for F ∈ D 0 is to show that F has an exponential tail, i.e. lim x→∞ F (x + k)/F (x) = e -βk where k ∈ R [START_REF] Shimura | Discretization of distributions in the maximum domain of attraction[END_REF]. Using L'Hôpital's rule and equation (3), the limit becomes

lim x→∞ F (x + k) F (x) = e -βk lim x→∞ C(x + k) C(x) .
It remains to show that the latter limit is equal to 1. Because C(•) is slowly varying, we can use the Karamata representation

C(x) = c(x) exp x 1 t -1 η(t)dt ,
where c(•) and η(•) are both functions from R + to R + , lim x→∞ c(x) = c > 0 and lim x→∞ η(x) = 0 [START_REF] Bingham | Regular Variation[END_REF]. Then the limit equals

lim x→∞ C(x + k) C(x) = lim x→∞ exp x+k x t -1 η(t)dt .
Because η(x) → 0, for any ε > 0 then 0 < η(x) < ε for x large enough. Then

0 < x+k x t -1 η(t)dt < ε x+k x t -1 dt = ε log x + k x < ε log(1 + k)
which implies the limit is equal to 1 and establishes the sufficient condition.

Theorem 1. Let F M be a Poisson mixture with λ distributed according to a gamma type distribution F . Then for any integer k ≥ 1, lim n→∞

1-F M (n+k) 1-F M (n) = (1 + β) -k ∈ (0, 1). In particular, F M is not long tailed (k = 1).
Proof. Let P M and F M be the probability and survival functions of a Poisson mixture using a gamma type mixing distribution, then [START_REF] Willmot | Asymptotic tail behaviour of Poisson mixtures with applications[END_REF] showed that

lim n→∞ P M (n) C(n)n α (1 + β) -(n+α+1) = 1.
Using this result for integer k, we obtain

lim n→∞ F M (n + k + 1) -F M (n + k) F M (n + 1) -F M (n) = lim n→∞ P M (n + k + 1) P M (n + 1) = 1 1 + β k lim n→∞ C(n + k + 1) C(n + 1) ,
where the last limit converges to 1 using a similar proof as in Proposition 1. Because F M is monotonically decreasing, the proof can be concluded by applying the Stolz-Cesàro theorem.

This result allows to characterize a third category: the pseudo-Gumbel. It includes a broad class of mixing distributions among others gamma, gamma/Gompertz, exponential, exponential logarithmic, inverse-Gaussian and its generalization. Perline's result and Theorem 1 lead to consider now three categories for Poisson mixtures allowing the clarification of the mixing distribution choices. For examples, see Table 1 and the supplementary materials for details. Additionally, Theorem 1 also quantifies how 'close' those Poisson mixtures are with the quantity (1 + β) -1 involved in the limit. Indeed, if β → 0,

1-F M (n+1) 1-F M (n) → 1, i.e.
it approaches a long tailed distribution. Such property can blur the distinction between Gumbel and pseudo-Gumbel for some Poisson mixtures.

Impact of mixing distribution choice on goodness of fit

Are the categories previously defined useful to distinguish how the Poisson mixtures behaved? How can these categories be efficiently used when it comes to model selection? To answer these questions, we simulated 100 samples of different Poisson mixtures with size n = 250 using the gamma, Fréchet and lognormal distributions on λ, each one being a representative of the three categories. For each sample, the Poisson mixture is fitted with the same three distributions and the best model is kept using a Bayesian framework. This can be done with the rstan package for the language R [Stan Development Team, Mixing (λ) [START_REF] Meng | Simulating ratios of normalizing constants via a simple identity: A theoretical exploration[END_REF] and the dedicated R package Bridgesampling [START_REF] Gronau | Bridgesampling: An R package for estimating normalizing constants[END_REF]. All results are based on the following priors: a gamma(1, 1) distribution for positive parameters and a Normal(0, 1) for real parameters. Moreover, we simulated for each sample 4 MCMC with 10000 iterations each in order to ensure reasonable convergence for the parameter estimations and for the posterior model probabilities. Results are presented in Table 2 and, in general, the most selected model stood out. The only exception being the gamma [START_REF] Stephenson | evd: Extreme value distributions[END_REF][START_REF] Stephenson | evd: Extreme value distributions[END_REF] where the lognormal and the gamma Poisson mixtures are evenly selected throughout the simulations. This example shows the importance in the comparison of various mixing distributions. However, this approach based on systematic comparisons may suffer from computational limitations and is difficult to used in practice where a selection objective involves too much possibilities for the mixing distribution. Indeed, such an approach requires an appropriate choice of priors, a high number of MCMC iterations, and a study of convergence for each mixing distribution. Here, we assume the same priors and the same number of iterations for each mixing distribution in order to systematically simulate all our categories. In reality, each case must be studied with care and such an approach depends too much on the latent variables of the mixture model. As an alternative, we propose a simple strategy that uses directly the data and allows the user to focus on a specific family of mixing distributions. In doing so, the estimation of the latent variable can be done as a last step. The proposed alternative (see Section 3) relies on the use of a sequential approach by selecting first the most appropriate category (Gumbel, Fréchet or pseudo-Gumbel) and next by comparing only a few representative distributions belonging to this selected domain. For instance, such representative distributions are lognormal/Benktander-I for the Gumbel category, inverse-gamma/folded-Cauchy for the Fréchet category or a gamma/inverse-Gaussian for the pseudo-Gumbel case.

Poisson mixture (P M ) Category Fréchet(a, σ) Poisson-Fréchet Fréchet Folded-Cauchy(µ, σ) Poisson-folded-Cauchy Fréchet Inverse-gamma(a, b) Poisson-inverse-gamma Fréchet Beta-II(a, b) Poisson-beta-II Fréchet Gamma/Beta-II-mixture(r,a,b) Generalized Waring Fréchet Lognormal(µ, σ) Poisson-lognormal Gumbel Weibull(a, b) Poisson-Weibull Gumbel (if a < 0.5) Benktander-I(a, b) Poisson-Benktander-I Gumbel Benktander-II(a, b) Poisson-Benktander-II Gumbel b < 0.5) Exponential(a) Geometric Pseudo-Gumbel Gamma(a, b) Negative binomial Pseudo-Gumbel Inverse-Gaussian(µ, σ) Sichel Pseudo-Gumbel Generalized inverse-Gaussian(a, b, p) PGIG Pseudo-Gumbel

Strategy

This section proposes a strategy to choose a mixing distribution on λ using the categories defined in Section 2. As previously mentioned, an excess of zeroes and extreme values create overdispersion in count data which induces a particular tail behavior. The main idea is to choose mixing distributions among the three categories reflecting which ones best fit the empirical tail behavior. Peaks-over-threshold (POT) method [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF] is well adapted for this purpose. This technique analyses the distribution of the excesses defined by Y -u|Y > u. [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF], [START_REF] Balkema | Residual life time at great age[END_REF] showed that Y belongs to a domain of attraction if and only if the distribution of the excesses can be uniformly approached by a generalized Pareto distribution (GPD) as u tends to the right endpoint of the distribution of Y . The corresponding cdf is given by

H γ,σ (y) = 1 -1 + γ y σ -1/γ if γ = 0 1 -exp -y σ otherwise ( 4
)
with support R + if γ ≥ 0 or 0; -σ γ if γ < 0 and where γ ∈ R and σ > 0 are respectively shape and scale parameters. The sign of the γ parameter is intrinsically related to the domain of attraction. Indeed, the distribution of Y belongs to D -, D 0 or D + if γ < 0, γ = 0 or γ > 0 respectively. Therefore, fitting a GPD to the excesses of the count data can inform whether or not Poisson mixture distribution belongs to a known domain of attraction and, if so, which one.

Decision tree

Suppose overdispersed count data for which we need to fit a Poisson mixture. Our strategy to select appropriate mixing distributions is based on a decision tree (Figure 2) leading to the three categories defined in Section 2. The first step consists in the selection of a threshold u large enough for the data and to fit a GPD to the excesses. The choice of u can either be based on empirical quantile or on studying the mean residual life plot [START_REF] Scarrott | A review of extreme value threshold estimation and uncertainty quantification[END_REF] and the GPD parameters can be efficiently estimated using maximum likelihood [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF]. Two situations arise: first, if excesses are correctly fitted using the GPD model (left side in Figure 2), then we propose to use a distribution belonging to Gumbel or Fréchet such that the resulting Poisson mixture distribution remains in these domains. Indeed, if the excesses can be approximated by a GPD, then the Poisson mixture must be in a domain of attraction (see [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF], [START_REF] Balkema | Residual life time at great age[END_REF]). Therefore, the categories defined by [START_REF] Perline | Mixed Poisson distributions tail equivalent to their mixing distributions[END_REF] are an adequate set of mixing distributions to choose from. The choice of the specific category is directly based on the estimation of the shape parameter obtained at the first step. Testing whether γ = 0 or not can be done with the deviance statistic [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF]. If it is the case, we should use a Poisson mixture in the Gumbel category. Else if γ is positive, the Fréchet category should be prioritize.

Otherwise, if γ is negative, this strategy cannot assess which mixing distribution should be used. However, no such case has been relevant in our study.

The second situation (right side in Figure 2) corresponds to the case where GPD model is not well adapted. In this case, any mixing distributions such that the Poisson mixture belongs to Gumbel or Fréchet domains of attraction should be avoided. A distribution from the pseudo-Gumbel category should potentially be favoured. As demonstrated by [START_REF] Shimura | Discretization of distributions in the maximum domain of attraction[END_REF], such discrete random variables originate from an unique continuous distribution in D 0 that has been discretized. That is why we transform the excesses to continuous values thanks to a jittering technique consisting in the addition of a continuous random noise [START_REF] Nagler | A generic approach to nonparametric function estimation with mixed data[END_REF]. For an application, [START_REF] Coeurjolly | The median of a jittered Poisson distribution[END_REF] used this technique to study the Poisson's median and to construct an estimator for λ. In our case, the excesses Y -u|Y > u are discrete and greater or equal to 1. However the GPD with γ ≥ 0 is defined on (0, ∞). To transform the excesses into the same support, we jittered by subtracting a Uniform(0, 1). With these "jittered" data points, we fit a GPD again by fixing γ = 0 and testing again if it is adequate. If, in this case, the fit is adequate, we consider that the data are pseudo-Gumbel. Otherwise, one should proceed to another approach in order to choose a mixing distribution. However, we rarely encountered this situation in our simulations.

Evaluation of a sequential approach for mixing distribution selection

In order to study the performance of the proposed strategy, various Poisson mixtures samples have been simulated. The decision tree is then systematically applied using the evd package [START_REF] Stephenson | evd: Extreme value distributions[END_REF] for maximum likelihood estimation of GPD parameters, the modified Anderson Darling test for the goodness-of-fit and deviance statistic for testing the nullity of the shape parameter. Let X 1 , . . . , X m denote m i.i.d random variables ordered as

X (1) ≤ • • • ≤ X (m)
. The modified Anderson-Darling test statistic for a distribution H is defined by

T (X 1 , . . . , X m ) = n 2 -2 m i=1 H(X (i) ) - m i=1 2 - 2i -1 n log(1 -H(X (i) )).
As presented in [START_REF] Chu | A review of goodness of fit tests for Pareto distributions[END_REF], this statistic has an asymptotic distribution defined by a weighted sum of χ 2 1 . For this simulation study, the m random variables are the excesses and H is the GPD. We point that such a test works for any distribution H. However, some tests exist specifically when H is a GPD. See [START_REF] Falk | A LAN based neyman smooth test for pareto distributions[END_REF] or Villaseñor-Alva and González-Estrada [2009] for examples.

Finally, to test H 0 : γ = 0 versus H 1 : γ = 0, we fit the two models, that is the complete one and the restricted one, evaluate the corresponding log likelihoods namely L 1 and L 0 and we conclude with the deviance statistic D = 2 (L 1 -L 0 ) which follows a χ 2 1 under suitable conditions [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF].

The following simulation scheme is then applied in the language R:

1. For a fixed sample size n and a Poisson mixture F M with fixed parameters, simulate the mixed Poisson observations

Y := (Y 1 , . . . , Y n ).
2. For a threshold u based on the sample Y (example: 95th quantile), get the excesses X := Y -u|Y > u.

3. Calculate the MLE of γ and σ of the GPD for X using evd::fpot function using the Nelder-Mead optimization method.

4. Test the GPD for X with the modified Anderson Darling test (α = 0.05). The pvalues are calculated with a bootstrap approach using 250 iterations (see [START_REF] Chu | A review of goodness of fit tests for Pareto distributions[END_REF]).

5. Evaluate which category the sample is classified to using the decision tree with the following outcomes:

(a) If the test at step 4 for X is not rejected, use the deviance statistic to find which domain of attraction X belongs to. If γ < 0 is significant, the sample fails to have a category.

(b) Else repeat steps 3 and 4 for the jittered excesses X c := X -Unif(0, 1) and by fixing γ = 0. If the GPD is not rejected, the sample belongs to pseudo-Gumbel category. Otherwise the sample fails to have a category.

6. Repeat 1000 times the steps 1 to 5.

Distributions from the three categories are tested using these steps with n equal to 1000 or 2000 and the threshold fixed at the 95th and 97.5th empirical quantiles for n = 1000 and we add the 98.5th empirical quantile for n = 2000. For the Fréchet category, Fréchet and folded-Cauchy mixing distributions are simulated. For the Gumbel category, lognormal and Weibull mixing distributions have been simulated. Finally, gamma and inverse-Gaussian mixing distributions are simulated for the pseudo-Gumbel. Results are presented in Table 3 and also in the Supplementary Materials section (other sets of parameters, different sample sizes and different thresholds). Table 3: Average number of excesses, sample proportion of GPD rejection and of the most frequent category for the simulations with n = 1000 and u = 95th empirical quantile

These simulations aim to assess whether the decision tree adequately identifies the Poisson mixture categories. To do so, we calculate the proportion of samples where the GPD is rejected in the first branch and the proportion of the most frequent category where the Poisson mixture is identified. For cases like the Gumbel and Fréchet categories, we should see a low GPD rejection frequency. Conversely, the pseudo-Gumbel category should have a high GPD rejection. For all cases, we should have a high proportion of samples adequately identified to their category.

As presented in Table 3, most of the Poisson mixtures simulated can be adequately identified to the appropriate category. For the mixing distributions Fréchet(1,1), folded-Cauchy(0,1), lognormal(1,1) and Weibull(0.5,1), the GPD is mostly adequate for the excesses due to its low rejection rate. For pseudo-Gumbel mixtures Gamma(2,1) and inverse-Gaussian(1,2), the GPD is mostly rejected and, once the excesses are jittered, the Gumbel domain is found. Moreover, the simulation results considering a Weibull(a,b) as a mixing distribution are relevant with the theory. Indeed, as mentioned in Table 1, the Poisson-Weibull is in the Gumbel category if its parameter a is smaller than 0.5. In comparison to the limit case Weibull(0.5, 1) in Table 3, we also tested the Weibull(1,1) and found that this mixture does not belong to Gumbel category (Table 2 Supp. Material). Some interesting factors have been identified concerning the categorization. First, as noted by [START_REF] Hitz | Discrete extremes[END_REF], the discrete excesses need a certain amount of variability in order to have a smooth adjustment to the GPD. If the variance is not high enough, it can be difficult in practice to adequately identify the mixture category. For example, based on the simulations, the Fréchet(2,1) Poisson mixture belongs to the Fréchet domain only after the excesses are jittered (Table 1 Supp. Material). This distribution has its expectation defined which results in a more stable sample compared to the Fréchet(1,1). Here, both distributions do not have a defined variance because their parameter a ≤ 2. Still, because the Fréchet(2,1) is the limit case, it leads to a less volatile sample of excesses. This situation is identified in the lognormal case as well. Indeed, the Poisson mixture using the lognormal(0,1) has a smaller variance compare to mixture with a lognormal(1,1). In our simulations, the former case has some difficulty to be adequately identified to the Gumbel category compared to the latter (Table 2 Supp. Material).

Another important aspect is the choice of the threshold u. Indeed, the threshold affects the variance and the domain of attraction inferred by the data. For example, the lognormal(0,1) may have difficulties to be identified due to its variance, but once the u is large enough we do find the Gumbel category. For instance when n = 1000 and u = 95th quantile, 257 over 1000 samples are classified to the Gumbel category compared to 476 over 1000 samples in the pseudo-Gumbel category. However, for the same n and u = 97.5th, 790 over 1000 samples are classified to the Gumbel category compared to 64 over 1000 samples in the pseudo-Gumbel category. Therefore a larger threshold u is necessary in this case. However, it can be too large for some distributions like the gamma. In particular, when n = 1000 and u = 95th quantile there should be in average 50 excesses, but Table 3 indicates that the gamma(2,1) has 36.2 excesses. This can be explained by the lack of different values in the right tail of the Poisson mixture which leads to an underrepresented sample of excesses. Such discrepancy greatly affect the categorisation. For example, with n = 2000 and u = 98.5th quantile, the gamma(2,1) has in average 21.927 excesses and are mostly classified in the Gumbel domain (Table 3 Supp. Material). Clearly this inappropriate classification is due to the lack of excesses.

Finally, some pseudo-Gumbel Poisson mixtures can be very close to the Gumbel domain. For instance, the inverse-Gaussian(1, 2) is adequately classified, but the inverse-Gaussian(2, 1) is mostly classified in the Gumbel category (Table 3 Supp. Material). This can be explained by the 'closeness' property described in Theorem 1. Indeed, the density of a inverse-Gaussian(µ, σ) can be represented by

f (x) = C(x)x -3/2 exp - σ 2µ 2 x
where C(x) = C exp -σ 2x with C the normalizing constant. By equation (3), β = σ/2µ 2 and substituting the values µ = 2 and σ = 1 gives β = 1/8. Therefore the limit defined in Theorem 1 indicates that the resulting Poisson mixture has

lim n→∞ 1 -F M (n + 1) 1 -F M (n) = 8 9 .
Because this limit is pretty close to that of a long tailed distribution, i.e. the limit is near to 1, the distinction between pseudo-Gumbel and Gumbel gets blurred. To visualize how this 'closeness' affects the fit, let the parameter µ = 2 be fixed and vary σ from 0.1 to 8 for the inverse-Gaussian. For each value of σ, simulate 500 samples of size n = 2000 from the Poisson mixture, fix the threshold u to the 97.5th empirical quantile and calculate the proportion of samples where the GPD is rejected with p-value α = 0.05. In theory, the limit from Theorem 1 should approach 0 when σ gets larger, which implies getting further from the Gumbel domain and results into more rejection of the GPD. This result is reflected in Figure 3 when σ approaches 8. where the GPD has been rejected (α = 0.05) for the excesses (u = 97.5th quantile) according to σ.

Conclusion and perspectives

Overdispersed count data are commonly observed in many applied fields and Poisson mixtures are appealing to model such data [START_REF] Karlis | Mixed Poisson distributions[END_REF]. However, the choice of the appropriate mixing distribution is a difficult task relying mainly on empirical approaches related to modelers subjectivity or on intensive computational techniques combined with goodness-of-fit test or information criteria. In this paper, we proposed a new strategy based on the analysis of the tail behavior of the data. We extend the usual Gumbel and Fréchet domains of attraction introducing the pseudo-Gumbel category for Poisson count data. In particular, we show how tail behavior can provide a great amount of information to evaluate the mixing distributions. Based on a sequential strategy and decision tree, we proposed a useful and efficient approach to select the most appropriate category allowing to focus on a more restrictive set of potential candidates. The choice of the most appropriate distribution within a given category is not dealt with in this paper. Some strategies can be proposed helping the choice of such potential candidates. More specifically, it could be based on the simplicity either for the inferential step or for the inclusion of covariates or yet for biological interpretations. Moreover, recently, tremendous researches have been developed to jointly model count data. For instance, the joint species distribution models are proposed extending classical species distribution models in ecology [START_REF] Tikhonov | Joint species distribution modelling with the R-package Hmsc[END_REF] and are often based on the use of the multivariate lognormal distribution [START_REF] Aitchison | The multivariate Poisson-log normal distribution[END_REF]Ho, 1989, Chiquet et al., 2021]. Based on our approach, various and flexible models could be developed combining different mixing distribution belonging to different categories (Gumbel, Fréchet or pseudo-Gumbel) and the use of copulas to model dependencies structures between continuous mixing distributions [START_REF] Nelsen | An Introduction to Copulas[END_REF].

Figure 1 :

 1 Figure 1: Fitted Poisson, Poisson-Gamma, Poisson-Inverse-Gamma and Poisson-Lognormal models to data simulated from a Poisson-Beta type II distribution with a = 1 and b = 2.2 compared to its empirical distribution (Dark Gray)

Figure 2 :

 2 Figure 2: Decision tree for Poisson mixtures

Figure 3 :

 3 Figure 3: Proportion of inverse-Gaussian(2, σ) Poisson mixture samples (size n = 2000)where the GPD has been rejected (α = 0.05) for the excesses (u = 97.5th quantile) according to σ.

Table 1 :

 1 Examples of Poisson mixture and associated categories 2020] to estimate the parameters by MCMC. The best model is then kept using the highest posterior model probability. Those probabilities are approximated using the bridge sampling computational technique

Table 2 :

 2 Selected model frequencies for each Poisson mixture simulation with the highest frequency in bold.

	Poisson mixture Fréchet Lognormal Gamma
	Fréchet(1,1)	95	5	0
	Fréchet(2,1)	77	17	6
	Lognormal(0,1)	16	63	21
	Lognormal(1,1)	6	86	8
	Gamma(2,1)	0	23	77
	Gamma(2,2)	1	47	52

Acknowledgments

This research was supported by the GAMBAS project funded by the French National Research Agency (ANR-18-CE02-0025) and the french national programme LEFE/INSU. We also thank Éric Marchand for the helpful comments and fruitful discussions.