Are Tumor Cell Lineages Solely Shaped by Mechanical Forces?

This paper investigates cell proliferation dynamics in small tumor cell aggregates using an individual-based model (IBM). The simulation model is designed to study the morphology of the cell population and of the cell lineages as well as the impact of the orientation of the division plane on this morphology. Our IBM model is based on the hypothesis that cells are incompressible objects that grow in size and divide once a threshold size is reached, and that newly born cell adhere to the existing cell cluster. We performed comparisons between the simulation model and experimental data by using several statistical indicators. The results suggest that the emergence of particular morphologies can be explained by simple mechanical interactions. Electronic supplementary material The online version of this article (doi:10.1007/s11538-017-0333-y) contains supplementary material, which is available to authorized users.


Introduction
Cancer cells proliferate at a high rate and can be considered as a dynamic population of agents that grow and divide without constraints (Hanahan and Weinberg 2000) (at least in the early phase of avascular growth). In the present work, we aim to investigate this preliminary stage of a tumor growth. As such, this study could also help in understanding the development of 3D microtumours. Here, we examine a population initially composed of several tens of cells that proliferate to reach several hundreds of cells within a few days. In this situation, the experimental model consists of cancer cells grown in a culture medium containing all the necessary nutrients for their growth and division. The proliferation of the population is restricted to a single layer, which permits the use of simple two-dimensional simulation models for the comparisons. In the proposed simulation model, we do not include a detailed description of the cell cycle and mitosis events. Instead we focus on the role of the orientation of the cell division planes in the morphology of the tumor cells cluster. We aim to investigate how the orientation of cell divisions influences the structure of the cell lineages. In particular, we would like to determine whether cell lineages break up and whether they have different morphologies according to their initial position in the cell population. The impact of the orientation of the division plane on the organization of the population has been suggested from recent studies. In particular, the influence of the geometry of the cell, the influence of neighboring cells and the role of external mechanical forces on the determination of the orientation of the division plane have been studied in Gibson et al. (2011), Minc et al. (2011, Théry and Bornens (2006), Desmaison et al. (2013a), Fink et al. (2011), Thry et al. (2007), Nestor-Bergmann et al. (2014) and Ambrosi and Mollica (2004). In Höhme et al. (2010), the authors studied a model in which mitotic spindle has random orientation and they compare it with a model in which division is subsequent of an attraction mechanism. Recent results also suggest that the orientation of the division plane plays a role in the differentiation of stem cells (Lechler and Fuchs 2005;Wilcock et al. 2007). It has been experimentally observed that oriented cell division can generate cellular diversity and that tissue morphogenesis depends on the control of the direction of cell division. In other words, orientation shapes tissues. Moreover, recent discoveries tend to establish a link between cancer and disorientation of the division plane (Pease and Tirnauer 2011): dysfunction in cells can lead to disorientation and, conversely, disorientation can promote the development of cancer. The above observations lead us to focus on two main questions related to the organization of a growing tumor cell population: what is the impact of the multiplication of cells on the global organization of the entire cell population? Does the orientation of the division plane influence the evolution of the cell population and the organization of its lineages? The comparison of a mathematical model with biological experiments performed in this work shows that indeed a relation exists between division and organization and that lineages are strongly influenced by the initial position of the parental cell inside the population.
There exists a large number of mathematical models in the literature describing cell proliferation and tumor growth. Some recent review papers can be found in Araujo and McElwain (2004), Byrne and Drasdo (2010), Byrne et al. (2006) and Roose et al. (2007). Among the different ways of representing cells, one usually distinguishes between discrete models and continuous models. In a discrete model (Kam et al. 2012;Rejniak and McCawley 2010;Beyer and Meyer-Hermann 2007;Bonabeau 2002;Byrne 2009;Montalenti et al. 1998;Chiorino et al. 2001;Ingham-Dempster et al. 2017;Zhang et al. 2016;Jagiella et al. 2016;Fletcher et al. 2012;Mirams et al. 2012;Kurusawe et al. 2015), each element is treated as a separate entity. This makes the comparison with experiments easy, but the main drawback of this viewpoint is the huge computational cost when dealing with a large number of agents. On the other hand, continuous models (Johnston et al. 2007;Byrne and Preziosi 2003;Byrne et al. 2006;Araujo and McElwain 2006;Bresch et al. 2010;Sherratt and Chaplain 2001), which typically deal with an average density of cells, are more efficient when the system contains a large number of particles. However, it is difficult to establish direct links between the model parameters and the physical measures (Macklin et al. 2009). Comparisons between the two approaches can be found in Schaller and Meyer-Hermann (2006), Byrne and Drasdo (2009), Galle et al. (2009a), while hybrid models which employ both approaches at the same time can be found for instance in Ribba et al. (2006), Zhang et al. (2007), Patel (2001), Di Costanzo et al. (2015) and Di Costanzo et al. (2017).
Many models of the literature tend to be exhaustive in the description of the biological, physical and chemical phenomena. This leads to the introduction of many empirical parameters and makes the interpretation of the results difficult. Our approach is opposite: it relies on a simple mathematical model which focuses on few determinants and attempts to explore some specific questions through comparisons with experiments. This approach permits to explore the influence of any single modeling choice more efficiently and to formulate hypotheses about which mechanisms are associated with given observations. We follow a bottom-up strategy which starts from simple rules and gradually adds complexity into the model until a good fit with the experiments is reached. Rather than quantitative agreement, we look for similar trends between the model and the experiments when some key parameters are varied. Here, we show that the sole growth and division mechanisms are not sufficient to explain the observed lineage morphologies and that additional phenomena must be taken into account in the model in order to reproduce the experimental results. The biological situation we wish to investigate is a small population composed of 20 up to about 500 cells in which proliferation and movement are restricted to a two-dimensional plane. For this specific situation, the best modeling choice is an individual-based models (IBM) (see Drasdo and Höhme 2005;Galle et al. 2009b;Meineke et al. 2001;Palsson 2008;Rejniak 2007 for example). This permits to make the mathematical model reproduce the experiment and to consider objects moving only in a two-dimensional plane. In addition, since we seek to study the impact of the growth of cells and the influence of the orientation of division on the organization of the population, we need to be able to track individual entities. Indeed, a continuous model will not give access to such information. On the other hand, a discrete model on a grid (Anderson et al. 2007;Glazier and Graner 1993;Shirinifard et al. 2009) could introduce artificial bias on the organization of the cells and it will not allow us to explore orientation issues in depth, since the number of possible orientations is limited by the underlying grid (Byrne et al. 2006;Macklin et al. 2009;Deutsch and Dormann 2005).
The mathematical setting chosen is finally the following: cells aggregate spontaneously, grow and divide. After each growth or division event, a mechanical equilibrium between aggregation and cell-cell non-overlapping is supposed instantaneously reached and gives the instantaneous configuration of the population. Thus, motion arises from stresses between neighboring cells (Tambe et al. 2011;Chepizhko et al. 2016;Levine et al. 2001;Shao et al. 2012). This approach is different from more classical models based on introducing a repulsion potential between the cells (Höhme et al. 2010;Byrne and Drasdo 2009;Drasdo and Höhme 2005). Indeed, the temporal scale associated with cell (quasi)-incompressibility is much faster than that involved in the growth of the tumor as a whole. Modeling cell-cell non-overlapping via a repulsion potential requires making these two scales closer than they are in reality (to ensure numerical stability), thereby introducing a bias in the numerical solution. In a forthcoming work (Degond et al. 2017), we intend to document precisely the differences between these two approaches. Since cell-cell non-overlapping is associated with a faster scale than growth an approach based on realizing a mechanical equilibrium at every time step permit to bypass the numerical stability issue. During the time evolution of the system, different lineages are tagged and compared with experimental data as done for instance in Seplveda et al. (2013). Comparisons between the mathematical model and the biological experiments show that a relation exists between geometric determinants of cell, division, and the organization of the cell population and that lineage shapes are strongly influenced by the initial position of the parental cell inside the population.
The paper is organized as follows. In Sect. 2, we discuss the mathematical model and the numerical method adopted. We also discuss the statistical indicators used to measure the results of the numerical simulations, the experimental protocol and the image processing. In Sect. 3, we resume the principal results of the simulations and a first series of comparisons between the model and the experiments are presented and discussed. In Sect. 4, improvements introduced to the model and new comparisons with the data are analyzed. In Sect. 5, conclusions are drawn and future investigations are discussed.
Only the plasma membrane is described. The details of the intracellular phenomena are omitted. The cell shape is chosen to be a two-dimensional disk which is incompressible and continuously growing in time. The notion of preferred orientation, linked to the alignment of chromosomes during mitosis, necessary to define the division plane, is not inherently modeled by the shape of the agent. Instead it is an internal parameter owned by each agent which may vary with time depending on the chosen division strategy and detailed next. The cell evolution is determined only by the growth and the division laws, respectively, describing the interphase and the mitotic steps of a cell cycle. To reproduce the experimental setting, cells are free to move and have access to all the required nutrients. Cell cycle phases of growth and division alternate continuously; cell cycle checkpoints, G1 phase variability between cells and temporary cell cycle exits are excluded. The interphase is constituted solely of the cell growth phase: G1 , S and G2 phases are combined into a single no-division phase during which growth is assumed to be linear in time. The mitosis trigger occurs as soon as the cell reaches a critical size, which is the only control condition. Conservation of the volume is imposed during the division process. Agents interact by minimizing at each time the global mechanical energy of the system subject to a non-interpenetration constraint modeling the fact that living cells cannot intermingle. This dynamic causes global as well as individual cell movement, which is thus the product of the combined actions of growth and division on the one hand, and the non-overlapping constraint on the other hand. Growth is modeled as a continuous phenomenon except at the time of division, which occurs when a cell reaches approximately twice the volume of a newborn daughter cell. A uniform probability distribution is added to the growth increment over a time step to introduce some randomness in the cell division starting time. When this process starts the mother cell deforms itself in a dumbbell-shaped geometry to give birth to two identical daughter cells. Deformation occurs with total volume kept constant. The duration of mitosis is short compared to the interphase (around one over thirty units of time). Thus, when a division occurs, the other cells stop growing, i.e., we consider mitosis as an instantaneous phenomenon. During division an equilibrium between the mechanical adhesion forces and the non-overlapping constraint determines the state of the system at each instant of time. In order to explore the influence of the orientation of the division plane onto the lineage organization, we consider three different possibilities: divisions occur in (1) a random direction, (2) in the direction of the line joining the origin and the mother cell (radial direction) (3) in the direction orthogonal to the radial direction (tangential direction). Note that we only consider the case where all cells make the same orientation choices for division (i.e., either all choose random, or all choose radial, or all choose tangential direction). We leave the investigation of the more complex case where some cells choose a certain orientation and others a different one for future investigations. Furthermore, in addition to the different division planes, we consider two different strategies. (1) Free orientation strategy: orientation is chosen at the beginning of the division but is free to change during the deformation into a dumbbell-like shape and finally into two daughter cells. This change of orientation is only due to the interactions with neighboring agents through the energy minimization procedure which defines the new configurations. (2) Constrained orientation strategy: Orientation is chosen at the beginning of the division and remains fixed up to the end of the division process. The first strategy, freedom for the orientation plane to change during division, models the situation in which the orientation of the plane of division is a consequence of the instantaneous state of stress at which cells are subjected. The second strategy, fixed orientation during the division, models the fact that in some situations it has been proved that cells are capable to detect their position and consequently to choose their orientation (Lander 2013). At each time step, a minimum of this mechanical energy subject to the non-overlapping constraint is computed. At the beginning of the next time step, cells radii grow and divisions may arise. This induces a disruption of the mechanical equilibrium, and thus a new minimum of the energy is computed. The time unit is chosen in such a way that the mean duration of a cell cycle is 24 time units which corresponds to the 24 h of the real mean duration of a cycle for the type of cells experimentally considered. The final time of the simulation is 72 time units which corresponds to an average of three cell cycles. The time step between two new equilibrium configurations is 0.25 time units. This leads to a radius growth of 1/100 in one time step with respect to its initial size. This means that the change in configuration of the cellular aggregate changes very slowly in one time unit. Figure 1c, d shows a typical result of the simulation of the described model. To ease the interpretation of the results, the spherical cell shape is replaced by a polygonal one following a Voronoï tessellation of the cell centers, detailed in Sect. 2.4. This simplifies the definition of the geometrical characteristics of the cellular aggregates and facilitates the localization of the neighboring cells. Figure 1a, b illustrates the initial and the final phases of the proliferation and monitoring of the lineages. The cells tracking procedure is described in Sect. 3.

Detail on the Model
Rules for Cells Positioning As already stated, each cell is described by a 2D incompressible disk with a center positioned at a radius R i (t) > 0 and an orientation ω i (t) ∈ S 1 (the set of two-dimensional vectors with unit length) depending on time t. In this setting, we use ξ(t) to denote the vector whose elements are the positions of the cells, i.e., ξ(t) = (X 1 (t), X 2 (t), . . . , X N (t) (t)), while ρ(t) is the vector whose elements are the radii of the cells, i.e., ρ(t) = (R 1 (t), R 2 (t), . . . , R N (t) (t)). The number of cells at time t is denoted by N (t). Each cell belongs to a lineage i which defines the developmental history of a given initial mother cell and which does not evolve with time. What evolves in time is the number of cells N i (t), belonging to a given lineage i , due to mitosis.
The impenetrability condition between two cells i and j is expressed by an inequality constraint φ i j with a suitable function φ i j which expresses the fact that two cells should not overlap. Thus, an admissible configuration A(t) for the system is a set of positions ξ(t) such that φ i j (ξ(t), ρ(t)) ≤ 0 for all possible indices i and j: The global adhesion potential is expressed by a function is a convex function on R 2 . The instantaneous configuration at time t is then given by a minimum ξ * (t) of the potential W under the constraint that ξ * (t) belongs to the set of admissible configurations A(t), i.e., In this setting, the non-overlapping condition is defined by is the Euclidean distance on R 2 between cell located at X i (t) and cell located at X j (t).
The potential function models the trend of the cells to regroup themselves isotropically around a given position chosen to be the origin of the coordinate system. The potential τ 0 , Mother cell τ 1 τ f −1 τ f , Daughter cells Fig. 2 Different steps of the division process. From one mother cell on the left, at initial time t = τ 0 , up to two daughter cells on the right at the end of the process t = τ f function we consider is quadratic . This choice models the situation in which at the center of the cellular aggregate cells are more necrotic and consequently they are probably less able to build up connection bridges with their neighbors.
Growth Law We introduce the size of a new born cell R min , the size of a cell just before mitosis R max and T G the mean duration of the growth phase. Even though the model is two-dimensional, we consider cells as tridimensional structures whose volumes grow linearly in time. Thus the growth law of the i-th cell is given by where γ is a random variable sampled from an uniform distribution with support on [−α, α]. Once a cell reaches a radius R(t) ≥ R max , it starts to divide into two daughter cells. Equation (4) is discretized in small time steps t. After a time step cell growth leads to the violation of non-overlapping constraints. Thus a new energy minimum must be computed through (3) resulting in a repositioning of the cells. Then a new growth step is performed followed by a repositioning step. The cycle of growth and repositioning is repeated until one cell starts to divide.

Division Rules
The initial orientation ω i 0 of the division plane of the cell C i is random, radial or tangential. The radial and tangential directions are computed relative to the origin supposed to be the center of the tumor. The division process starts when a cell C i reaches a size R i 0 (t) ≥ R max at time t. The process is considered as discrete in time and at each time step the disk which describes the mother cell stretches apart in a peanut like shape until the final separation in two daughter cells as shown in Fig. 2. During this process, the volume is kept constant equal to the volume of the mother cell. At each discrete instant of time, τ k , k ∈ [1, f ] (where f is the total number of intermediate steps in the division process) a new equilibrium of the whole system is computed by solving (3) with a modified set of admissible configurations A(τ k ) (at step τ k ) as described below. When the mitosis comes to an end, each daughter cell has reached half the size in volume of the mother cell. Moreover, they have the same size and shape and their position is symmetric with respect to the division plane. The orientation of the division plane is described by the unit vector ω i (τ k ), ∀i. Since the division process is much faster than the cell cycle, we make the hypothesis that two cells cannot divide at the same time and that during division the other cells do not grow.
Let ω i 0 (t) the orientation of the division plane of the mother cell when the division starts, (x i 0 (t), y i 0 (t)) its coordinates and i 0 its lineage. Then, initially, the two daughter cells occupy the same location in space as the mother cell, i.e., , y i 0 (t)); they share the same orientation ω + (τ 0 ) = ω − (τ 0 ) = ω i 0 (t) and they belong to the same lineage + (τ 0 ) = − (τ 0 ) = i 0 , where the upper indices + and − refer to the two daughter cells and the i 0 index to the mother cell. During division the lineage of the two daughters remains unchanged while the two radii R + (τ ) and R − (τ ) are functions of the time during division τ (which is rather a degree of completion of the division process), and are such that the initial volume of the mother cell is preserved in time. During the division process the real time variable t is kept constant. In particular, at the end of the process the two radii are such that which expresses the conservation of volume during the division process since the volume of a daughter cell is the volume of the mother cell at time t before the division starts. This value then defines the new positions through since the two new born cells are placed along the normal vector direction to the plane of division at the distance R ± (τ k ) −h(τ k ) from this plane. Once the new positions are computed, the non-overlapping constraint is likely to be violated. A new minimal energy configuration ξ * (τ k+1 ) must be computed at step τ k+1 solving (3). Here the definition of the set of admissible configuration is different from (1) and incorporates equality constraintsφ i j associated with the maintenance of the peanut shape when the pair (i, j) corresponds to two daughter cells, i.e., (i, j) = (i + , i − ). In addition in the case of fixed orientation strategy another constraint is added to the system which imposes ω ± (τ k+1 ) = ω ± (τ k ), ∀k, which means that the dividing cells do not change their orientation during the repositioning. By contrast, in the free orientation case, ω ± (τ k+1 ) = ω ± (τ k ), i.e., no constraint is imposed on the new orientation at step τ k . This new constraints are defined in 2.3.

Numerical Solution of the Model
We now detail the numerical method used to solve our model. The general structure of the algorithm is the following The computation of the statistical quantifiers is detailed in the next section.

Positioning
Step We discuss now step b)-(iv). In order to find a solution to the mini- where W Q is the global adhesion potential relative to the quadratic choice of the potential function V Q , we construct a method based on the Uzawa algorithm (Arrow et al. 1958). Given N (t) cells, the number of constraint functions φ i j (ξ(t), ρ(t)) due to the non-overlapping condition is M = N (t)(N (t) − 1)/2. Then, the algorithm consists in finding a saddle point of the Lagrangian function L Q (ξ(t), λ(t)) : where the λ i j are called the Lagrange multipliers. The algorithm constructs a sequence of approximate values (ξ(t) ( p) , where β and μ are numerical parameters and where the dependence on t has been omitted for simplicity and will also be omitted in the sequel of this paragraph if not strictly necessary for comprehension. After some computations, the first equation of the above system can be rewritten for k ∈ [1, N ] as which clarifies the role of the numerical parameter β in the scheme; it is related to the displacement of the cells during the search of an equilibrium position. Two stopping criteria, which need to be satisfied at the same time, are used in order to advance to the next step. They are based on measuring the following quantities where tol φ and tol W are two tolerances the values of which are given below. These criteria permit to control the largest overlapping permitted between the cells and to exit the algorithm when two consecutive values of the total mechanical energy of the system are very close to each other, indicating that a saddle point is likely to have been reached. Finally, the parameter μ is related to the speed at which the constraints are updated.
In order to reach a solution to the minimization problem as fast as possible, an adaptive β has been chosen which depends on the number of cells considered. In practice, β = 3 10 −4 for 1 ≤ N ≤ 100, β = 3 10 −5 for 100 ≤ N ≤ 300 and β = 6 10 −6 for 300 ≤ N ≤ 500, while μ is kept fixed to μ = 100. This reflects the observation that the Lagrange multipliers values grow with the number of cells N . Consequently, the value of β should diminish when N grows in order to avoid too large displacements of the cells which may lead to saddle points very far from the initial configuration and thus unrealistic. However, it may happen that when constraints are strongly violated, these choices for β are not sufficient to prevent ejection of cells from the aggregate. This is measured by computing the distance traveled by a cell between two consecutive steps ( p) and ( p + 1) of the minimization algorithm. If this distance goes beyond a fixed tolerance tol X , this is repaired by repeating the positioning algorithm with a choice of β which avoids too large displacements. Details on the value used for tol X are given below.

Growth Step
Step (i) consists of the simple implementation of the growth law discussed in the previous section. Given the parameters R min , R max , γ , T G and the time step t, we just sample a random number u between [−α, α] and we compute . After the growth, in general, an overlapping between cells is produced which is resolved by the repositioning step described in Sect. 3.1.

Division
Step We assume that the cell C i 0 is ready to start the division, i.e., R i 0 (t) ≥ R max . For each simulation, we fix the number of steps of the division process k = [1, f ], the initial direction ω i 0 (t) of the division plane and the division strategy. As soon as the cell begins its division, the cell C i 0 is replaced by two new cells. The algorithm can be summarized as follows.
where the variable t is fixed during all the division process. This series of actions causes the cells to partially overlap with their neighbors. This is corrected by a new application of the positioning algorithm.
The additional constraints imposed by the mitosis are, for both free and fixed orientation strategies, the change in the non-overlapping constraint between the cells is changed into an equality constraint between the two daughter cells (indexed by i + , i − ) which, for the iteration ( p) of the minimization algorithm relative to the generic division step τ k , reads as follows: while the corresponding Lagrange multiplier is updated accordingly to λ ( p) i + i − and where the dependences on τ k and t have been omitted for simplicity. In the constrained strategy case, two additional constraints should be added to the positioning algorithm to take into account that ω ± (τ k ) remains constant equal to while the new positions of the two daughters cells take into account these constraints through the corresponding Lagrange multipliers λ ( p) 1 and λ ( p) 2 as follows: are the positions computed during the iteration ( p + 1) of the positioning algorithm without divisions and where once again dependence on τ k and t are omitted.

Initialization and Numerical Parameters
The initialization is done by inserting N 0 cells in the computational domain with radius R min , at a random location, each cell defining a different lineage i , i = [1, N 0 ]. Then, a positioning step finds a first saddle point of the Lagrangian function L Q (ξ(t = 0), λ(t = 0)) which furnishes the initial positions of the cells in the tumor.
We list now all numerical values given to the parameters. We distinguish the model parameters listed in Table 1 from the numerical parameters listed in Table 2. In particular, considering Table 1, the choice of N 0 represents the effective number of initial cells used in average in the experiments. In the same spirit, T G and T max are also chosen to be as close as possible to experimental values. The HCT116 line used in the experiments has a cycle of an average duration of 24 h and they are typically tracked The ratio between the minimal and the maximal radii assumed by the cells corresponds to the average ratio of the HCT116 line. Finally, the value f (number of steps during the division process) is chosen to avoid too large cell overlapping during the division process. Concerning the numerical parameters, the time step is chosen to guarantee small enough cell size increments which avoids too large cell overlapping before repositioning. The values assigned to β have been already discussed in Sect. 3.2 together with μ. The value chosen for tol W is directly related to the values of β. The tolerance tol φ is chosen to permit only very small overlaps (of the order of 5 × 10 −2 for values of the radius R between R min and R max ). The value of tol X detects too large displacements of cells (of the order of twice the radius R min ). For the choice of the numerical parameters, we refer to Peurichard et al. (2017) where an extensive study is reported. Since minima of the mechanical energy subjected to non-overlapping constraints are not unique, the numerical parameters are chosen such that the algorithm selects the closest configuration to that of the previous time step.

Statistical Indicators
To get insight into both the experimental and the numerical results, we develop several indicators to measure the characteristics and the morphologies of the single lineages or of the whole population. The Voronoï Diagram In Fig. 3 is reported a typical result of a simulation. The left picture shows the initialization after the first positioning phase, with N 0 = 50 cells, while the right picture shows the solution at T max which corresponds to a situation with about N = 400 cells. This representation of the solution has several limitations due to the difficulty in defining the notions of cell neighborhood, perimeter and area. In order to overcome this problem, we use a modified Voronoï diagram representation. This representation is introduced for statistical analysis purposes only. It does not play any role in the cell dynamics. This approach is frequently used in the context of growing cell populations, see for instance (Osborne et al. 2010;Fletcher et al. 2014;Bi et al. 2016) where this method is used to determine the interaction forces between two cells.
We recall some basics about the Voronoï diagram in the present setting. A Voronoï site is defined as a point p i belonging to a predefined subset Thus, given a set of sites S, the Voronoï diagram partitions the plane by which site is the closest. Two sites p i and p j are considered as neighbors if R i and R j share a common edge. The intersection of three regions, if not empty, is called a Voronoï vertex or node.
In order to use this representation in our model, we define the Voronoï sites as being the cell centers, i.e., S = {X i |1 ≤ i ≤ N }. But this leads to a problem at the boundary of the cellular aggregate, since the tumor does not occupy the entire plane. The Voronoï regions corresponding to these outer sites are consequently unbounded. To overcome this problem, we add fictitious sites on the borders of the population. These additional sites are located on the polygon whose boundaries are the edges of the convex hull of the tumor slightly enlarged by dilation. More precisely, given a small δ > 0 and P the polygon obtained by enlarging the boundaries of the convex hull by an increment δ, we place n additional sites equally spaced along the segment of P, where n corresponds to the number of outer sites. The result of this procedure leads to Fig. 3 on the middle right and right where the Voronoï diagram has been traced for the same situation as on the left and middle left. In this figure, the Voronoï regions related to the fictitious sites are not represented. Now, the definitions of area, perimeter and neighborhood of cells or groups of cells are easier even if it should be stated that these quantities in the passage from the sphere representation to the Voronoï one do not remain constant. We will use the same Voronoï approach for the experimental data to have similar definitions of neighborhoods, areas and perimeters in the experimental and numerical results.
Diagnostic Definitions One of the main questions we address is the influence of the orientation of the division plane during mitosis on the morphology of the lineages. We consider all descendants of a given ancestor cell present at the initial time and store this information in the lineage i , i = [1, N 0 ]. Thus, if cell C i is present at the initial time, then i = i; otherwise, when a cell is created its i becomes the lineage of its mother cell. To define diagnostics, we use several concepts of graph theory. We briefly recall them. A graph G is an ordered pair comprising a set of points or vertices and a set of edges or links where an edge is related with two vertices and two vertices of the graph may or may not be connected by a link. A chain in G is a finite sequence of vertices connected by links with possible repetitions. A cycle is a closed chain. We are now ready to introduce some definitions: Neighboring Cells Two cells C i and C j are said to be neighbors if their Voronoï regions R i and R j share a common edge.

Graph of N Cells
The set of all cells is a graph where the vertices are the cells centers.
Connected Set of Cells The cells C i 1 , C i 2 , . . . , C i p are said to be connected when the subgraph induced by these cells is a connected graph.

Connected Component of a Lineage A connected component of a lineage is a connected component of the subgraph generated by the cells of the same lineage.
Cell Cycle A cell cycle is a finite sequence of cells corresponding to a cycle in the graph. In other words, for all k ∈ [1, p + 1], C i k and C i k +1 are neighbors and C i p+1 = C i 1 .
Cell Polygon Associated with a Cell Cycle Let p ≥ 3 and a cell cycle C i 1 , C i 1 , . . . , C i p+1 . The line passing from the centers of the cells X i 1 , X i 2 , . . . , X i p+1 forms a polygon which is called the cell polygon associated to C i 1 , C i 1 , . . . , C i p+1 .
Boundaries of a Set of Connected Cells Given a set of connected cells, the boundary is defined as a cell cycle whose associated cell polygon contains all cells of the connected set. This polygon is called the boundary polygon.
Perimeter and Area of a Set of Connected Cells Let us consider a set of connected cells with p ≥ 1 and its boundary denoted by C i 1 , C i 1 , . . . , C i p+1 . The perimeter of this array of cells is defined as the perimeter of the cell polygon associated to C i 1 , C i 1 , . . . , C i p+1 . It is given by P = p k=1 d(C i k , C i k+1 ). The area of this set is defined as the area of the cell polygon. It is given by

Convex Hull of a Set of Cells
The convex hull a set of cells is defined as the convex hull of the point cloud consisting of all the cell centers.
We are now ready to define the statistical indicators used for studying the morphology of the cellular aggregate.
(1) Sphericity of the Population ratio between the area and the perimeter squared of the entire cell population: R 1 = 4π A P 2 where R 1 ∈ [0, 1]. R is one when the boundary of the population is a perfect circle. (2) Convexity of the Population ratio between the area A of all cells and A conv , the area of their convex hull: Its value is one when the boundary of the population coincides with the boundary of its convex hull.
(3) Sphericity of a Lineage for each connected component (of at least two cells) of a lineage, ratio between the area and the perimeter squared of this connected component: R 3 = 4π A P 2 where now P and A are, respectively, the perimeter and area of the boundary polygon of the considered connected component.
(4) Lineage Fragmentation number of connected components in which a lineage is split: R 4 . It permits to understand whether the cells of the same lineage tend to remain grouped or to scatter. (6) Lineage Orientation Orientation direction for a given lineage R 6 . The main direction of orientation is computed by using the inertia matrix of the cells composing the lineage.
Inertia Matrix Let p ≥ 2 and for k ∈ [1, p], X i k = (x i k , y i k ). Denoting X G = (x G , y G ) the barycenter of X i k : X G = p k=1 X i k p , the inertia matrix of the cloud is defined by Assuming that the X i k are not aligned, this matrix is symmetric and positive definite with strictly positive eigenvalues. Thanks to these values we can measure the angle θ 2 formed by the semi-major axis of the inertia ellipse and a reference direction, and the angle θ 0 formed by the line joining the origin to the barycenter of the lineage this a reference direction. The angle θ 2 is the angle which measures the direction of the eigenvector v 2 associated to the larger eigenvalue λ 2 . We then measure the quantity R 6 = θ 2 − θ 0 − π 2 where the shifting of π 2 is done in order to have an angle always between − π 2 and π 2 . In Fig. 4 is reported an example where the main direction of the ellipse (continuous line) representing the inertial matrix together with the ellipse (dotted line) is shown for three different lineages.

Experimental Protocol
In order to generate reference experimental data to which mathematical modeling can be confronted, we set up the following protocol. The experiments are performed on cells of the colon adenocarcinoma HCT116 cell line, modified by lentiviral transduction to express a histone H 2 B fused with the mCherry fluorescent protein. This allows visualizing by fluorescence microscopy the nuclei of the cells. The cells are seeded in culture chambers (Lab-Tek, Dutscher) at a density of 7500 cells/cm 2 in an OPTIMEM medium supplemented with 3% of fetal calf serum (FCS) and penicillin/streptomycin. After 48 h, the chosen cell density provides in the bottom of the chamber a culture of islets composed of about 50 cells, isolated from one another. This allows following the individual cell evolution in real time by an inverted fluorescence microscope. Before making microscopy acquisitions, isolated groups of cells are selected and the bottom of the culture chamber manually processed, to prevent neighboring cells to join the main group. Acquisitions are performed on an Axiovert microscope (Zeiss) fitted with a CoolSnap HQ camera (Roper Scientific) and piloted by the MetaView software. For every acquisition, several individual groups are followed in parallel by videomicroscopy (1 image every 10 min). The images are processed using Metamorph and Image J before being analyzed.
The data post-processing consists of the following steps repeated for each single cell culture. Cells are selected at the initial time and the filiation tracked manually through direct labeling. Figure 1a, b shows an example of this procedure. Segmentation is first performed automatically thanks to the level of fluorescence, then checked and, if necessary, results are corrected by post-processing the data manually. This is possible thanks to the fact that the analysis is made on a relatively small number of cells. A segmentation result is shown in Fig. 5. A post-treatment is performed after the segmentation procedure. It consists of identifying each single lineage. The result is shown in Fig. 6 where the same experiment as in Figs. 1a, b and 5 is considered.
We describe now how the experimental data are processed. The areas detected during segmentation correspond to the nuclei of the cells. For each of these areas, The colored disks correspond to the cells whose lineage is tracked. The same experiment as for Figs. 1a, b and 5 is reported the coordinates of the center of mass and the equivalent circular diameter, i.e., the diameter of the disk that has the same surface as the area of interest, are computed. This permits to have a first representation of the data by disks (Fig. 7), similar to that used in the simulations. However, this gives large gaps between some disks and overlapping between others and it does not correspond to the observed experiments. Indeed, this is only an artifact of the representation because cells are stuck together and do not overlap. We then choose to represent the population by the same adapted Voronoï diagram as that used for the simulations. This result is reported in Fig. 8.
The lineages used for the comparisons with the simulations are either at the periphery of the aggregate or in the central region. A lineage is considered as peripheral if, at initial time, the progenitor cell is situated at the periphery of the aggregate, while it is The colored Voronoï regions correspond to the cells whose lineage is tracked. The same experiment as in Fig. 1a, b and 5 is reported considered as central if, at initial time, there are at least two cells between the progenitor cell and the boundary of the tumor. The data at our disposal are based on fourteen cellular cultures. This corresponds directly to the size of the samples which has been used in the first and the second diagnostic. In these fourteen experiments, we have followed twenty-two central and fifty-four peripheral lineages which corresponds to the size of samples used in diagnostics 4 and 6. The twenty-two central lineages split up into fifty-three connected components, while fifty-four peripheral lineages divided into one hundred and fifteen connected components which correspond to the size of samples for diagnostics 3 and 5.

Identification of the Different Morphologies
Since we have chosen three possible orientation planes and two division strategies, for each of the six possible situations we perform 100 numerical simulations and we compute the values of the statistical indicators R 1 up to R 6 . Concerning the study of the lineages, we extract two lineages, one in a central position and one on the boundaries of the cell population for which the different diagnostics are computed. This choice of the different position of the lineages is done in order to highlight the strong disparity between the behaviors of internal and external cells. This disparity which is one of the major results of this work is also confirmed by in vitro experiments as detailed next and not known before. The difference between the central and peripheral lineages is an information which is extrapolated from diagnostics R 3 up to R 6 . In the appendix, we report the detailed results of the simulations. In this part, we summarize the principal results.
(1) Presence of larger numbers of lineage fragments at the periphery, typically of smaller sizes. This can be explained by the fact that, as cells are farther from the center, their contribution to the total energy of the system grows. Thus, when they are pushed from the interior of the cellular aggregate because of the growth, they preferably move in tangential direction not to excessively increase the energy of the system. This result in an intercalation of new cells coming from the interior between originally neighboring cells. (2) Preferred tangential orientation of the lineages at the periphery. An interpretation of this phenomenon is that far from the center, the potential energy obliges cells to place themselves on an equipotential energy curve, i.e., along a single shell, all other equilibrium states being unstable. (3) In the central zone, except for the constrained radial division strategy, there are no privileged directions for the lineages. This suggests that in the central area, the division rule, when radial, plays a more important role than the positioning rule, this feature being opposite for the other division rules. (4) In general all strategies which leave freedom to the division plane orientation to change during division have a very small impact on the overall population shape as well as on the shape of a single lineage. On the contrary, the constrained strategy has much more influence on the shape of the entire population and of single lineages. (5) Concerning the global population, only the constrained strategy permits to detect differences between the radial and tangential directions of the division plane. The boundaries of the growing cell population are smoother for this second choice. (6) With the constrained radial orientation strategy for the division plane, the lineages in the center of the population are radially oriented with few fragments of large sizes.

Results of the Experiments
The same diagnostic analysis done for the mathematical model is performed on the experimental data. The results can be summarized as follows.
(1) Circularity of the lineages is very low compared to that of the entire population (the average value of R 1 ∈ [0, 1] is R 1 = 0.56. (2) Lineages are filamentous (the average value of R 3 ∈ [0, 1] is very small, around 0.12). Central lineages are more filamentous than peripheral lineages: there is a larger number of lineages with R 3 = 0 in the central part compared to the periphery. (3) There is a strong disparity (measured by R 4 ) between the central and the peripheral zones which is expressed by a stronger fragmentation of the central lineages compared to the peripheral ones. In the central zone, the number of lineages divided in one, two or three pieces is almost the same while on the periphery the number of lineages divided in two and in three parts is lower than the number of lineages which did not divide. (4) The central lineage fragments, measured by R 5 , are much smaller than those of the periphery. 2/3 of the fragments are composed of one or two cells in the center while on the periphery they are around 40%. (5) Half of the central lineages have a direction very close to the radial one, while in the periphery no preferred orientation seems to arise (this value is measured by R 6 ).

Comparisons with the Experimental Data
The comparisons between the simulation results and the experimental outcomes lead us to the following conclusions shown in Fig. 9. In this figure, the results for all the (1) Regarding the first two diagnostics, i.e., R 1 and R 2 , which are not reported in the figure and which study the entire cell population, the experimental values are on average slightly lower than those obtained by the numerical simulations.
In experiments cells organization is less regular. In particular boundaries are less smooth and the shapes are less round. This is likely due to the fact that real cells are less regular than the perfect disks chosen in the mathematical model. However, the results for these two diagnostics are qualitatively comparable and the choice of perfect disks does not seem to affect other results.
(2) Both simulations and experiments highlight a difference between center and boundaries of the growing cell population. However, this disparity is expressed differently (indicator R 6 ): in the experiments, no preferred tangential orientation of the lineages appears at the periphery of the population (or more precisely only a slight preference for the tangential direction), while a clear preferred tangential orientation is obtained in the numerical simulations for all the different orientation strategies considered. On the other hand, a preferred radial orientation is visible for the central lineages in the experiments, which is very close to the results obtained when the constrained radial orientation strategy is used for the numerical simulations.
(3) In the mathematical model, the peripheral lineages are more filamentous more fragmented into smaller fragments than the corresponding central lineages, while in the experiments the situation is reversed: peripheral lineages are less fragmented with larger size fragments than those in the center (from diagnostic R 4 and R 5 ). Globally we can conclude that the initial choices done for constructing the mathematical model permit to reproduce some of the observed features, but they are not sufficient to correctly describe all the cells behaviors measured in the experiments. In particular, the difference observed at center of the aggregate between the data and simulations especially in diagnostic R 4 (lineages fragmentation) suggests that some additional mechanism is responsible for the displacement of the central cells to the periphery. This mechanism disrupts the organization of the cells and makes central lineages more fragmented and filamentous. In the next section, we propose such mechanism and show that it results in simulations being closer to experimental data.

Improved Model: Bounded Confinement Force and Cell-Cell Interchange
Here, we propose additional mechanisms to reconcile the simulation results with the experimental data. The section is divided into three parts; in the first part, we discuss two improvements to the model. In the second part, we detail the modifications of the algorithm necessary to take into account these modifications. In the third part, we discuss comparisons between the new results and the experimental data.

Model Improvements
The first modification consists of the possibility for two adjacent cells to switch their positions. The second modification consists of modifying how the interaction potential depends on the distance from the center of the tumor. The first mechanism permits to switch the position between a new born cell and a neighboring cell. Indeed we hypothesize that compression by the other cell may induce deformations of the cell membrane and that the so deformed cell may be able to migrate into the extra-cellular medium. The random switch between two cells is a way to model this migration favored by cell deformation. We only allow newborn cells to shift their position with a neighboring cell because newborn cells have the smallest size and are more likely to find a migration path in the extra-cellular medium than mature cells. Finally, we only allow cells from the central region to perform this switch because cells from outer regions are subject to lower compression and weaker deformations which decreases their ability to migrate. We consider the possibility for a cell C i to switch if its position is inside a disk defined by with C int ∈ [0, 1] a modeling parameter. The permutation rule is then as follows: after a division of a progenitor which lies inside this disk, one of the two daughter cells has the possibility to switch its position with one of its neighbors. The decision to switch is modeled by a probability following a Bernoulli distribution of parameter p: with probability p a newborn cell inside the sphere switches its position with a neighboring cell.
The second additional mechanism consists of applying a weaker attractive potential to the cells which are far from the center of the aggregate. With this aim, the quadratic potential is replaced by a linear potential when the tumor reaches a critical size and the linear potential applies only to cells which are farther than a critical distance. This reflects the hypothesis that cells are submitted to a lower mechanical compression at the periphery of the aggregate. However, for small-sized aggregates, we hypothesize that the adhesion forces between the cells are stronger (since a strong grouping enhances the survival chances of the cells) which motivates the use of a quadratic potential. The positioning rule is consequently modified as follows: as soon as the number of cells reaches a critical value N C , a linear potential is used for the remote area while the same quadratic potential as defined in Sect. 2 is used for the other cells. The modified global adhesion potential is where R c = C L max 0≤ j≤N (t) (d(0, X j (t))) while the remote area is defined as the set containing the remote cells and consequently a cell at position X (t) is said remote if with C L ∈ [0, 1] a modeling parameter.

Algorithm Adaptation
The new structure of the numerical algorithm is the following We discuss the modified positioning step and the permutation step. The modified positioning step consists in finding a saddle point of the new Lagrangian function Thus, starting from an initial guess (ξ(t) (0) , λ(t) (0) ), the method reads, as in the previous case, as where β and μ are the same numerical parameters as those discussed in Sect. 3. After some computations, the first equation of the above system can be rewritten for all cells in the remote region, i.e., |X j | > R c , as The same stopping criteria are used for this new positioning algorithm. During the research for a saddle point, it may happen that a cell close to the boundary between the central and the remote regions changes zone. Since this displacement is typically very small, we choose not to change the potential energy to which this cell is submitted during the minimization procedure. The permutation algorithm consists simply in choosing with probability p if a newborn cell performs a switch. In order to do that, denoting by C I the cell which has decided to perform a switch, we first determine all the neighboring cells of C I and then we pick randomly one of them with uniform distribution and we perform the switch. The values of the new parameters added to the model are summarized in Table 3. The value of C L is chosen so that the central region coincides with the definition of a central lineage in the experiments.

Comparisons Between the Experiments and the Improved Numerical Model
We first analyze separately the impacts of the two modifications introduced in the model. We start discussing the effect of the switch. The results for diagnostic R 4 (lineage fragmentation) and R 5 (size of the fragments of a lineage) are reported, respectively, in Fig. 10 and 11 for three different values of the permutation probability p: p = 0.1, p = 0.15, p = 0.2. They show that the larger the switching probability, the more fragmented the lineages are and the smaller the fragments are. Cell switching even though it is performed only in the central region has also an influence on the peripheral region. However, the impact of these changes is low and we do not report it. The direction of lineages computed with diagnostic R 6 does not change noticeably with the switching. Now, we consider the influence of the change in the attractive potential. The results for diagnostics R 4 and R 5 for two different values of the parameter N s , i.e., mechanism at the periphery of the tumor. In particular, we observe that the smaller N s is, the less segmented the lineages are and the bigger the fragments are at the periphery. The change in the potential energy has also an effect in the central zone, where the lineages are becoming less fragmented. Finally, regarding the direction of the lineages with diagnostic R 6 , we do not observe noticeable compared to the quadratic potential. We finally observe the results obtained while combining the two modifications. As shown in Fig. 14, the simulated results obtained with the two model improvements are more similar to experimental results, both in the central zone and at the periphery. In the figure, the results of the Diagnostics 4 and 5 are reported. For these tests a switching probability for a daughter cell of p = 0.2 and a threshold of N = 200 to pass from a quadratic to a linear potential have been chosen. These parameters permit the best fit between the simulations and the data. We conclude that these modifications improve the match between the model and the experiments. The distributions of fragments number and size of lineages obtained with the modified model match those of the experimental data remarkably well.

Conclusions
This work represents a step toward understanding the impact of division parameters on the growth of a cell population via comparison of mathematical models and experiments. The approach consists of a bottom-up strategy where the behavior of a growing population of cells and the structure of the associated lineages are modeled through simple interaction rules of mechanical type between cells. After observing that these simple rules could not explain the morphology of the cell population alone, we introduced additional phenomena, which, at first were considered negligible. This process of gradually increasing complexity was repeated until a good fit between experiments and simulations is obtained. This approach permits to assess which key mechanisms are more likely to be underlying the observed phenomena without introducing too many empirical parameters. From the mathematical point of view, the energy minimization considered here is motivated by the observation that physical principles are often expressed in variational forms. Our results suggest that this variational approach seems also at play here. Our numerical and experimental results show a wide disparity between central and peripheral lineages. Peripheral lineages are more fragmented and slightly tangentially oriented. With the simplest model, some differences between simulations and experimental data are found: experiments show that cells in the central region move over larger distances than in the model. We thus introduced the possibility for cells to switch positions inside the cell population. We also lowered the aggregation force at the periphery of the cellular aggregate to model weaker aggregation of peripheral cells. The corresponding numerical results match the experimental data very convincingly. Another interesting experimental observation is that in the central area, a preferred radial direction of the lineages emerges and that a similar feature can be found in the model if a radial division plane orientation is imposed. The analysis suggests that the disparity between the center and the periphery, found experimentally and verified numerically, could be explained by the simple hypotheses made in the model. More quantitative work is needed to understand the role of the division orientation on the observed emerging structures and on lineage shape and orientation. In the future, larger populations should be considered. On the other hand, considering very small populations permits to ignore the role played by nutrients, growth factors and oxygen, and to consider simple two-dimensional settings for both the model and the experiments. This also made the image processing easier and more reliable avoiding problems related to the tracking of cells in a three-dimensional structure and it allowed for the use of CPU-effective agent-based model which would be too costly in three dimensions. However, these simplifying choices lead to a model with restricted validity and questionable applicability to three-dimensional structures. In the experimental setting considered, cells move on a plane. Thus, our two-dimensional model is a valid representation of the biological situation. However, we may expect the results to be different in a fully three-dimensional experiment and not well accounted for by a two-dimensional model. In particular, we may expect that the differences between the internal and the peripheral lineages that we observe in the two-dimensional case will be even stronger in the three-dimensional situation as the scale ratio between the surface of the periphery of the spheroid to its internal part will change. To understand the role of the orientation of the division plane, future experiments which track cells during mitosis would permit to further explore the biology of tumor aggregates and provide a better benchmark for the validation of three-dimensional models. Another direction of research is to explore the impact of mechanical confinement on cell proliferation. Recent experiments (Desmaison et al. 2013b) showed that proliferation gradients within mechanically confined spheroids are different from those in spheroids grown in suspension. This discovery strengthens the hypothesis of mechanical forces playing a central role in the morphologies of the lineages and requires further studies. Diagnostic 2. Convexity of the Population R 2 : results are consistent with the ones of the previous diagnostic. Whatever the strategy and division direction chosen, the values of R 2 are very similar, around a mean value of 0.94, as reported in Fig. 16. Giving a closer look at the results, we see that the differences between the three division orientations are more pronounced when the constrained strategy is employed. In particular we can state that if the free strategy is adopted, no difference between the three orientations of the division plane is observed. In the constrained case, results show larger R 2 values for the tangential direction, which can be interpreted by saying that the shape of the boundaries is more regular when the division plane has tangential direction.
Diagnostic 3. Sphericity of the Lineage R 3 : in Fig. 9a, c we plot the distribution of connected components of a given lineage as a function of R 3 . On the left picture, the results for a lineage situated close to the center of the tumor are shown while on the right picture those for a lineage at the periphery of the tumor are displayed. All the different adopted strategies are averaged together. The results clearly evidence differences between the center and the peripheral zones on the morphology of the lineages. The sphericity of the lineages is much lower at the periphery than in the  , d, g, j). Middle radial division orientation (b, e, h, k). Right tangential division orientation (c, f, i, l). a-c Free orientation, d-f constrained orientation. g-i Free orientation strategy, j-l constrained orientation. a-f Lineages which lie at the center of the tumor, g-l lineages which lie at the boundary of the tumor central region. In particular, almost half of connected components have R 3 = 0 at the boundary of the aggregate while only one third of the components have R 3 = 0 in the central region. We can conclude that the lineages are more likely to remain grouped in the central region. This is true for all the tested division orientations and strategies. For each of these, the results are reported separately in Fig. 17. Diagnostic 4. Lineage Fragmentation R 4 : the number of connected component for central and peripheral lineages R 4 is reported in Fig. 9d, e for all the different strategies averaged together. The main result we can draw is that a strong disparity occurs between the central and peripheral lineages. In particular, almost all lineages  , d, g, j). Middle radial division orientation (b, e, h, k). Right tangential division orientation (c, f, i, l). a-c free orientation, d-f constrained orientation, g-i free orientation strategy, j-l constrained orientation. a-f Lineages which lie at the center of the tumor, g-l lineages which lie at the boundary of the tumor are divided from one up to four connected components in the central region and from one to five connected components when at the periphery. In the central region around 57% of lineages remain connected while 33% of the lineages are separated into two connected components. On the periphery of the tumor, the proportions are 33 and 34%, respectively. We can conclude that lineages are more fragmented at the periphery. This is probably due to the fact that divisions inside the aggregate push cells farther from the central region which results in their intercalation between cells of the lineages of the periphery. The complete set of data for all the studied different strategies is  , d, g, j). Middle radial division orientation (b, e, h, k), right tangential division orientation (c, f, i, l). a-c Free orientation, d-f constrained orientation, g-i free orientation strategy, j-l constrained orientation. a-f Lineages which lie at the center of the tumor, g-l lineages which lie at the boundary of the tumor reported in Fig. 18. We can notice that two cases give clearly remarkable results: in the central region, the radially constrained division orientation presents lineages that are divided almost always in one or two components. On the other hand, at the periphery, the constrained tangential orientation of the division plane presents many lineages divided in four or five components.  Fig. 14g, l. The figure shows the cumulative distribution of connected components as a function of the number R 5 of cells of the connected component. The figure reports also the profile of this distribution between five zones: from zero to one cell, from one to two, from two to seven, from seven to eight cells and from eight to larger connected components. This profile is obtained by computing the linear interpolation lines of the cumulative distribution between the ends of these five zones. On the periphery of the cellular aggregate, there are 50% of small connected components and 16% of large connected components. By contrast, in the central region, there is a higher percentage (almost 30%) of small connected components and a similar percentage (34%) of large connected components. In the central and peripheral regions, the percentage of lineages of medium size (three to seven cells) is almost identical (around 7.5%). There is a sort of phase transition from small to large connected components. It can be concluded that, at the periphery, the connected components are more numerous and composed of a fewer number of cells. In other words, more cells are isolated from their lineage. Full data, strategy per strategy, are shown in Fig. 19. From this figure, we can observe that the free division strategy for the three considered division directions gives very similar results, whereas the fixed strategy enhances the differences between the different possible orientations of the division. Two situations are remarkable. The first one is the radially constrained strategy which causes the formation of many large connected components in the central region. The second remarkable situation is obtained for the constrained tangential orientation strategy. For this situation, on the boundary of the aggregate we obtain many isolated cells separated from the rest of their lineage.
Diagnostic 6. Lineage Orientation R 6 : in Fig. 20 we report the cumulative distribution of lineages as a function of |R 6 |. A concentration of the frequencies near 0 indicates that the main direction of the lineages is tangential while a concentration close to π/2 indicates a radially oriented lineage. At the periphery, there is a high percentage (almost 50%) of lineages the orientation of which is in the tangential direc-  , d, g, j), middle radial division orientation (b, e, h, k), right tangential division orientation (c, f, i, l). a-c free orientation, d-f constrained orientation, g-i free orientation strategy, j-l constrained orientation, a-f lineages which lie at the center of the tumor, g-l lineages which lie at the boundary of the tumor tion. In the central region, no direction seems privileged. The tangential direction is strongly favored by the lineages at the periphery probably because of the influence of the positioning rule. The complete data are shown in Fig. 21. From this figure, we can observe that for the free strategy the three orientations of the division plane give very similar results whereas for the constrained strategy differences become more visible. Moreover, for the radially constrained division strategy the results depart significantly from those of the other strategies. In the central region, this orientation gives 50% of lineages whose direction angle is less than 25 • . At the periphery, for tangential or random division orientations, the three quarters of the lineages exhibit a difference with respect to the tangential direction of less than 35 • , while a similar result is only reached for the radial division strategy when considering all the lineages until the value of 60 • .