Daria S Roshal

Marianne Martin

Kirill Fedorenko

Ivan Golushko

Virginie Molle

Stephen Baghdiguian

Sergei B Rochal

Random nature of epithelial cancer cell monolayers

Although the polygonal shape of epithelial cells has drawn the attention of scientists for several centuries, only a decade and a half ago, it has been demonstrated that distributions of polygon types (DOPTs) are similar in proliferative epithelia of many different plant and animal species. In this study we show that hyper-proliferation of cancer cells disrupts this universal paradigm and results in randomly organized epithelial structures. Examining non-synchronized and synchronized HeLa cervix cells, we suppose that the spread of cell sizes is the main parameter controlling the DOPT in the cancer cell monolayers. To test this hypothesis, we develop a theory of morphologically similar random polygonal packings. By analyzing differences between tumoral and normal epithelial cell monolayers, we conclude that the latter have more ordered structures because of the lower proliferation rates and, consequently, more effective relaxation of mechanical stress associated with cell division and growth. To explain the structural features of normal proliferative epithelium, we take into account the spread of cell sizes in the monolayer. The proposed theory also rationalizes some specific highly ordered post-mitotic unconventional epithelia.

Introduction

Symmetry and topology determine the structure and laws of motion for relatively simple abiotic systems studied by physics and chemistry. In living systems, gene expression is usually considered as the fundamental mechanism controlling development and homeostasis [START_REF] García-Bellido | Genetic control of wing disc development in Drosophila[END_REF][START_REF] Silva | The biological chemistry of the elements: the inorganic chemistry of life[END_REF]. Nevertheless, the polygonal (prismatic in 3D) shape of cells, and their highly ordered packing in epithelia, clearly demonstrate that the organization of these cell monolayers directly follows basic physical and topological rules [START_REF] Azzag | The non-proliferative nature of ascidian folliculogenesis as a model of highly ordered cellular topology distinct from proliferative epithelia[END_REF][START_REF] Farhadifar | The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing[END_REF][START_REF] Nagpal | Epithelial topology[END_REF]. Epithelial growth is achieved by intercalary cell divisions within a constrained volume [START_REF] Rauzi | Cell intercalation in a simple epithelium[END_REF]. Dividing cells change their mediolateral neighbors, thus maintaining the apico-basal architecture and tightness of the intact layer [START_REF] Gibson | Cell Topology, Geometry, and Morphogenesis in Proliferating Epithelia[END_REF], which retains both robustness and plasticity. The most striking phenomenon is the so-called topological invariance observed in almost all proliferative epithelia within phylogenetically distant organisms harboring different global architectures. In fact, during their formation, different epithelial structures converge to polygonal packings with very similar DOPTs [START_REF] Lewis | The effect of cell division on the shape and size of hexagonal cells[END_REF][START_REF] Korn | The geometry of plant epidermal cells[END_REF][START_REF] Kwiatkowska | Growth and morphogenesis at the vegetative shoot apex of Anagallis arvensis[END_REF][START_REF] Gibson | The emergence of geometric order in proliferating metazoan epithelia[END_REF][START_REF] Patel | Modeling and inferring cleavage patterns in proliferating epithelia[END_REF]. Therefore, in very different epithelia the probabilities to observe cells with the same number of sides are approximately equal. This topological invariance seems to be closely associated with the physiological invariance of epithelium: in all eumetazoans, epithelial structures form a selective paracellular barrier, which controls fluxes of nutrients, regulates ion and water movements, and limits host contact with antigens and microbes. This universal function of epithelia is achieved by the maintenance of the epithelial tightness throughout various morphogenetic processes, like embryonic development, organogenesis, or continuous cell renewal [START_REF] Axelrod | Cell shape in proliferating epithelia: a multifaceted problem[END_REF][START_REF] Chelakkot | Mechanisms regulating intestinal barrier integrity and its pathological implications[END_REF][START_REF] Higashi | Maintenance of the epithelial barrier and remodeling of cell-cell junctions during cytokinesis[END_REF].

Topology of the cellular borders was first studied back in 1928 [START_REF] Lewis | The effect of cell division on the shape and size of hexagonal cells[END_REF]. In this pioneer work, the DOPT in the proliferative epithelium of cucumber was investigated: out of 1000 epithelial cells studied, 474 were hexagonal, 251 were pentagonal, and 224 were heptagonal. Later, in an excellent article [START_REF] Gibson | The emergence of geometric order in proliferating metazoan epithelia[END_REF], the hypothesis of topological invariance was formulated. It was established that the similar distributions of polygons are observed in epithelia of several more animal species, and a Markov-type theory was proposed to explain such invariance. This theory assumed that the junction that appears during cell division is always located in such a manner that it forms two additional polygonal vertices, which are necessarily on non-adjacent sides of the original polygon. It was also supposed that the probability of the junction formation is independent of the ratio between the parts in which the cell is divided. Nevertheless, the theory [START_REF] Gibson | The emergence of geometric order in proliferating metazoan epithelia[END_REF][START_REF] Patel | Modeling and inferring cleavage patterns in proliferating epithelia[END_REF] cannot explain the existence of 4-valent cells, for which the observed fraction varies from 2 to 3% [START_REF] Lewis | The effect of cell division on the shape and size of hexagonal cells[END_REF][START_REF] Gibson | The emergence of geometric order in proliferating metazoan epithelia[END_REF]. In addition, this approach considers neither the possibility of cellular motility nor mechanical interactions between cells, which are crucial for the the epithelium properties [START_REF] Puliafito | Collective and single cell behavior in epithelial contact inhibition[END_REF][START_REF] Angelini | Glass-like dynamics of collective cell migration[END_REF][START_REF] Alt | Vertex models: from cell mechanics to tissue morphogenesis[END_REF].

These facts together with new findings regarding the effect of mechanical stress on the reorganization of focal adhesion, adherent junctions, and cytoskeleton, as well as cell division and growth [START_REF] Kwiatkowska | Growth and morphogenesis at the vegetative shoot apex of Anagallis arvensis[END_REF][START_REF] Puliafito | Collective and single cell behavior in epithelial contact inhibition[END_REF][START_REF] Ingber | Cellular mechanotransduction: putting all the pieces together again[END_REF][START_REF] Yeung | Effects of substrate stiffness on cell morphology, cytoskeletal structure, and adhesion[END_REF][START_REF] Janmey | Cell mechanics: Integrating cell responses to mechanical stimuli[END_REF] motivated development of the microscopic models of epithelia. In the last 15 years, several approaches, which consider cells as polygons and define energy of the system as a function of their areas and lengths of their edges, were proposed to explain the influence of mechanics on the topology and collective behavior in cellular monolayers [START_REF] Farhadifar | The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing[END_REF][START_REF] Fletcher | Vertex models of epithelial morphogenesis[END_REF][START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF][START_REF] Bi | A density-independent rigidity transition in biological tissues[END_REF][START_REF] Staple | Mechanics and remodelling of cell packings in epithelia[END_REF][START_REF] Giavazzi | Flocking transitions in confluent tissues[END_REF][START_REF] Li | Mechanisms of regulating cell topology in proliferating epithelia: Impact of division plane, mechanical forces, and cell memory[END_REF]. One of the major successes of these models is a discovery of the solid-liquid transition controlled by the effective parameter, so-called target shape index [START_REF] Farhadifar | The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing[END_REF][START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF][START_REF] Bi | A density-independent rigidity transition in biological tissues[END_REF] that is defined by the competition between active contractility of actin-myosin subcellular cortex, cortical tension, and cell-cell adhesion [START_REF] Bi | A density-independent rigidity transition in biological tissues[END_REF]. It has been demonstrated that decrease of the adhesion between cells leads to the "jamming" of the soft liquid-like phase which is accompanied by the increase of the share of sixvalent cells and overall ordering of the structure. Experimental observations confirm that epithelial tissues indeed exhibit glassy behavior essential for such processes as embryonic development, cancer metastasis, and wound healing [START_REF] Bi | A density-independent rigidity transition in biological tissues[END_REF]. Normal development of epithelium involves mesenchymal-to-epithelial transition (MET) leading to the dramatic decrease of cell mobility and mitotic rate. At this transition, cells take more regular shapes with smaller average perimeter. The reverse process, leading to the appearance of elongated cells with high mobility, is called the epithelial-to-mesenchymal transition (EMT) [START_REF] Chaffer | Mesenchymal to epithelial transition in development and disease[END_REF]; cells undergoing oncogenic EMT are believed to drive metastasis process [START_REF] Thompson | Carcinoma Invasion and Metastasis: A Role for Epithelial-Mesenchymal Transition?[END_REF]. In this context, new data on relatively easily measurable geometric parameters of cell monolayers as DOPTs and distribution of cell areas can be helpful to distinguish pathological states and develop theoretical models of cancer and normal epithelia.

In this paper, we investigate the structural characteristics of non-synchronized (conventional) and synchronized confluent monolayers obtained from HeLa epithelial cancer cells. The latter monolayers are of particular interest, since most of the cells constituting this epithelial model belong to the second generation, and the Markov-type theory used in ref. [START_REF] Gibson | The emergence of geometric order in proliferating metazoan epithelia[END_REF] goes beyond the scope of its applicability. Moreover, despite numerous previously published data [START_REF] Gómez-Gálvez | The complex three-dimensional organization of epithelial tissues[END_REF], a putative influence of the cell cycle duration and cell synchronization on the epithelial structure remains unclear [START_REF] Axelrod | Cell shape in proliferating epithelia: a multifaceted problem[END_REF].

In order to rationalize the observed epithelial structures, we propose a new geometrical model generating polygonal packings that are very similar to the structures observed in non-synchronized and synchronized HeLa monolayers. Testing and applying our approach to several normal proliferative epithelia, we propose the physical mechanism underlying the topological difference between normal and hyperproliferative epithelia. As we demonstrate, in the epithelia with lower proliferation rates the relaxation of mechanical stresses associated with cell division and cell growth results in more ordered structures and maintains the topological invariance.

Results

Structural characterization of cancer cell monolayers.

Before carrying out the structural characterization of the epithelial monolayers that we obtained, it is important to discuss some general properties of DOPT. Note that the cell polygons in the epithelium are generally convex and form a tessellation of the monolayer surface that is similar to the Voronoi one [START_REF] Delaunay | Sur la sphère vide: A la mémoire de Georges Voronoi (French)[END_REF], which, in turn, is a dual of Delaunay triangulation [START_REF] Okabe | Spatial tessellations: concepts and applications of voronoi diagrams[END_REF]. In this case, for an infinite arbitrary flat monolayer, the average number of nearest neighbors equals 6. Indeed, if the Gaussian curvature is absent, then the equality Δ=0 takes place, where

Δ = ∑ 𝑃 𝑖 • (𝑖 -6) ! ,
(1) 𝑃 " = 𝑁 " /Σ # 𝑁 # is the concentration of cells with i nearest neighbors and 𝑁 # is the total number of cells with j nearest neighbors. Recall that the quantity 𝑄 = ∑ 𝑁 " • (𝑖 -6) " , proportional to Δ, is called the topological charge [START_REF] Conway | Sphere packings, lattices and groups[END_REF][START_REF] Roshal | Crystal-like order and defects in metazoan epithelia with spherical geometry[END_REF][START_REF] Roshal | Extended topological defects as sources and outlets of dislocations in spherical hexagonal crystals[END_REF]. For a triangulation of the sphere, Q=12, for a torus, as well as for an infinite plane, Q=0.

In fact, the non-local equality Δ=0 or the equivalent statement about the average number of neighbors means that the DOPT must be balanced: the number of cellular n-gons, with n<6, must balance the number of n-gons with n>6. In this equilibrium, the weights of 5-and 7-gons are equal to one, the weights of 4-and 8-gons are twice as large, etc. equation ( 1) can be used to estimate the error in the experimental determination of DOPT. For a finite monolayer, or when averaging over several samples, the obtained value of Δ can deviate from zero, characterizing the error of the experimental method for calculating the DOPT. In particular, due to the relatively large number of cells considered in Cucumis epithelium [START_REF] Lewis | The effect of cell division on the shape and size of hexagonal cells[END_REF], the error (1) for this case is ~0.004, which is almost 10 times lower than for the data in ref. [START_REF] Gibson | The emergence of geometric order in proliferating metazoan epithelia[END_REF]. Nevertheless, in the geometrically correct model of cell division [START_REF] Gibson | The emergence of geometric order in proliferating metazoan epithelia[END_REF], with an increase in the number of successive divisions, Δ tends to 0, reaching Δ≈0.001 at 10 th division. Since both the left and right sides of the distribution contribute to Δ, it is reasonable to estimate the maximum error in determining the probabilities 𝑃 " in the DOPT as Δ/2. Also, the condition Δ=0 severely restricts the possible shapes of the DOPT. In real epithelia, probabilities Pi, where i>7 or i<5, are small. The critical probability is 𝑃 % , and other probabilities follow it conserving the condition Δ=0. For example, if in a hypothetical planar epithelium consisting only of 5-,6-and 7-valent cells, the percentage of hexagonal cells is P6 then concentrations of 5-and 7-valent cells should be equal to (1-P6)/2.

HeLa cells are human malignant epithelial cells derived from an epidermoid carcinoma of the cervix. The growth of confluent HeLa cell monolayers and synchronization procedure are described in Methods, see also electronic supplementary material. In synchronized monolayers, most of the cells belong to the second generation with similar time elapsed after the division process. The characterization of the HeLa and HeLa synchronized epithelia (see figure 1) was carried out using 9 and 7 assembled images obtained by the juxtaposition of contiguous microscope fields. The studied epithelial areas contained from 404 to 933 cells. The first line of figure 1 shows typical non-synchronized and synchronized HeLa epithelial cells (respectively samples HeLa9 and HeLasyn5 in Table I). Due to visualization specificities (see Methods), the cell nuclei in micrographs are clearly visible, while the cell boundaries are indistinguishable in most cases. Therefore, in order to determine the number of nearest neighbors and obtain additional structural data, we used the Voronoi tessellation (see Methods) with the nodes located at the centers of the nuclei as shown in the second line of figure 1. That is, strictly speaking, we analyzed not the valency of cells, but the valency of their nuclei. The areas of the epithelial cells were also calculated as the areas of the cells of the Voronoi tessellations. I. In all cases, the averaged topological error Δ is less than 0.01.

To study correlation between DOPTs and the spread of cell sizes we analyzed histograms of cell area distribution for the considered samples (see figure 1(e,f)). All the histograms (including the samples not shown in figure 1) are wide and similar to the Gauss type, which can be explained by the spread of the cell sizes before the mitosis, and the possibility of cell division into substantially unequal parts. Asymmetry of distributions, in our opinion, is due to the fact that in the considered cellular structures, the minimum cell area is limited, not by zero, but by a specific positive value. We have also noticed that the maximum spread of areas strongly fluctuates and can noticeably differ in structures with similar morphology. Therefore, it is reasonable to characterize a cellular structure with the average area 𝑆 &' of its cells and the average spread ΔS, which we define as the difference between the smallest and the largest cell areas among the half of the cells with the areas closest to the 𝑆 &' value. We then introduce a dimensionless effective spread Δ𝑆/𝑆 &' .

Structural data on the investigated epithelia are presented in Table I. The data in the last two lines are averaged for all the, respectively, synchronized (<Helasyn>) and non-synchronized (<HeLa>) HeLa samples studied. Note, that in 2-4 th columns, the values are weighted arithmetic means. Namely, the Sav, Δ𝑆/𝑆 &' and P6 values are weighted by the numbers of Nvor in each line.

The average cell area in the HeLa non-synchronized epithelium is larger than in the synchronized one (see Table I), since most of the cells in the latter belong to the second generation with similar time elapsed after the division process. In all the examined specimens of both types, the Δ𝑆/𝑆 &' value was substantial. The spread in the 𝑆 &' value is apparently due to the growth of samples upon coverslips with an uneven surface. Therefore, the cellular monolayer undergoes a strain, which we consider to be homogeneous in the image size scale. The spread in the average cell areas between the samples of the same type can be also related to the following experimental feature. To prevent the formation of multilayers, the growth is stopped just before the total confluence of the cells. As a result, small and relatively sparse empty areas appear (see Methods). Note also that the averaged value < Δ𝑆/𝑆 &' > is 7% larger in the HeLa non-synchronized data set. However, the averaged DOPTs remain very close, with the differences in probabilities Δ𝑃 " ≲ 0.01 (see figure 1g).

Table I. Characterization of the studied samples. The columns contain: sample name, average area of Voronoi cells Sav, dimensionless effective spread of areas Δ𝑆/𝑆 "# , probability P6, number of identified Voronoi cells Nvor, total number of cells Ntot. The last two lines in bold correspond to the averaged structural data (including standard deviations) for HeLa and HeLa synchronized samples.

Code

Sav (μm As we have already mentioned, synchronized cell monolayers are beyond the scope of applicability of the Markov-type theory [START_REF] Gibson | The emergence of geometric order in proliferating metazoan epithelia[END_REF]. Its inapplicability for the non-synchronized case is evidenced by our finding that DOPTs in both types of monolayers are very similar. Differences in structural parameters of different samples, substantial spread of the cell sizes, and the specific asymmetry of the cell area distributions lead to the hypothesis that in the hyperproliferative epithelia the cell packings are close to random but, nevertheless, satisfy the geometric constraint associated with the existence of minimal cell size. Below we develop the theory of random polygonal packings and then test our hypothesis on epithelial cancer cell monolayers.

To construct the random polygonal packings, we start from random distributions of points with minimal allowed distance between them, dmin. We use periodic boundary conditions (which ensure Δ=0) and follow the random sequential adsorption algorithm [START_REF] Feder | Random sequential adsorption[END_REF]. Points are randomly and sequentially inserted into the fundamental region with the area Ar. If the distance between the point that is currently being placed and any of the points placed earlier is less than the distance dmin, then this point is deleted, and the random insertion is repeated. When the desired N number is achieved, the insertion is stopped. The subsequent Voronoi tessellation yields polygons with the average area 𝑆 &' = 𝐴 ( /𝑁.

Obviously, there exists an averaged maximum limit on the number of points Nmax that can be randomly placed in the region with area Ar for the given minimal distance dmin. Thanks to the discovery of new types of disordered structures and the continuous increase of available computing power, this limit, also known as the jamming limit [START_REF] Feder | Random sequential adsorption[END_REF], has been extensively studied in the literature over the past 30 years for various shapes and dimensions of adsorbed particles [START_REF] Feder | Random sequential adsorption[END_REF][START_REF] Penrose | Random Parking, Sequential Adsorption, and the Jamming Limit[END_REF][START_REF] Evans | Random and cooperative sequential adsorption[END_REF][START_REF] Talbot | From car parking to protein adsorption: an overview of sequential adsorption processes[END_REF]. In the case of the random sequential adsorption of equivalent disks, the average ratio of the area occupied by them to the total surface area tends to 𝐿 ≈ 0.547 [START_REF] Feder | Random sequential adsorption[END_REF]. Naturally, in our model the average ratio 𝜋𝑁 )&* 𝑑 )"+ , /𝐴 ( tends to the same value.

If a packing is not random, its surface coverage can differ significantly. For example, for the densest hexagonal packing of equivalent disks the surface coverage equals -, √ / ≈0.9064. Let us consider the densest hexagonal packing with N disks per area Ar. In this packing dmin=dhex, where dhex= ? 01 ! 2 √ / is the distance between the disk centers. Then, in our packing algorithm the jamming limit corresponds to dmin= ? Note that to model the hyperproliferative epithelia we use the unsaturated structures, where the ratio 𝜂 =dmin/dhex is smaller than the above critical value 0.776. Figure 2a-f shows examples and area distribution histograms of random polygonal packings with different ratios 𝜂 and different degrees of hexagonality. The plots in figure 2g are calculated up to η=0.75, since at larger η the calculation time increases sharply. Note also, that packings with the same η can have slightly different morphology, so the averaging is needed. In particular, when N~5000, in the region where η>0.4, the ratio Δ𝑆/𝑆 &' and probabilities Pi are reproduced with standard deviations smaller than 0.01 (see figure 2g).

The proposed random packing model is in perfect agreement with the averaged structural data for the non-synchronized and synchronized HeLa epithelial cell monolayers (see the last two lines of Table I (<Hela> and <Helasyn>) and the histograms in figure 1g). Indeed, as presented in figure 2g, the change in Δ𝑆/𝑆 &' from 0.4 to 0.43 corresponds to the variation of η from 0.5 to 0.47, respectively. In this region, the slope of the plot P6(η) is small, and the values of η correspond to very close DOPTs with 𝑃 % ≈ 0.35 -0.37. The deviations of other averaged probabilities (see figure 1g) from their theoretical values also do not exceed 0.02, which is within the spread of model calculations at N~5000. The random packing model can also explain the scatter of Δ𝑆/𝑆 &' and P6 values between different samples presented in Table I. Note that these monolayers contain, on average, slightly less than 500 cells, and such a small value of N increases the morphological inhomogeneity of the generated packings. To justify this, one can perform a series of computations with the appropriate inputs: the number of calculations should not be less than the total number of monolayers, and N is equal to Nvor in the sample. As a result, it is highly probable that a packing that deviates from the mean by about (or even more) than the considered sample will be generated.

As we already mentioned, small and relatively sparse empty domains are present in nonsynchronized HeLa monolayers. We decided to evaluate the impact of these domains on the obtained results. For this purpose, on the assembled images, we selected 40 smaller regions without empty areas. These confluent rectangular regions contained from 36 to 76 cells. When averaging over these regions, the values of 𝑆 &' were found separately, and cells for which it was impossible to determine the number of nearest neighbors were not considered. This treatment resulted in Δ𝑆/𝑆 &' ≈0.394 and 𝑃 % ≈ 0.38. Comparing these values with those for averaged <HeLa> and using the graphs shown in figure 2g, we observe that the decrease in Δ𝑆/𝑆 &' corresponds to the small increase in 𝑃 % , and the order is still properly described by the developed theory of random polygonal packings. So, we conclude that both types of considered hyperproliferative epithelia represent the random packings, while the main structural difference between the non-synchronized and synchronized HeLa monolayers consists in different values of 𝑆 &' .

Discussion

In the literature, one can find many biophysical approaches describing the correlation between epithelial properties and averaged cell area and perimeter (see, for example, [START_REF] Farhadifar | The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing[END_REF][START_REF] Puliafito | Collective and single cell behavior in epithelial contact inhibition[END_REF][START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF][START_REF] Bi | A density-independent rigidity transition in biological tissues[END_REF][START_REF] Giavazzi | Flocking transitions in confluent tissues[END_REF][START_REF] Li | Mechanisms of regulating cell topology in proliferating epithelia: Impact of division plane, mechanical forces, and cell memory[END_REF][START_REF] Classen | Hexagonal packing of Drosophila wing epithelial cells by the planar cell polarity pathway[END_REF]). The spread in these geometrical parameters is often considered as a 'biological noise' and its importance for the epithelium is discussed in only handful of papers [START_REF] Roshal | Crystal-like order and defects in metazoan epithelia with spherical geometry[END_REF][START_REF] Classen | Hexagonal packing of Drosophila wing epithelial cells by the planar cell polarity pathway[END_REF][START_REF] Atia | Geometric constraints during epithelial jamming[END_REF][START_REF] Ballard | The energy landscapes of bidisperse particle assemblies on a sphere[END_REF]. In particular, the recent paper [START_REF] Atia | Geometric constraints during epithelial jamming[END_REF] demonstrates that the relation between average aspect ratio of the cells comprising epithelium and its standard deviation governs processes as diverse as maturation of the pseudostratified bronchial epithelial layer cultured from non-asthmatic or asthmatic donors, and formation of the ventral furrow in the Drosophila embryo. Interestingly, in the paper [START_REF] Classen | Hexagonal packing of Drosophila wing epithelial cells by the planar cell polarity pathway[END_REF] a possibility of a similar connection between cell area variability and DOPT was only briefly mentioned. In our paper, as far as we know, we have for the first time proposed a model that establishes the relation between the cell size spread Δ𝑆/𝑆 &' and DOPT in the cancer epithelial monolayers. Examining non-synchronized and synchronized HeLa cervix cells, we have shown that this spread is the main parameter controlling the DOPT in both types of the monolayers that represent a real-life example of random packings.

Below we discuss how a healthy proliferative epithelium distinguishes from random cancer monolayers and point out a physical mechanism underlying this difference. As one can see from figure 2g, at η~0.62-0.7, the model of random polygon packing reproduces perfectly the DOPTs [START_REF] Lewis | The effect of cell division on the shape and size of hexagonal cells[END_REF][START_REF] Korn | The geometry of plant epidermal cells[END_REF][START_REF] Kwiatkowska | Growth and morphogenesis at the vegetative shoot apex of Anagallis arvensis[END_REF][START_REF] Gibson | The emergence of geometric order in proliferating metazoan epithelia[END_REF] typical of many plant and animal proliferative epithelia: P4≈ 0.02 -0.03, P5≈ 0.25 -0.27, P6≈0.41-0.47, P7≈ 0.22, P8≈ 0.03 -0.05, and P9≈ 0.001 -0.006. A more detailed analysis, however, shows that the proposed approach leads to a value of Δ𝑆/𝑆 &' that is slightly smaller than the one observed experimentally. In particular, after analyzing the images of the Cucumis proliferative epithelium [START_REF] Lewis | The effect of cell division on the shape and size of hexagonal cells[END_REF], we have estimated the value of Δ𝑆/𝑆 &' as 0.29, which in the random polygon packing model corresponds to η~0.63 and P6≈0.43 instead of the observed value P6≈0.47 [START_REF] Lewis | The effect of cell division on the shape and size of hexagonal cells[END_REF].

Let us consider the situation in more detail using our experimental data on healthy proliferative Human Cervical Epithelial Cells (HCerEpiC). Figure 3a-b demonstrates a typical sample of HCerEpiC confluent monolayer with the corresponding Delaunay triangulation and Voronoi tessellation. Forty analogous experimental images, previously obtained in our laboratory, capture relatively small separate fragments showing simultaneously from 30 to 54 cells.

Note that borders of the HCerEpiC cells are more clearly visible than those of HeLa cells, and thus can be determined directly, without using Voronoi tessellation. First, we analyzed the images following the same procedure we used for HeLa cells. In this way we have obtained the following structural data: 𝑆 &' ≈3.0x10 3 μm 2 , Δ𝑆/𝑆 &' ≈ 0.43, P4≈0.04, P5≈0.30, P6≈0.41, P7≈0.20, P8≈0.04, and P9≈0.01. It is interesting to note that these Pi values differ from those obtained for the Xenopus frog [START_REF] Gibson | The emergence of geometric order in proliferating metazoan epithelia[END_REF] by less than 2% and the HCerEpiC epithelium, as expected, has a structure typical of other normal proliferative epithelia [START_REF] Lewis | The effect of cell division on the shape and size of hexagonal cells[END_REF][START_REF] Korn | The geometry of plant epidermal cells[END_REF][START_REF] Kwiatkowska | Growth and morphogenesis at the vegetative shoot apex of Anagallis arvensis[END_REF][START_REF] Gibson | The emergence of geometric order in proliferating metazoan epithelia[END_REF]. We then analyzed the same forty images directly, without using Voronoi tessellation. Although cells with common borders are not necessary linked in the Delaunay triangulation (discrepancy ≈15% of cases), the resulting difference in the calculated average structural parameters is several times smaller indicating compatibility of these two methods. Therefore, below we refer only to the structural data obtained using Voronoi tessellation.

We have also considered the fact that the 𝑆 &' value can vary from one HCerEpiC monolayer fragment to another, because of the influence of the substrate it was grown on. With the averaging method taking this difference into account, the ratio Δ𝑆/𝑆 &' decreases from 0.43 down to 0.36, corresponding in the random packing model to P6≈0.37 (instead of 0.37). Nevertheless, this slightly increased probability is still smaller than the observed value 𝑃 % ≈0.41. Consequently, based on the Cucumis and HCerEpiC examples, we can conclude that in normal proliferative epithelia the 𝑃 % value is greater than the one predicted by the random packing model for the same Δ𝑆/𝑆 &' ratio. Note that the mitosis rate in HeLa cells is ~5.5 times higher than in HCerEpiC cells, while the rate of apoptosis is the same in both cases (see electronic supplementary material, figure S1). Due to this difference, for the same number of seeded cells, the considered monolayers are at confluence after 2 days for HeLa cells compared to 4 days for HCerEpiC ones. Thus, we can assume that the 𝑃 % increase in healthy epithelium is associated with the lower rate of cell division, although, intercellular interactions are also very important for this phenomenon.

To justify the latter statement, we recall that ordering of equivalent particles, retained on a planar surface and interacting with each other via very different pair potentials, readily leads to the formation of a simple hexagonal order [START_REF] Roshal | Extended topological defects as sources and outlets of dislocations in spherical hexagonal crystals[END_REF]. In our opinion, the mechanism of the 𝑃 % increase in healthy proliferative epithelia is similar and caused by minimization of the elastic energy associated with the mechanical interaction between cells. Namely, an internal local stress caused by the cell division and nonuniform growth can relax through cell motility, which increases, in average, the hexagonal coordination.

Below, we consider a simple way to demonstrate that mechanical interactions between cells can order the epithelial structure by increasing the number of cells with six neighbors. Following the works [START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF][START_REF] Bi | A density-independent rigidity transition in biological tissues[END_REF], the elastic deformation energy E of a monolayer containing N cells can be written as:

𝐸 = G 𝛽 "( 𝐴 " -𝐴 " 6 ) , + 𝜁 "( 𝑃 " -𝑃 " 6 ) , , 2 "78 (2)
where the subscript i labels each cell; Ai and Pi are the cell area and perimeter, respectively. The origin of the elastic moduli 𝛽 " is associated with a combination of the cell volume incompressibility and resistance to height differences between nearest cells [START_REF] Bi | A density-independent rigidity transition in biological tissues[END_REF]. The second term including the cell perimeters 𝑃 " results from active contractility of the actomyosin subcellular cortex (quadratic in 𝑃 " ) and effective cell membrane tension due to cell-cell adhesion and cortical tension, which are linear in perimeter. Usually, different cells of the same type are assumed to have the same elastic moduli 𝛽 " and 𝜁 " and thus index 'i' is omitted. The same goes for the parameters 𝑃 " 6 and 𝐴 " 6 , which are usually substituted with 𝑃 599 and 𝑆 &' , correspondingly. The resulting energy reads:

𝐸 = G 𝛽 ( 𝐴 " -𝑆 &') , + 𝜁 J 𝑃 " -𝑃 599 K , . 2 "78 (3)
In [START_REF] Bi | A density-independent rigidity transition in biological tissues[END_REF], using the Surface Evolver software [START_REF] Brakke | The Surface Evolver[END_REF], the energy (3) was minimized with respect to the position and shape of the intercellular boundaries, while in [START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF] an overdamped dynamic model was built. The latter model takes into account the effective temperature of the monolayer, which corresponds to the value of additional (not related to intercellular interaction (3)) random forces simulating active cell migration. In the model [START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF], the shape of the cells coincides with the shape of Voronoi cells therefore the energy (3) depends only on positions of the Voronoi cells centers.

The results of the works [START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF][START_REF] Bi | A density-independent rigidity transition in biological tissues[END_REF] are similar, and the properties of both models depend on so-called target shape index 𝑞 = : "##

; < $% [START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF][START_REF] Bi | A density-independent rigidity transition in biological tissues[END_REF]. To clarify the meaning of this quantity, note that for the regular honeycomb packing the quantity 𝑞 = % ; 8.> √ / ≈ 3.72, while q = 4 and 𝑞 ≈ 4.56 minimize the energy (3) for regular square and triangular lattices, respectively. As shown in the works [START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF][START_REF] Bi | A density-independent rigidity transition in biological tissues[END_REF] at q>q0, where q0=3.81-3.813, the cells acquire the possibility of relative movement without overcoming the potential barrier, and the monolayer demonstrates the fluid-like behavior [START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF][START_REF] Bi | A density-independent rigidity transition in biological tissues[END_REF]. At the same time, as can be seen from the model monolayers presented in these works, at q>q0 the number of 6-valent cells decreases, and many cells acquire an elongated shape.

Since our analysis of experimental data is also based on the Voronoi tiling, let us discuss the results [START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF] in more detail. Image analysis (see figure 2b in [START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF]) reveals that Δ𝑆/𝑆 &' are 0.1 and 0.12, respectively, for the solid and liquid states of the model monolayers shown. However, even in Cucumis epithelium, the most ordered proliferative monolayer considered in our work, Δ𝑆/𝑆 &' ≈0.29, which is almost 3 times more than in the model structures [START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF]. Also note that the probability of 𝑃 % in the solid phase [START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF] (its comparing with normal epithelium is reasonable) is 0.595, whereas for ordinary proliferative epithelia, the value of 𝑃 % lies in the range of 0.41-0.47 [START_REF] Lewis | The effect of cell division on the shape and size of hexagonal cells[END_REF][START_REF] Korn | The geometry of plant epidermal cells[END_REF][START_REF] Kwiatkowska | Growth and morphogenesis at the vegetative shoot apex of Anagallis arvensis[END_REF][START_REF] Gibson | The emergence of geometric order in proliferating metazoan epithelia[END_REF].

To assess the influence of cell-to-cell interactions on the DOPT in the model [START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF], we performed a series of numerical experiments, where the energy (3) was minimized using the ordinary coordinate descent method. As the initial positions of N=400 cell centers, we used random coordinates obtained within the framework of the above algorithm for constructing random polygonal packings. The analysis showed that mechanical intercellular interactions can strongly order the system. In the region q<q0, depending on the initial positions of the particles and the ratio of the coefficients 0.01< 𝛽𝑆 &' /𝜁<100, minimization of (3) led to the probabilities 𝑃 % = 0.6 -0.8 (see figure 3e). The scatter of the obtained 𝑃 % values is associated not only with the change in the model parameters, but also with the existence of numerous local energy minima. Note that the anomalous epithelia with 𝑃 % ≈0.6-0.9 are also observed in nature. For example, after the mitosis in the wings of a Drosophila stops (and before the start of hair growth), the equilibrium changes and a sharp increase in the proportion of 6-valent cells occurs [START_REF] Classen | Hexagonal packing of Drosophila wing epithelial cells by the planar cell polarity pathway[END_REF]. As our preliminary estimates show, the use of Lennard-Jones potential [START_REF] Martinand-Mari | Topological control of life and death in non-proliferative epithelia[END_REF] to describe the intercellular interaction and assumption of the cell size dispersion [START_REF] Roshal | Crystal-like order and defects in metazoan epithelia with spherical geometry[END_REF] allow to obtain epithelium-like structures with 𝑃 % up to the level of 0.9.

Rationalization of normal proliferative epithelia within the framework [START_REF] Bi | Motility-driven glass and jamming transitions in biological tissues[END_REF] is problematic not only due to the significant effective temperature and monolayer disordering required to reduce the probability 𝑃 % from 0.6-0.8 to 0.4-0.47. The assumption of the same equilibrium area 𝑆 &' for all cells (made in equation ( 3)) limits the ability of the model to generate structures with realistic values of Δ𝑆/𝑆 &' . In this context, transition from energy (2) to (3) by replacing 𝐴 " 6 , 𝑃 " 6 with corresponding averages seems to be an oversimplification; thus, we use the following energy instead:

𝐸 = G 𝛽 ( 𝐴 " -𝐴 " 6 ) , + 𝜁 U 𝑃 " -𝑞 ? 𝐴 " 6 V , , (4) 2 "78
where the equilibrium areas 𝐴 " 6 are defined by our random polygonal packing model. Methods to analyze an energy landscape are known [START_REF] Ballard | The energy landscapes of bidisperse particle assemblies on a sphere[END_REF][START_REF] Wales | Energy landscapes: Applications to clusters, biomolecules and glasses[END_REF], nevertheless, for energy (4) it is a separate complex problem, since the energy (4) has a huge number of local minima with different DOPTs. We have only verified that minimization of this energy approximately preserves the initial value of Δ𝑆/𝑆 &' , and in the region 𝑞 < 𝑞 6 , leads to an increase in 𝑃 % that is smaller than the one occurring during the minimization of energy [START_REF] Azzag | The non-proliferative nature of ascidian folliculogenesis as a model of highly ordered cellular topology distinct from proliferative epithelia[END_REF].

Using the energy (4), one can more realistically simulate normal proliferative epithelia, in particular HCerEpiC epithelium. T he DOPT in HCerEpiC epithelium corresponds to the random polygonal packing with 𝜂 ≈0.62, while the observed Δ𝑆/𝑆 &' corresponds to 𝜂 ≈0.53. To correct the DOPT, we generated a random polygonal structure with 𝜂 ≈0.53 (see figure 3d), that yields a distribution of cell areas very close to that of the HCerEpiC epithelium (see figure 3c). At q=3.81, we minimized energy (4) at different ratios 𝛽𝑆 &' /𝜁 in the region 0.01 -100. This gave rise to structures with 𝑃 % = 0.38 -0.43. One of these structures, closest to the HCerEpiC epithelium, is shown in figure 3f.

Finally, we note that in epithelia there are multiple processes brining system out of equilibrium such as mitosis, apoptosis, varying cell growth rate. Their influence on DOPT cannot be taken in account by the simple minimization of energy [START_REF] Farhadifar | The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing[END_REF]. The construction of such a more complex dynamic model is clearly beyond the scope of this work and is a task for the next studies.

Conclusion

In this paper, we have demonstrated the significant difference in topology of cancerous and normal epithelial monolayers. In epithelial cancer monolayers, the cell arrangement is close to random but, nevertheless, satisfies the geometric constraint associated with the existence of minimal cell size. This type of the cell order can be described with a good accuracy by a single dimensionless control parameter, namely the normalized halfwidth Δ𝑆/𝑆 &' of the cell area distribution. Cancer cells divide so quickly that the relaxation processes associated with the minimization of free energy are unable to effectively rearrange cells and order the monolayer by increasing the number of 6valent cells. Consequently, the growing disorder predominates over the intercellular interactions, and this yields structures very similar to random polygonal packings. In contrast, normal proliferative and non-proliferative epithelia are morphologically different from them. Overall, our study of epithelial monolayers provides novel mathematical tools for a topological analysis of epithelial morphogenesis in a developmental and/or pathological context. In particular, elaborated tools could be used to analyze the change in cell topology during the 2D/3D transition observed throughout epithelial tumorigenesis. In the near future, we hope to apply our approach to studing the zebrafish embryonic epidermis development before and after tumor cell injection in vivo.

Methods

Cell line growth and synchronization procedure conditions.

Human cervical cancer cell line HeLa was purchased from the American Type Culture Collection (ATCC; Manassas, VA, USA) and maintained in DMEM, high glucose (Dulbecco's Modified Eagle Medium) containing 5% heat-inactivated foetal bovine serum and supplemented with GlutaMAX TM (Gibco Life Technologies), penicillin (100 units/mL), and streptomycin (100 µg/mL). Normal primary cervical epithelial cells (HCerEpiC) isolated from human uterus were purchased from ScienCell Research Laboratories (Clinisciences S.A.S., Nanterre, France). HCerEpic cells were grown in Cervical Epithelial Cell Medium according to the manufacturer's instructions. Cells were incubated in a humidified incubator at 37°C in 5% CO2.

For confocal microscopy analysis of non-synchronized cells, HeLa or HCerEpiC cells were seeded on glass coverslips (12mm diameter round) coated with 10µg/ml of poly-L-Lysine (P4707, Sigma) at 7•10 4 cells/coverslip in a 24-well culture plate. Confluency of the monolayer was achieved 48 h and 4 days later for HeLa and HCerEpiC cells, respectively. HeLa cells were synchronized in G0/G1 phase of the cell cycle using the double thymidine block procedure. The day before the first thymidine block, cells were seeded on poly-L-Lysine treated glass coverslip at a density of 7.5x10 4 cells/ coverslip. The next day, 2.5 mM thymidine was added for 16 h (first block). Then cells were washed with phosphate-buffered saline (PBS) and incubated during 8 h without thymidine. Lastly, the second thymidine block was applied for 16h. At the end of the second thymidine block, cells were washed and incubated for 2 additional hours, after which cells were treated with antibodies and analysed by confocal microscopy. Following this procedure, cells are fully confluent and cell synchronisation in G0/G1 is about 92%. This was confirmed by cytofluorometry analysis (FACS) (see electronic supplementary material, figure S2).

Immunocytochemistry and fluorescence microscopy.

Cells cultured on glass coverslips were fixed with 4% paraformaldehyde in PBS for 20 min, and washed in Tris-buffered saline (25mM Tris pH7.4, 150mM NaCl) (TS) for 10 min. After permeabilization with 0.2% Triton X-100 in TS for 4 min, non-specific binding was blocked with 0.2% gelatin from cold water fish skin (#G7765 Sigma-Aldrich Chimie, Lyon) in TS for 30min. Cells were incubated with primary antibodies in blocking buffer for 1h and then washed 3 times with 0.008% TritonX-100 in TS for 10 minutes. Rabbit anti-ezrin antibody [START_REF] Andréoli | Ezrin has properties to selfassociate at the plasma membrane[END_REF] was used to visualize cell body and membrane. Cells were incubated for 30 minutes with Alexa-Fluor 488labelled secondary antibodies (P36934-Molecular Probes, InVitrogen) in blocking buffer. After rinsing in washing buffer, cell nuclei were stained with 1 µg/ml Hoechst 33342 (62249-Thermo Scientific Pierce) in TS for 5 minutes. Finally, coverslips were mounted with Prolong TM Gold Antifade (P36934-Molecular Probes, InVitrogen) and examined under a Leica TCS SPE confocal microscope equipped with a 25X/0.75 PL FLUOTAR oil objective (HCerEpiC) and a 40X/1.15 ACS APO oil objective (HeLa) (figure 3 and electronic supplementary material, figure S1). In the case of HeLa cells (figure 1a) and synchronized HeLa cells (figure 1b), the analysis was performed with a Zeiss LSM880 FastAiryScan confocal microscope equipped with a 40X/1.4 Oil Planapochromat DIC objective. Between 7 to 9 'z-stacks' (0.457µm thickness each) were acquired per field and a 2D image was generated by applying Maximum Intensity Projection processing. For each coverslip, the acquisition pattern was 6 neighbouring images per row for a total of 2 or 3 rows. The resulting images (12 or 18 images) were adjusted for brightness, contrast and color balance by using ImageJ and assembled side by side in PowerPoint to reconstruct a cell monolayer consisting of N>500 cells.

Image analysis.

After determining the geometric centers of the cell nuclei, triangulation was performed by the Delaunay method [START_REF] Delaunay | Sur la sphère vide: A la mémoire de Georges Voronoi (French)[END_REF]. Next, Voronoi tiling was constructed, and the areas of the epithelial cells were calculated as the areas of Voronoi cells. Obviously, for a correct statistical analysis, it is necessary to discard the cells (located too close to the image border), the number of neighbors for which cannot be determined. Note that even if it is possible to construct a closed Voronoi cell, then it is also necessary to check whether the cell polygon boundary can be changed by additional hypothetical nuclei lying directly outside the image border. Therefore, the center of a reliably constructed Voronoi cell should be located at least twice as far from the image border as any of the vertices of this cell. However, this method leads to the appearance of an excessive total positive topological charge, which is localized at the image border. On one hand, 4-and 5-valent cells have a smaller area [START_REF] Lewis | The effect of cell division on the shape and size of hexagonal cells[END_REF], while on the other hand, the smaller the cell located near the image border, the more chances its nucleus has to satisfy the selection criterion formulated above. This fact, when processing images with a small number of nuclei (about 40), leads on average to the formation of a 5% preponderance of the total positive topological charge (which is carried by 4-and 5-valent cells) over the total negative topological charge. This, in turn, leads to errors when constructing DOPTs diagrams and determining the value of 𝑃 % and Δ𝑆/𝑆 &' . To avoid preferential selection of small cells, we used additional cutting of the image borders. In the statistical analysis, we took into account only cells whose nuclei centers fall within the rectangle, which has maximum possible size and does not contain any nucleus with an uncertain coordination. Thus, it is possible to significantly reduce the total topological charge of the images and, accordingly, the error in the values of Δ𝑆/𝑆 &' .

Figure 1 .

 1 Figure 1. Structural characterization of HeLa cell monolayers. (a) Non-synchronized and (b) synchronized cellular structures. White scale bars are 100 μm. The triangulation with nodes in the centers of cell nuclei is imposed on the monolayers. (c,d) Voronoi tessellations for the monolayers (a) and (b), respectively. Histograms (e,f) show a probability P to observe the certain ranges of cell areas S in the samples (a,b). Polygon types and corresponding contributions to the histograms are color-coded as shown in the legend. (g) DOPTs for the samples (a,b) and averaged data. Dark blue, blue, red, and pink colors correspond to HeLa9, <HeLa>, HeLasyn5, and <HeLasyn> lines from TableI. In all cases, the averaged topological error Δ is less than 0.01.

𝑑 45 *

 45 ≈0.776𝑑 45* .

Figure 2 .

 2 Figure 2. Examples of random packings and their geometric characteristics. (a-d) Packings obtained at η values equal to 0.1, 0.4, 0.6, and 0.758, respectively. It can be observed how the fraction of hexagonal cells grows with increasing η (these cells are shown in yellow; the coloring is the same as in figures 1c-f). (e-f) Area distribution histograms obtained for η=0.1 and 0.758. For both cases, the value of 𝑆 "# is renormalized to 1. (g) Dependences Δ𝑆/𝑆 "# and Pi on η. Plots of Δ𝑆/𝑆 "# , P6, P5, P7, P4, P8, P9, calculated with the step Δη=0.05, are colored with black, yellow, orange, green, red, light blue and dark blue colors, respectively. Probabilities P3, P10 and P11 are too small to be shown in the chosen scale. Each center of vertical bars represents the averaging of 10 calculations at N=5000. The bar sizes denote the standard deviations obtained for the calculations.

Figure 3 .

 3 Figure 3. Structural characterization of HCerEpiC epithelium. (a) Typical sample of the epithelium with the superimposed triangulation; white scale bar is 50 μm. Vertices of the Delaunay triangulation coincide with the cell nuclei. The borders between the cells are shown in green. (b) Dual Voronoi tessellation for the same sample. A cell was taken into account in the statistical analysis if its shape could be determined unambiguously and the cell vertices were not too close to the border of the image (see Methods). Such cells are colored. (c) Averaged histogram of cell areas in forty samples. (d) Voronoi tiling for a random polygonal structure with the same Δ𝑆/𝑆 "# value and very similar to (c) histogram of the cell area distribution. The polygon coloring is the same as in figures 1c-f (e,f). More ordered structures obtained from (d) by minimizing the model energies of elastic intercellular interaction (see main text). (e) Minimization of energy (3) at 𝛽𝑆 &' /𝜁=1 and q=3.81 results in a highly ordered structure with 𝑃 & ≈ 0.7. (f) Minimization of energy (4) at 𝛽𝑆 &' /𝜁=10 and q=3.81 yields the packing structurally similar to HCerEpiC epithelium: Δ𝑆/𝑆 "# ≈ 0.45 and 𝑃 & ≈ 0.41.

. Model of random polygon packing.

			𝚫𝑺/𝑺 𝒂𝒗	P6	Nvor	Ntot
	HeLa1	839.0	0.444	0.318	402	578
	HeLa2	1024.3	0.388	0.373	263	404
	HeLa3bis	951.2	0.431	0.349	373	561
	HeLa4	942.8	0.501	0.336	318	532
	HeLa5	831.0	0.462	0.348	399	623
	HeLa6bis	1057.9	0.388	0.343	539	742
	HeLa8	987.9	0.411	0.371	582	780
	HeLa9	813.1	0.441	0.352	691	915
	HeLa10	923.9	0.444	0.343	577	828
	HeLasyn1	499.8	0.452	0.379	688	933
	HeLasyn3	628.9	0.418	0.337	591	826
	HeLasyn4	853.9	0.391	0.325	409	593
	HeLasyn5	696.5	0.386	0.355	512	714
	HeLasyn6	970.3	0.367	0.371	353	525
	HeLasyn8	815.9	0.378	0.316	396	624
	HeLasyn9	693.9	0.391	0.372	521	739
	<HeLa>	924.9±100.0	0.433±100.0	0.348±100.0	460.4	662.6
	<Helasyn>	705.6±100.0	0.403±100.0	0.353±100.0	495.7	707.7
	2.2

2 )

Commenté [РСБ1]: Поднять масштаб (a-d) процентов на 30-40, в 3-ю строчку в два столбика добавить исправленную по цветам легенду в 2 столбца. Убрать не наблюдаемые типы из легенды

Commenté [РСБ2]: Посчитай и впиши в таблицу.

Acknowledgments. M. M., V.M., and S.B. acknowledge the imaging facility MRI, member of the national infrastructure France-BioImaging infrastructure supported by the French National Research Agency (ANR-10-INBS-04, «Investments for the future»). S.R. and D.R. acknowledge the financial support from the Russian Foundation for Basic Research (Grant № 19-32-90134). I.G. acknowledge Russian Science Foundation for the financial support of the research (Grant № 20-72-00164). We thank A. Parmeggiani for useful discussions.

Author contributions. To explain the topological universality of normal epithelia, S.B. proposed to study the cancer cell monolayers. S.R. developed the model of random polygonal packings and controlled the theoretical part of the project. M. M. performed the biological experiments and the confocal microscopy acquisition as well as assembly of images. K.F. realized the initial analysis of images. D.R. did most of the calculations and analysis. I.G. reproduced the SPV computer model and found the relation between our approach and RSA algorithm. S.R. wrote with the essential help of S.B. I.G. and M.M. the main part of the text. V.M. proofread and corrected the article extensively. All authors participated in discussion and writing of the article.

Competing interests. The authors declare no competing interests.


Hussein El Khansa
email: hussein.el-khansa@ird.fr

Carmen Gervet
email: carmen.gervet@umontpellier.fr

Audrey Brouillet
email: audrey.brouillet@ird.fr

Hussein Al Khansa

PROMINENT DISCORD DISCOVERY WITH MATRIX PROFILE : APPLICATION TO CLIMATE DATA INSIGHTS

Keywords: Prominent discord discovery, Large time series, Matrix profile, Climate data

The definition and extraction of actionable anomalous discords, i.e. pattern outliers, is a challenging problem in data analysis. It raises the crucial issue of identifying criteria that would render a discord more insightful than another one. In this paper, we propose an approach to address this by introducing the concept of prominent discord. The core idea behind this new concept is to identify dependencies among discords of varying lengths. How can we identify a discord that would be prominent? We propose an ordering relation, that ranks discords and we seek a set of prominent discords with respect to this ordering. Our contributions are 1) a formal definition, ordering relation and methods to derive prominent discords based on Matrix Profile techniques, and 2) their evaluation over large contextual climate data, covering 110 years of monthly data. The approach is generic and its pertinence shown over historical climate data.

INTRODUCTION

The analysis of climate data towards the extraction of global climate trends using ensemble mean approaches is receiving a wide interest. To this date, the search for pattern anomaly or outlier patterns has received less attention in climate data studies. Such data is contextual due to the geo-localization and timestamps of the data series relative for instance to soil humidity, temperature or rainfall. It comes from historical sources or complex simulation impact models (physical processes of atmosphere and ocean) such as Earth System Models ( [START_REF] Randall | Climate models and their evaluation, Climate Change 2007: The Physical Science Basis[END_REF], [START_REF] Knutti | Climate model genealogy: Generation CMIP5 and how we got there: CLIMATE MODEL GENEALOGY[END_REF]). Existing approaches mainly focus on seeking long-term trends, and the study of abnormal behaviour tend to focus on the search of extreme values.

An important element in climate data analysis is the observation window. For instance, a thirty years window has been commonly accepted for climate studies, and is now being reconsidered given the impact it can have on change detection (e.g. [START_REF] Sippel | Quantifying changes in climate variability and extremes: Pitfalls and their overcoming[END_REF]). Thus to our knowledge, in climate data analysis, we note 1) a lack of robust outlier pattern detection, and 2) a need to consider very large data sets to minimize the bias induced by the limited observation window. We propose to address those issues in this paper, by introducing a novel concept of pattern outlier, and evaluating our approach to the field of climate data.

In the realm of data mining, outlier detection is receiving much interest, and has shown its benefits in a wide range of applications including fraud detection [START_REF] Hilal | Financial Fraud: A Review of Anomaly Detection Techniques and Recent Advances[END_REF], cybersecurity [START_REF] Hopkins | Sub-Erroneous Outlier Detection of Cyber Attacks in a Smart Grid State Estimation System[END_REF] and the health sector [START_REF] Xia | Outlier knowledge management for extreme public health events: understanding public opinions about COVID-19 based on microblog data[END_REF]. Identifying outliers through data sets contributes towards decision support, risk and impact studies. Various definitions of what constitutes an outlier have been proposed, along with associated detection methods. A survey [START_REF] Ayadi | Outlier detection approaches for wireless sensor networks: A survey[END_REF] specifies twelve different interpretations of outliers from the perspective of different studies. Overall, an outlier can be commonly defined as "an observation that is significantly dissimilar to other data observations or an observation that does not behave like the expected typical behavior of the other observations" [START_REF] Hodge | A survey of outlier detection methodologies[END_REF]. The observations can specify a single point outlier, or a shape/pattern denoting an abnormal sub-sequence over a time series data. The latter form our outliers are also called discords.

In this article, we investigate discords over contextual time series. Furthermore, we are interested in very large data series to mitigate the potential impact of the length of the data set at hand upon the results. We apply our approach to large data series coming from monthly data between 1902 and 2005. A scalable and exact approach that has proven its computational and space efficiency to detect discords is the Matrix Profile method [START_REF] Yeh | Matrix profile I: all pairs similarity joins for time series: a unifying view that includes motifs, discords and shapelets[END_REF]. It requires the length of the sub-sequence (window) to be set as a parameter. The chosen window has a strong impact to detect meaningful discords. In existing studies, the detection algorithm is run with different window [START_REF] Madrid | Matrix profile xx: Finding and visualizing time series motifs of all lengths using the matrix profile[END_REF], leading potentially to multiple discords. The ranking of such discords is challenging since it requires meaningful criteria to prefer a discord to another. If the data is labelled such ranking is possible since it can be related to an event, but in case of unlabelled data it usually requires expert knowledge.

We propose a novel concept and approach to identify relevant discords over different windows automatically. To do so, we introduce the concept of prominent discord that specifies the most significant discord as an anomalous pattern over the longest continuous period, from a shared starting position. A core benefit of the prominent discord is to gain insight onto discords that coincide over their start date, while searching for the ones with longer and subsuming anomalous pattern. For instance if a drought is found through a window of four months, we seek whether it belongs to a longer dry period that may last six months, one year, ten years. By doing so, we search for the longest span of a discord that would subsume other discords and relates to similar occurrences. In contrast to point outliers that are likely to indicate extreme events, resulting from potential long term changes, prominent discords would reveal the anomalous patterns that cause such events.

We formally define the concept of prominent discord, propose a detection algorithm and present an application to large date sets of so called climate impact runoff data. This means that they are historical data relative to surface and subsurface runoff observations, a variable that provides information about flood and drought risks depending on values being high or low.

Figure 1 gives an intuition of what a prominent discord is. We show three discords of respective lengths 13, 37 and 58; all starting at the same position in the series. The reading of the plots corresponds to daily runoff data over five years in millimetres. The Xs covers the daily timestamp while the Ys the height of the runoff in the Sahelian region. Peaks clearly indicate the rainy season. For each window length the prominent discord is highlighted in a continuous line over the time spam. We notice that when the window size changes we have three prominent discords, all starting at the same position date and the one corresponding to the window size of 58 months covers the other two. It is the most prominent one.

The longest one discovered is the most prominent discord: it subsumes the other two. It actually gives insight into lasting changes and anomalous behaviour that could lead towards a change of system state altogether when studying global change and climate data.

The main contributions of this paper are: 1) the novel concept, formalization and method to compute prominent discords and extract the most prominent ones applied to large-scale time series, 2) its application to climate-related data to detect and evaluate the relevance and insights of such discords from a climate point of view.

The paper is organized as follows: section 2 gives a background and related work in the field, section 3 presents our conceptual and methodological contributions, and section 4 its evaluation before concluding in section 5.

BACKGROUND AND RELATED WORK

In this section we review previous work relative to time series discord discovery, more specifically with variable length discords and very large time series such as climate models data.

Climate data analysis

In climate data studies, long-term trends, defined as a tendency towards a climate change pattern, are often characterized by basic statistical measures, such as the average rate of variables increase/decrease over a given time period [START_REF] Hartmann | Observations: Atmosphere and Surface[END_REF]. In the field of weather extreme events (e.g. droughts, floods, heatwaves, storms), a common approach consists in quantifying how a given climate indicator jumps out, against the background of former climate records (in intensity, frequency or duration; ; [START_REF] Sillmann | Climate extremes indices in the CMIP5 multimodel ensemble: Part 2. Future climate projections: CMIP5 PROJECTIONS OF EXTREMES INDICES[END_REF]). These approaches are conducted under an arbitrary choice of a base time window since a climate reference is inherently defined to assess whether a long-term change and extreme event occurrences can be considered as emergent and/or anomalous. This choice greatly affects the meaning and the robustness of climate studies outcome when this reference is shifted [START_REF] Sippel | Quantifying changes in climate variability and extremes: Pitfalls and their overcoming[END_REF].

Time series discord discovery

Time series discord detection is receiving an increasing interest in data mining since its formalization ([13]- [START_REF] Kleijnen | Response surface methodology for constrained simulation optimization: An overview[END_REF]). Efficient and exact methods have been proposed to discover discords in data series [START_REF] Yeh | Matrix profile I: all pairs similarity joins for time series: a unifying view that includes motifs, discords and shapelets[END_REF], [START_REF] Yankov | Disk aware discord discovery: Finding unusual time series in terabyte sized datasets[END_REF]. These approaches require some parameter settings including the size of the observation window. The window size is fixed and needs to be specified as an input parameter (HOT SAX [START_REF] Keogh | Hot sax: Efficiently finding the most unusual time series subsequence[END_REF], QUICK MOTIF [START_REF] Li | Quick-motif: An efficient and scalable framework for exact motif discovery[END_REF]). As a result, recent works have drawn on the challenges and insight limitations of a fixed set window size leading to research towards computing all possible discords within a size range, using different methods such as quadratic regression [START_REF] Leng | Variable length methods for detecting anomaly patterns in time series[END_REF], dynamic time warping (DTW) [START_REF] Vy | Detecting variable length anomaly patterns in time series data[END_REF], or a graph-based approach [START_REF] Boniol | Series2graph: Graph-based subsequence anomaly detection for time series[END_REF].

PanMatrix [START_REF] Madrid | Matrix profile xx: Finding and visualizing time series motifs of all lengths using the matrix profile[END_REF] and VALMOD [START_REF] Linardi | Matrix profile goes MAD: variable-length motif and discord discovery in data series[END_REF] compute variable-length top k discords, using the Matrix Profile method for different window sizes, given a value k. VALMOD considers an interval of possible subsequence lengths as initial parameter; whereas PanMatrix computes exact distances for subsequences of all lengths.

These approaches address the issue of multiple discord computations, but there has been no attempt to order the discords of variable length, and seek those that would have more impact in terms of revealing a potential change of state. Also, the role played by the data series time span has received little attention with respect to its link to discord discovery. Both issues are important for climate model data towards insightful impact studies. A key element is to investigate coincident variable length discords to be able to extract actionable insight through the identification of prominent discords, longest ones sharing a starting position with smallest discords, and thus subsuming them. This is the goal of our approach that can be considered as a meta discord discovery problem over very large data series, with no a priori interpretation of outlier patterns. In other words, we seek discords for which all the subsequences within a resulting length interval are also discords sharing their starting position.

We use the Matrix profile approach, because it is an exact method to compute discords for a given window length. It is also computational and space efficiency. Let us now the key notions at hand that will be used to formalize our concept of prominent discord.

Matrix profile and discords

The matrix profile is a data structure computed to discover discords and motifs using similarity search algorithms [START_REF] Yeh | Matrix profile I: all pairs similarity joins for time series: a unifying view that includes motifs, discords and shapelets[END_REF]. Many algorithms have been proposed with different space and time performance (e.g. STOMP and GPU-STOMP [START_REF] Zhu | Matrix profile ii: Exploiting a novel algorithm and gpus to break the one hundred million barrier for time series motifs and joins[END_REF], SCRIMP and SCRIMP++ [START_REF] Zhu | Matrix profile XI: SCRIMP++: time series motif discovery at interactive speeds[END_REF]). The data structure builds on the following notions [START_REF] Zhu | Matrix profile ii: Exploiting a novel algorithm and gpus to break the one hundred million barrier for time series motifs and joins[END_REF], recalled hereafter for further usage and completeness : In other words, the discord of length m, is the subsequence in the data series (specified by its starting position j), such that among the shortest Euclidean distance, it is the one with the maximal value, ie. with largest anomalous pattern among all.

Definition 1. A time
The following example gives the Distance profiles for two susbsequences (in red) of window size 4 over a time series of length 13. The distance profile vector gives the z-Euclidian distance between the chosen subsequence and all the other ones (the next one is highlighted in blue). There are 10 of them, one per subsequence of size 4, sliding over the time series.

The resulting matrix profile extracts for each subsequence the smallest distance in each distance profile (yellow). Finally, we extract the discord from the matrix profile that is the subsequence corresponding to the greatest distance value in the profile. In this example, the biggest distance value is 14.1,distance between subsequence 7 and subsequence 3 (𝑉 7 = 3). Thus the discord is the seventh subsequence [START_REF] Chiu | Probabilistic discovery of time series motifs[END_REF][START_REF] Kleijnen | Response surface methodology for constrained simulation optimization: An overview[END_REF][START_REF] Randall | Climate models and their evaluation, Climate Change 2007: The Physical Science Basis[END_REF][START_REF] Knutti | Climate model genealogy: Generation CMIP5 and how we got there: CLIMATE MODEL GENEALOGY[END_REF] in the time series.

PROMINENT DISCORDS: DEFINITION AND ALGORITHMS

In this section we present our approach and contributions, that include the concept of prominent discord and the method we developed to extract a set of prominent discords. The overall problem we address is the following:

Problem addressed: Given a data series 𝑇, a large range of possible lengths and a unit step, compute all the variable length discords, and find the discords, such that all of their subsequences within a length interval and shared starting position in 𝑇 are discords. From these prominent discords, identify a relevant ordering that relates to their length, and number of subsumed discords.

Definitions

We work at a meta level with respect to the matrix profile and the discord since we compute matrix profiles over sequences of variable lengths, and seek the longest discords that contain discords with a given ordering relation. To formalize our problem, we introduce three concepts: discord profile, discord subsumption ordering and prominent discord. Note that to allow a reliable reasoning over a large number of subsequences, the interval for the variable lengths is set to [4, Our approach will compute a set of prominent discords for a given time series. We propose an ordering that accounts for 1) the number of subsumed discords (of different lengths of course), and 2) the relative length of the shortest subsumed discord. The idea behind this ordering is to exploit discords for insight studies on the anomalous patterns that can have a lasting impact, and pertained change of behavior in the time series. Intuitively, a point outlier can be the consequence of an existing change of behavior (e.g. rising number of extreme weather events), whereas a discord of droughts of length 4, also found in subsequences of lengths 8 and 15 for example, can indicate a first anomalous pattern, that pertains as an anomalous pattern in longer discords. The longest discord subsuming a much shorter one, can indicate a potential important weather change, and impact on soils, agriculture etc.

The ranking function sorts the prominent discords in decreasing order of 𝑠𝑜𝑟𝑡((𝑙 𝑒𝑛𝑑 -𝑙 𝑠𝑡𝑎𝑟𝑡 ) 𝑙 𝑠𝑡𝑎𝑟𝑡 ⁄ ), where the top value corresponds to the longest set of subsumed discords (𝑙 𝑒𝑛𝑑 -𝑙 𝑠𝑡𝑎𝑟𝑡 ), and the lower one the length of the first subsumed discord (𝑙 𝑠𝑡𝑎𝑟𝑡 ). As illustrated in Figure 3 the prominent discord A will outrank prominent discord B even though they subsume the same number of discords, because A builds upon a shorter outlier (𝑙 𝑠𝑡𝑎𝑟𝑡 ).

Prominent discord A Prominent discord B Figure 3. Prominent discord ordering: A is preferred to B with higher ratio function value

Algorithms

To derive the set of prominent discords, we first derive the discord profile (Algorithm 1) over variable length discords and extract the prominent ones (Algorithm 2) through a counting method based on shared starting positions. Note that Algorithm 1 makes use of an efficient matrix profile computation algorithm (line 5), the STOMP algorithm [START_REF] Zhu | Matrix profile ii: Exploiting a novel algorithm and gpus to break the one hundred million barrier for time series motifs and joins[END_REF] omitted for space reasons. This algorithm, like other matrix profile computation methods (eg, STAMP) derives the z-normalized Euclidean distance to measure efficiently the distance between subsequences.

Algorithm 1 takes as input the whole time series, a maximum subsequence length and list of variable lengths (line 2--3), called windows (from the terminology of the matrix profile approach) that specifies the subsequence lengths considered. For each window size (lines 4--7), we compute the matrix profile, extract the discord 𝛥 𝑗,𝑙 to be stored in the Discord profile list 𝛥𝑃.

The main algorithm, Algorithm 2, computes the list of prominent discords 𝛥𝐶 and returns the sorted list of prominent discords (in decreasing value of respective (𝑐𝑜𝑢𝑛𝑡 𝑙 𝑠𝑡𝑎𝑟𝑡 ⁄ ) value. It takes as input the time series and returns the sorted list of prominent discords 𝛥 ‾ 𝑗,𝑙 𝑠𝑡𝑎𝑟𝑡 ,𝑙 𝑒𝑛𝑑 including its starting position, first discord length and longest one. Line 3 initializes a counter of discords having identical starting positions. In line 4 the discord profile 𝛥𝑃 is computed from Algorithm 1 and contains all the discords 𝛥 𝑗𝑙 one per length 𝑙 considered in Algorithm 1. Line 5 and 6 define the variables used to extract the length 𝑙 and starting position 𝑗 of a discord in 𝛥𝑃. Line 7 extracts the length of the first discord. Lines 9--11 increment the count as long as the next discord has the same starting position as the current one denoted by 𝑖. Lines 12--14 create a new prominent discord with starting position 𝑗, starting length 𝑠, last length 𝑠 + 𝑐𝑜𝑢𝑛𝑡. It is added to the prominent discord list. Lines 15--16 reinitialize the count, and new starting 𝑠 position to the one of the next discord in 𝛥𝑃. Line 19 sorts the prominent discord list in decreasing order according to the proposed function (𝑙 𝑒𝑛𝑑 -𝑙 𝑠𝑡𝑎𝑟𝑡 )/𝑙 𝑠𝑡𝑎𝑟𝑡 ; and finally line 20 returns the final sorted list 𝛥𝐶.

Worst case time complexity. Algorithm 2 (calling algorithm 1) overall runs in the worst case in 𝑂(𝑛 3 ) where 𝑛 is the length of the time series. To decompose, we have: Algorithm 1 calls 𝑛/2 times STOMP, thus runs in the worst case in 𝑂(𝑛 2 × 𝑛) = 𝑂(𝑛 3 ). The For loop in Algorithm 2 runs in the length of the discord profile list thus in the worst case 𝑂(𝑛) since there is one discord per length (𝑛/2 variable length discords), and the list of prominent discords is sorted in the worst case in 𝑂(𝑛𝑙𝑜𝑔𝑛).

EXPERIMENTAL EVALUATION AND COMPARISONS

We now present an application of our method to the analysis of large datasets relative to runoff historical climate data. Runoff data correspond to measures of waters in terms of distance (in mm) above the land surface to reach a stream but also to infiltrate the soil surface. All experiments were run on an Intel(R) Xeon(R) Bronze 3106 CPU processor at 1.70GHz with 8 core with 64 GB of RAM. We also compare our proposed approach to other discord discovery methods.

The climate data. We consider observed monthly runoff data, defined in climate data science as an impact variable analyzed to quantify flood and drought risks at regional and global scales (e.g. [START_REF] Sillmann | Climate extremes indices in the CMIP5 multimodel ensemble: Part 2. Future climate projections: CMIP5 PROJECTIONS OF EXTREMES INDICES[END_REF], [START_REF] Arnell | The global-scale impacts of climate change on water resources and flooding under new climate and socio-economic scenarios[END_REF], [START_REF] Gosling | A comparison of changes in river runoff from multiple global and catchment-scale hydrological models under global warming scenarios of 1 °C, 2 °C and 3 °C[END_REF]). These monthly runoff observations are obtained from the Global Runoff Reconstruction dataset (GRUN) that covers the 1902-2014 time period (113 years), with a 0.5° × 0.5° spatial grid resolution [START_REF] Ghiggi | GRUN: an observation-based global gridded runoff dataset from 1902 to 2014[END_REF]. We focus our analysis on the Sahel region, a particularly soil water vulnerable area, and we spatially average monthly runoff over the corresponding grid box [5°W-25°E ; 10°N-18°N]. Our prominent discord approach is then applied to the obtained Sahel time series between 1902 and 2005 (i.e. 104 years, 1248 months), a period commonly considered in historical climate analysis.

Prominent discords discovery

For the dataset of monthly runoff observation data over 1248 months, we applied our approach and derived the list of all prominent discords, including their ranking based on our proposed ratio function, after calculating the discords for all variable length windows in the interval [4,..,1248/2], with a monthly step increment. The set of prominent discords was derived in 3.5 minutes.

Figure 4 shows for five prominent discords, including the most prominent one, their respective subsumed discords. The X-axis indicates window lengths up to 130 months, and the Y-axis the discords starting date. Each blue dot represents a discord with its window length (Xs) and its starting date (Ys). The arrows show the prominent discords with all the subsumed discords of coinciding starting dates. The length of the arrows illustrates how many discords the prominent one subsumes. Here we have four prominent discords. The most prominent discord starts at the position date 1903-10-15, with a lower window size of 13 months for its first subsumed discord, and upper length of 58 months.

Comparison with alternative approaches

We compared our approach with existing discord discovery methods, illustrating mainly the importance of considering an exact approach for the prominent discord extraction, and the need of variable length discord computation without parameterized length settings.

We considered HOT-SAX, an extension of the SAX algorithm [START_REF] Lin | A symbolic representation of time series, with implications for streaming algorithms[END_REF]. SAX discretizes time series into words and detects motifs in time series but not discords. HOT-SAX algorithm was developed to detect discords. It builds a suffix tree that stores the words generated by SAX. A word with the least number of occurrences is a discord. A requirement is that the number of extracted discords is predefined.

Rare Rule Anomaly (RRA) [START_REF] Senin | Time series anomaly discovery with grammar-based compression[END_REF] is an algorithm that uses grammar-based compression able to detect motifs and discords in time series. Similar to HOT-SAX, RRA uses SAX algorithm to discretize the time series. A grammatical induction algorithm (ex: Sequitur [START_REF] Nevill-Manning | Identifying hierarchical structure in sequences: A linear-time algorithm[END_REF]) is used to generate the grammar. These generated grammars are used to detect the discords.

One of the main parameters for HOT-SAX and RRA is the window size. To be able to compare them with our approach, we use the value 𝑙 𝑒𝑛𝑑 of each prominent discord as the window parameter. We can see that both HOT SAX and RRA lead to different results in terms of starting and end date of the prominent discord for a given window size. This comes from the fact that they are not exact methods and given an input window length, lead to a different discord.

We also compared with the PAN MATRIX algorithm. It calculates variable length discords, using the SKIMP algorithm to compute the matrix profile for all motif lengths.

Analysis and insights for climate data analysis

To analyze those results, as well as the relevance and insights of our prominent discord and proposed orderings, we compared with statistical analysis of those historical climate data (e.g. long-term trends, seasonal cycles, standard-deviation). It is worth noting that the prominent discord approach is unsupervised, and does not consider any climate behaviour nor known physical processes.

Figure 6 shows the annual average runoff data over the Sahel region illustrating the maximum, minimum and average annual values, providing a general idea of the global fluctuations and extremes. According to climate studies, a well known soil drying trend mostly resulting from a rainfall decrease is observed between 1900 and 2013 in the Sahel ((e.g. [START_REF] Ghiggi | GRUN: an observation-based global gridded runoff dataset from 1902 to 2014[END_REF])). We find consistent trends using runoff data, particularly in the annual maximum runoff time evolution between 1902 and 2005 ( [START_REF] Gosling | A comparison of changes in river runoff from multiple global and catchment-scale hydrological models under global warming scenarios of 1 °C, 2 °C and 3 °C[END_REF]). In our work, four of the five first prominent discords coincide with starting dates and length intervals during the first 22 years of the time series (Table 5 and Figure 7). The corresponding averaged monthly runoff over 1902-1924 also shows a mean of 15e -07 kg/m 2 /s and a standard-deviation (i.e. an interannual variability indicator; ([31]) of 0.14 kg/m 2 /s, whereas the entire 1902-2005 time series is characterized by a mean of 12e -07 kg/m 2 /s and a standard deviation of 0.12 kg/m 2 /s. These four prominent discords may thus illustrate the specific time pattern between 1902 and 1924 with higher soil water amount and larger inter-annual variability, before the upcoming continuous long-term drying trend observed within Sahel.

The second prominent discord is detected for a starting date at 1982-09-15 and window sizes in [109 ; 120 months]. This period corresponds to the time period with the smallest mean runoff values of the entire time series (Fig. 6 and Fig. 7). The associated 1982-1992 mean runoff is 9e -07 kg/m 2 /s compared to a mean runoff of 12e -07 kg/m 2 /s over 1900-2005. This prominent discord illustrates the temporal pattern resulting from intense droughts that occurred in Sahel in the 1980s. During that decade, the most severe drought ever recorded over the African continent occurred during 1983-1984 (Figure 2 in [START_REF] Masih | A review of droughts on the African continent: a geospatial and long-term perspective[END_REF]).

In this study, we demonstrated two major usefulness of our proposed prominent discord approach and ordering relations, in climate data analysis in terms of real insights towards the emergence of a longterm change, and the detection of recurring anomalous events.

First, we showed that the prominent discord concept does capture a 20-30 years pattern illustrating a different (former) climate regime compared to the rest of the considered time series (emergence).

Second, we showed that the subsumption ordering and prominent discords ranking capture the time pattern at a decade scale of the driest recorded yearly event (recurring anomaly detection).

CONCLUSION AND FUTURE WORK

In this paper we proposed a new concept in the realm of variable length discords, the prominent discord. It focuses on identifying anomalous patterns that last through time and subsume sets of discords. The main contributions are the formalization and method to compute prominent discords and extract the most prominent ones applied to large scale time series, and the application to climaterelated data.

Our ordering and results show their relevance in the field of climate data series. In particular, they show that through such ordering we gain insights on the anomalous patterns that have a lasting impact, and pertained change of behavior in the time series. This is new to our knowledge.

Future works include further experimental studies on different impact climate data and other data sets, to evaluate the thematic insights of our approach, as well as some optimization of the algorithm to ensure scalability.

Figure 1 .

 1 Figure 1. Prominent discord and subsequent discords

Figure 2 .

 2 Figure 2. Distance and matrix profile, sequences of length 4

Fig. 4 .

 4 Fig.4. Prominent discords derived from observed runoff monthly data.

Figure 7

 7 Figure6shows the annual average runoff data over the Sahel region illustrating the maximum, minimum and average annual values, providing a general idea of the global fluctuations and extremes. Figure7relates the monthly runoff values together with the top five prominent discords we discovered, cf Figure5. For each color, the dotted lines indicate the starting position of a prominent discord and the vertical lines the respective lstart and lend.

Figure 6 .Figure 7 .

 67 Figure 6. Annual mean runoff values

  series 𝑇 is a sequence of real-valued numbers 𝑡 𝑖 : 𝑇 = [𝑡 1 , 𝑡 2 , . . . , 𝑡 𝑛 ] where n is the length of 𝑇. A subsequence 𝑇 𝑖,𝑚 of a time series 𝑇 is a continuous subset of values in 𝑇, of length 𝑚 and starting at position 𝑖. Formally, 𝑇 𝑖,𝑚 = [𝑡 𝑖 , 𝑡 𝑖+1 , … , 𝑡 𝑖+𝑚-1 ], where 1 ≤ 𝑖 ≤ 𝑛 -𝑚 + 1. A distance profile 𝐷 𝑖,𝑚 of time series 𝑇 and length 𝑚 is a vector of the z-Euclidean distances between a given query subsequence 𝑇 𝑖,𝑚 and all subsequences of length 𝑚 in the time series 𝑇. Formally, 𝐷 𝑖,𝑚 = [𝑑 𝑖,1 , 𝑑 𝑖,2 , … , 𝑑 𝑖,𝑛-𝑚+1 ], where 𝑑 𝑖,𝑗 (1 ≤ 𝑖, 𝑗 ≤ 𝑛 -𝑚 + 1) is the distance between 𝑇 𝑖,𝑚 and 𝑇 𝑗,𝑚 with 𝑖 ≠ 𝑗. The z-normalized Euclidean distance 𝑑 𝑖,𝑗 between two subsequences 𝑇 𝑖,𝑚 and 𝑇 𝑗,𝑚 of length m, is defined by: where 𝜇 𝑖 and 𝜎 𝑖 are respectively the mean and standard deviation of 𝑇 𝑖,𝑚 , 𝜇 𝑗 and 𝜎 𝑗 the mean and standard deviation of 𝑇 𝑗,𝑚 . A matrix profile 𝑃 𝑚 of time series 𝑇 and given length 𝑚 is a meta series of the Euclidean distances vector between each subsequence 𝑇 𝑖,𝑚 of given length 𝑚 where 𝑖 varies, and its nearest neighbor (closest match) in time series 𝑇, together with the corresponding position vector for each closest neighbor associated with 𝑚𝑖𝑛(𝐷 𝑖,𝑚 ) . We denote it 𝑃 𝑚 = [𝑚𝑖𝑛(𝐷 1,𝑚 ), . . . , 𝑚𝑖𝑛(𝐷 𝑛-𝑚+1,𝑚 )], where 𝐷 𝑖,𝑚 (1 ≤ 𝑖 ≤ 𝑛 -𝑚 + 1) is the distance profile 𝐷 𝑖,𝑚 of time series 𝑇 for subsequences of length 𝑚.Pm = [min(D1,m),...,min(Dn-m+1,m)], where Di,m(1 ≤ i ≤ n-m+1) is the distance profile Di,m of time series T for subsequences of length m. A matrix profile index vector 𝑉 𝑚 , associated with a matrix profile 𝑃 𝑚 denotes the vector of starting position 𝑗 of the subsequence corresponding to the minimal distance. It is specified by the vector: 𝑉 𝑚 = [𝑉 1 , 𝑉 2 , … , 𝑉 𝑛-𝑚+1 ], such that 𝑉 𝑖 = 𝑗 if 𝑑 𝑖,𝑗 = 𝑚𝑖𝑛(𝐷 𝑖,𝑚 ).

	Definition 2. Definition 3 (Distance Profile). 𝑑 𝑖,𝑗 = √2𝑚 (1 -Definition 4 (Matrix Profile). Definition 5 (Index vector).	𝑇 𝑖,𝑚 . 𝑇 𝑖,𝑚 -𝑚𝜇 𝑖 𝜇 𝑗 𝑚𝜎 𝑖 𝜎 𝑗	)	(1)

Definition 4 (z-Euclidian distance).

Definition 6 (Discord).

A discord denoted 𝛥 𝑗,𝑚 is a subsequence 𝑇 𝑗,𝑚 of length 𝑚 starting at the position 𝑗 in 𝑉 𝑚 , that corresponds to the maximum distance value in 𝑃 𝑚 .

  Note that this ordering relation is a partial order since we assume that two discords with different starting position in the time series are not comparable. The motivation behind this ordering is that such discords might relate to very different events, whereas discords that co-occur in their starting position are more likely to share the root event for the outlier pattern. With identical start position, two discords can relate to the same anomalous pattern. This is not necessarily true for different starting positions. A Prominent discord 𝛥 ‾ 𝑗,𝑙 𝑠𝑡𝑎𝑟𝑡 ,𝑙 𝑒𝑛𝑑 of a time series 𝑇 is the top discord 𝛥 𝑗,𝑙 𝑒𝑛𝑑 of a lattice of all discords 𝛥 𝑗,𝑙 in 𝛥𝑃 such that 𝛥 𝑗,𝑙 𝑠𝑡𝑎𝑟𝑡 ≼ 𝛥 𝛥 𝑗,𝑙 ≼ 𝛥 𝛥 𝑗,𝑙 𝑒𝑛𝑑 .

	Definition 9 (Prominent discord).

..,n/2]. Definition 7 (Discord profile). The discord profile 𝛥𝑃 of a time series 𝑇 is a set of discords 𝛥 𝑗,𝑙 of variable lengths 𝑙 and starting positions 𝑗, such that 𝛥𝑃 = {𝛥 𝑗,𝑙 | 𝑗 ∈ 1. . 𝑛/2 -1,0 < 𝑙 ≤ 𝑛/2}. Definition 8 (Discord subsumption ordering). A discord 𝛥 𝑗,𝑚 subsumes a discord 𝛥 𝑖,𝑙 , specified as 𝛥 𝑖,𝑙 ≼ 𝛥 𝛥 𝑗,𝑚 if and only if, 𝑖 = 𝑗 and 𝑙 < 𝑚..

Table 1

 1 shows the top five prominent discords in ∆C according to our proposed ratio ordering. The first and second ones are found in Figure4.

	date	starting window size	ending window size ratio
	1903-10-15	13	58	3.46
	1982-09-15	109	120	0.10
	1902-09-15	242	264	0.09
	1903-09-15	193	209	0.08
	1902-01-15	80	86	0.075

Table 1 .

 1 Top 5 prominent discords sorted by the ratio function

Table 2 .

 2 HOT SAX results using windows extracted from the 𝑙 𝑒𝑛𝑑 of each prominent discord

		HOT SAX
	window size	start	end
	58	1902-01-01	1906-11-01
	120	1943-08-01	1953-08-01
	264	1972-02-01	1994-02-01
	209	1973-06-01	1990-11-01
	86	1902-07-01	1909-09-01
		RRA
	window size	start	end
	58	1902-02-01	1907-04-01
	120	1944-08-01	1954-11-01
	264	1953-01-01	1977-08-01
	209	1953-12-01	1971-12-01
	86	1943-04-01	1951-01-01

Table 3 .

 3 RRA results using windows extracted from the 𝑙 𝑒𝑛𝑑 of each prominent discord

Table 4

 4 shows the top five discords with the PAN MATRIX. Compared with our results in Table1, the top 5 discords are different, and the top 5 discord of PAN MATRIX are not relative to subsequences since they are not based on an ordering among variable length discords. PAN MATRIX is efficient to calculate variable length discord, but is not designed to order discords.

	Pan Matrix
	Date	Window size
	1952-07-01	619
	1952-07-01	618
	1952-07-01	617
	1952-07-01	616
	1952-07-01	615

Table 4 .

 4 PAN MATRIX top 5 discords

ACKNOWLEDGEMENTS

The authors would like to thank the Occitanie Region, who partially funded this research, and the readers of previous versions for the tremendous guidance.