Self-organization in brain tumors: How cell morphology and cell density influence glioma pattern formation

Modeling cancer cells is essential to better understand the dynamic nature of brain tumors and glioma cells, including their invasion of normal brain. Our goal is to study how the morphology of the glioma cell influences the formation of patterns of collective behavior such as flocks (cells moving in the same direction) or streams (cells moving in opposite direction) referred to as oncostream. We have observed experimentally that the presence of oncostreams correlates with tumor progression. We propose an original agent-based model that considers each cell as an ellipsoid. We show that stretching cells from round to ellipsoid increases stream formation. A systematic numerical investigation of the model was implemented in R2. We deduce a phase diagram identifying key regimes for the dynamics (e.g. formation of flocks, streams, scattering). Moreover, we study the effect of cellular density and show that, in contrast to classical models of flocking, increasing cellular density reduces the formation of flocks. We observe similar patterns in R3 with the noticeable difference that stream formation is more ubiquitous compared to flock formation.

Introduction Primary brain tumors are one of the most lethal cancers. In spite of surgery, radiotherapy and chemotherapy, median survival remains at 14-18 months. The key to develop successful cancer therapy is to understand the essential mechanisms by which individual cancer cells proliferate, grow as a tumor, and invade normal brain. It is of particular importance to understand how individual cell morphology relates to collective macroscopic behaviors (e.g. stream formation, diffusion behavior). As our data indicate that glioma oncostreams promote tumor growth, this raises the question of whether cell morphology influences pattern formation and therefore the overall dynamics of growing tumors. This question is difficult to answer as it requires access to the time evolution of the positions of the cells in vivo.
We wished to explore the relationship of cell morphology to collective microscopic behavior patterns using mathematical modeling as it has already been successful in the exploration of various biologically relevant scenarios [1][2][3][4]. In particular, agent-based models (modeling at the microscopic level) are convenient as they incorporate essential features of cell behavior (i.e. motility, cell-cell interactions, etc.) and have been exploited to understand various self-organizing dynamical systems (e.g. pedestrians [5], birds [6], fish [7], bacteria [8,9]).
To investigate the influence of the shape of the cell on tumor dynamics, we modeled cells as ellipses or ellipsoids. This assumption is motivated by experimental observations (see Fig 1) where cells within oncostreams display a length to width ratio of 2.7:1. Using ellipsoid shape is common in the study of bacteria, for instance viscoelastic ellipsoids have been used in [10] or self-propelled spheres in [11] (see also [12][13][14][15][16][17]). We were particularly interested in studying the dynamics in a regime of high cellular density where cells are always in contact with each other. 'Stretching' the cells' in this regime could potentially increase the formation of streams since streams would reduce overlapping of elongated cells. Indeed, in the context of soft-mater with elongated cylinders (e.g. nail, log, rice), stream formations are ubiquitous [18][19][20].
We propose an agent-based model that utilizes two mechanisms: i) self-propulsion, ii) cellcell avoidance due to non-overlapping constraints. Since the cells have an ellipsoid shape, cellcell avoidance leads to two possible effects: repulsion (i.e. cells move away from each other) and steering (cells turn to avoid collision). The larger the eccentricity of the cell, the larger the effect on steering. In contrast to classical models of flocking [21,22], in our model cells do not take into account the velocity of their neighbors.
We first investigated numerically in R 2 how eccentricity influences flock formation (i.e. all the cells moving in the same direction) using as an indicator the polarization of the configuration. We observed that increasing eccentricity increases polarization. Surprisingly, this effect saturates and even becomes counterproductive as flock formation becomes less likely when eccentricity exceeds a threshold (eccentricity e �.7). Then, we studied how cellular density affects the dynamics by increasing the number of cells while maintaining the same size of the domain. Since we do not suppose a mean-field type interaction (there is no averaging in the interaction), increasing slightly the density could lead to drastic changes [23]. In our dynamics, we observed the emergence of streams when the density becomes large, meaning that cells are aligned but not necessarily moving in the same direction. We measure streams using the nematic average where we identify a vector ω and its opposite −ω.
Beside the influence of cell morphology, the other key component of the dynamics is the strength of both the repulsive effect and the steering effect, as each determine the two coefficients α and β, respectively. One could speculate that increasing the steering effect (i.e. larger β) would enhance alignment and therefore lead to flocks or streams. Our numerical investigation revealed this not to be the case. Particularly at large densities, it is only when β is small that a flock or a stream emerge. This result seems counter-intuitive. However, we need to emphasize that the alignment in our dynamics is only indirect, as cells do not perceive each other's velocity. Thus, it is an interplay between spatial constraint and steering that leads to the emergence of a stream or flock. Increasing a single parameter (repulsion or steering) does not necessarily enhance alignment.
Stream formation is more challenging to observe in R 3 since cells avoiding each other no longer move aligned or in opposite direction as in R 2 . However, our simulations show that flock and stream formation do still occur in R 3 providing that we maintain a large density of cells in the domain.
The complexity of the dynamics uncovered shows that it is difficult to predict a priory the effect of each mechanism. Therefore, it would be of great interest to develop a multi-scale approach to study the dynamics from a macroscopic viewpoint [24][25][26][27]. Moreover, this will facilitate data-model comparison [28,29], as much of the experimental observations are made at a macroscopic scale. Investigating the partial-differential equation associated with the dynamics [30][31][32] could provide a way to bridge this gap.
The manuscript is organized as follows: we first present the agent-based model in section 1, then we study how the cell morphology influences the dynamics in section 1. A systematic numerical investigation of the model in R 2 varying two key parameters is performed in section 1 which produces several phase diagrams of the dynamics at various densities. We explore the model in R 3 in section 1 and draw our conclusions and future work in section 1.

Material and methods
We propose an agent-based model to describe the motion of individual glioma cells. The dynamics combine cell-motility (i.e. self-propulsion) and cell-cell interaction (e.g repulsion or adhesion). Specifically, we consider N cells described with a position vector x i 2 R d with d the spatial dimension (d = 2 or 3), moving with velocity cω i where c > 0 is the speed (supposed constant) and o i 2 S dÀ 1 the velocity direction. The main novelty of the model is to consider an elliptic or ellipsoid shape for each cell. Thus, we consider two axes denoted a and b for (respectively) the major and minor axis (see . As two cells cannot occupy the same spatial position, cells will push each other if they are too close. Thus, we define an interaction potential V i between cells that measures the tension exerted on cell i generated by the surrounding cells: The quantity r ij is referred to as the normalized distance between the centers of the cells i and j. For instance, if a = b we recover that r ij is simply the norm kx j − x i k/a. The modification takes into account that the cell is more sensible to an obstacle in front rather an obstacle on the side. The model is defined in R 2 (i.e. d = 2) and can be generalized to R 3 by defining r ij as follows: where e 2 (0, 1) is the eccentricity of an ellipse defined as e ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi To prevent overlapping, the function F has to be singular at the origin. We choose the following smooth function with compact support (see As r ij decreases, F increases resulting into repulsion. We have now defined all the quantities required to define our agent-based model.  where α, β and c are positive constant, V i is the tension defined in (1) and P o ?
the projector operator onto the the normal plane to ω i (it ensures that ω i stays of norm 1).
In order to reduce the tension generated by neighboring cells, a cell can either move away (i.e. repulsion effect) or change its direction (i.e. steering effect). Both maneuvers are pondered by the coefficients α and β representing the strength of each effect. Using the expression of V i (1), one can deduce an explicit expression of the dynamics (see S1 Text). Notice that if the cell has a circular shape (i.e. a = b and the eccentricity e = 0), its orientation will remain constant i.e. _ o i ¼ 0. Indeed, in that case, turning will have no effect on the reduction of the tension V i . Thus, steering effects can only occur when if a 6 ¼ b.
Remark 2 Notice that r ij cannot be considered a distance between cells i and j as it is not symmetric (i.e. r ij 6 ¼ r ji in general). Thus, the influence of the cell i on j is in general different from the cell j on i, i.e Fðr 2 ij Þ 6 ¼ Fðr 2 ji Þ.

Eccentricity effect on the dynamics
Eccentricity induces alignment. Our first investigation of the agent-based model (4) and (5) is concerned with the impact of the morphology of the cell on the global behavior of the dynamics. As stated above, cells have perfect rounded shape when the two parameters a and b are equal (i.e. eccentricity e is zero) whereas they have elliptic or ellipsoid shape when a > b (i.e. e > 0). We varied the eccentricity e and measured how this change affects the cells spatial configuration.
Before varying the morphological parameters a and b, there are several other parameters to be determined in our dynamics. When possible, we use experimental values that have been quantified in vivo. For instance, it has been observed that glioma cell size varies in between 5μm to 20μm for their diameter and their speed varies around 10μm/h [33]. However, some parameters cannot be inferred from experimental observations such as the strength of the repulsion α and the steering β. A more detailed investigation of these two parameters will be conducted in the next section. For now, we fix their values as indicated in Table 1.   (Fig 3-right). The full simulation is also available (see S1 Video).
Statistical characterization. To further investigate the dynamics, we introduce several statistics to characterize the emergent behavior.
Definition 3 Consider a velocity distribution fo i g i¼1::N 2 S dÀ 1 . We denote by ψ the polarization: Similarly, we define the nematic polarization [12]: g ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where θ i is the angle between the direction ω i and the horizontal axis and hi denotes the averaging over the indices i (i.e. h cos . We define the configuration as a flocking configuration (i.e. cells moving in the same direction) when the polarization ψ � 1 and the nematic polarization γ � 1. We define the configuration as a streaming configuration (i.e. cells' directions are parallel but not necessarily moving in the same direction) when the nematic polarization γ � 1 and ψ < 1.

Remark 4
The nematic polarization can be generalized in higher dimensions (see S1 Text). PLOS COMPUTATIONAL BIOLOGY Self-organization in brain tumors: Influences of cell morphology and cell density The previous statistics only involve the velocity of the cells ω i . We propose a third statistics to also characterize the spatial configuration.
Definition 5 Consider a spatial configuration fx i g i¼1...N � R d and a fixed radius R. We say that cell i is linked to cell j if the distance between the two particles is less than R. This defines a relationship (i.e. i * j) with i 6 ¼ j defined: Clusters C k are defined as the connected components for this relationship: two cells i 0 and j 0 belong to the same cluster if there exists particles i 1 , . . ., i k (a path) such that The cluster size jC k j denotes the number of cells belonging in the cluster k. The average cluster size � jCj is defined as the expected cluster size jC k j over all the positions: where Cðx i Þ denotes the cluster containing the cell i. We illustrate the three statistics used in Fig 4. In Fig 5-left, we measure the value of the polarization ψ over time for different shapes of the cells by varying the coefficients a and b. When the cells have a circular shape (a = b = 4μm, e = 0), the polarization ψ remains close to zero for all time whereas it increases up to a maximum close to 1 when the eccentricity is greater than zero. The relation between eccentricity and polarization is however non-trivial: increasing the eccentricity does not necessarily lead to large polarization. For instance, the polarization with eccentricity e = .89 is significantly smaller than with eccentricity e = .84.
To further investigate the relationship between polarization and eccentricity, we plot in Fig  5-right the polarization at the final time for several experiments (changing the seed for the initial condition) and various eccentricities e. We then perform a local regression ('loess' method) to estimate the expected polarization ψ as a function of e. We observe that increasing the eccentricity e leads to larger polarization up to e �.7 but then the polarization quickly decays for larger eccentricities.
Indirect alignment. In classical models of flocking [21,22], each individual has access to the velocity of its neighbors. By relaxing its own velocity toward the average velocity of its neighbors, a flock emerges. This begs the question on how individual agents communicate. However, in the agent-based model we propose (4) and (5), the cell i has no knowledge of the velocity of any of its neighbors, i.e. ω j is not used to define the evolution of (x i , ω i ). Therefore, it is unclear at first why the dynamics proposed could generate similar flocking patterns.
To address this question, we provide a linear perturbation analysis of the model with respect to the eccentricity in the case of only two cells i and j. Let's denote θ ij the angle between the direction of the cell ω i and the relative position vector x j − x i (see Fig 6). Thus, one can

PLOS COMPUTATIONAL BIOLOGY
Self-organization in brain tumors: Influences of cell morphology and cell density write ω i = (cos θ ij , sin θ ij ) T using the basis fo i ; o ? i g. In particular, We deduce that: Moreover, elementary geometry shows that x ij = | x j − x i | cos θ ij and y ij = |x j − x i | sim θ ij and thus: Notice that C � 0 since F 0 ðr 2 ij Þ � 0. As a consequence, if i and j stay close enough (e.g. r 2 ij < 1), there are two stable equilibria for θ ij at ±π/2. Sketching the phase portrait in Fig 6 indicates that ω i will rotate toward the stable equilibrium to align with (x j − x i ) ? . By a similar argument, ω j will be orthogonal to (x j − x i ) as well. Therefore, instead of a direct alignment between ω i and ω j , we have an indirect alignment: both vectors will align with (x j − x i ) ? . Notice that this indirect form of alignment allows for the two vectors ω i and ω j to be negatively aligned, i.e. ω i = −ω j which could lead to streaming formation. In dimension larger 2, (x j − x i ) ? is an hyperplane, thus even if ω i and ω j become orthogonal to (x j − x i ) ? , it is insufficient to conclude that ω i and ω j will become parallel. Indeed, as we will show numerically in the next section, one has to consider also the spatial configuration (e.g. density) to predict whether the dynamics will generate flock or stream formations.

Density effect
In the previous section, we investigated how cell morphology (i.e. a, b) promotes the emergence of flocking patterns (i.e. cells moving in the same direction). Our formal analyses show that we could also observe stream formation (i.e. cells moving in opposite directions). In this section, we will define the conditions under which streams emerge. To define the conditions which allow the emergence of streams we will study the dynamics of our system as we vary the parameters α (strength repulsion), β (strength steering) and N (density). We will fix the shape of the cells with a = 5.5μm and b = 3μm as they are the most common values experimentally. The range of the parameters are given in Table 2.

PLOS COMPUTATIONAL BIOLOGY
Self-organization in brain tumors: Influences of cell morphology and cell density steering). In Fig 7, we illustrate snapshots of three simulations where we increase the density from N = 1000 to N = 2000 at t = 1000 unit time. Cells are color coded by orientation: we determine the nematic average direction O nem (see S1 Text) and then color each We notice that when the number of particles is low with N = 1000 (Fig 7a), almost all the cells are perfectly aligned in the same direction leading to a flocking configuration. The evolution of the polarization ψ given in Fig 8 confirms this observation since ψ becomes close to 1 for N = 1000. As we increase the density with N = 1500 (Fig 7b), we observe that the number of cells moving in the opposite direction becomes larger making the polarization decay to only ψ = .2. Finally in the case where N = 2000 (Fig 7c), the number of cells moving in opposite direction becomes balanced and we observe the formation of a stream where the flow inside the domain is bidirectional. Indeed, the polarization ψ is close to zero for N = 2000 where the nematic polarization γ is around.9 (see S2 Video).
Local minimum for the energy. We conclude that increasing the density of cells is the underlying mechanism for stream formation. However, one has to notice that we always use as initial configuration random configurations for the velocities ω i . If one would start from a

PLOS COMPUTATIONAL BIOLOGY
Self-organization in brain tumors: Influences of cell morphology and cell density perfect flock with no overlapping (i.e. ω i = O for all i and r ij > 1 for all i, j), then the configuration will remain in this configuration, it will be simply transported with a constant velocity. Thus, we will not observe the formation of a stream even at large cell density (e. g. N = 2000). In other words, the flocking configuration can be a globally stable configuration. In contrast, in a stream formation, the cells at the border between regions moving in opposite direction are in an unstable equilibrium (see Fig 6). Therefore, if the steering coefficient β remains small enough, the streaming configuration will be maintained, the non-overlapping physical constraint (through the repulsion α) prevents the cells from turning.
Another formal justification for the emergence of the stream configuration comes from the total potential energy V (1): with r ij given by (1). The dynamics (4) and (5) is a gradient descent of the potential V (i.e. The gradient descent decays the potential V whereas the free transport could either increase or decrease V. But as we increase α, the dynamics become more likely to become fixed in a local equilibrium (i.e. stream). Perturbations to the free transport component of the dynamics will be insufficient to move the configuration away from a local equilibrium (see Fig 9). However, on a large time scale, it is still possible that a stream configuration would eventually become a flock. The reverse situation, a flock becoming a stream, is unlikely as it would require an increase in the potential energy V. Since the dynamics is not conservative, we cannot rule out this scenario but numerically we haven't observed such transition.
Phase diagram. We have identified two configurations: flocking when the cells are aligned (i.e. ψ � 1, γ � 1), stream when the cells are moving in opposite directions (i.e. γ � 1). The convergence of the dynamics toward one of these configurations depends on the density N  (Figs 7 and 8) and on the parameters α and β. We would like to study the effects of repulsion and alignment (i.e the coefficients α and β resp.) for the three distinct cases of N.
With this aim, we fix the shape of the cells (a = 5.5μm and b = 3μm) and make a systematic analysis by varying continuously α 2 [0, 200] and β 2 [0, 10]. For each value of α and β, we perform 5 simulations and compute several statistics after t = 1000 unit of time. For instance, in Fig 10 (top-left) we plot the polarization ψ depending on α and β in the case N = 1000. The scatter plot represents all the data points (α j , β j , ψ j ). Notice that for a given set (α, β) we find varying polarization ψ due to random initial conditions. We represent the average polarization hψi as a surface computed from the 5 final configurations with similar parameters (α, β). To reduce the fluctuation, we estimate a local averaging ('loess') in Fig 10 (top-right) which makes possible the estimation of a smooth region where the polarization ψ is higher than a given threshold. Regarding the use of averaging or smoothing, we observe that the polarization is surprisingly small when β is large and α small. There is also another region where the polarization decays when α is large and β is small.
A better visualization is to draw the polarization ψ as a heat-map depending on α and β ( Fig  10 center-left). Through the use of smoothing, we can also estimate contours at different thresholds (ψ = .5 and ψ = .8). We proceed similarly with the nematic polarization γ (7) in Fig  10 center-right. We notice that in contrast to the polarization ψ, the nematic polarization remains large even when β is small and α is large: indeed, this is the regime where we observe the formation of streams.
Finally, we also use a third statistic to characterize the configuration using the average cluster size � jCj (10). We use the radius R = 10μm to define the clusters (i.e. two cells are connected if their distance kx j − x i k is less than 10μm). The average size cluster � jCj is then estimated in Fig 10 (bottom). We observe two regions: cluster sizes are (relatively) smaller when α is small and independent of β. Thus, the repulsion α governs the formation of clusters.
We combine the three statistics (polarization ψ, nematic polarization γ, cluster size) to create a phase diagram in the parameter space (α, β). Three regions are delimited: The results are given in Fig 11. For most of the parameters α and β, the dynamics converge to a flock.
Performing a similar investigation for N = 1500 and N = 2000 lead to drastically different results. The regions where flocking occurs are more narrow (Fig 12a). But surprisingly stream formation is still occurring for all values of α as long as the steering coefficient β is small enough (Fig 12b). Only the cluster formation through the statistic � jCj remains similar (see Fig  12c) as in the case N = 1000. As a result, the phase diagrams for N = 1500 and N = 2000 contain a large region not identifiable as either flock or stream (Fig 13). Notice that increasing density does not penalize the formation of streams in the region where β is small and α is large.

Dynamics in 3D
Finally, we would like to study the dynamics (4) and (5) in R 3 . There are several key differences between R 2 and R 3 for the dynamics. Our formal discussion in see section 1 showed that the dynamics enforce that nearby cells (denoted i and j) must have their velocity (ω i and ω j ) orthogonal to their relative position (x j − x i ). In R 2 , we concluded that nearby cells must be aligned at equilibrium, i.e. ω i = ω j or ω i = −ω j . This is no longer the case in R 3 : ω i and ω j could The estimation has been smoothed using local regression ('loess'). We then deduce the region when the cluster size is below a certain threshold. be orthogonal to (x j − x i ) without being aligned. Therefore, it is uncertain if one should observe the emergence of flock or stream formations in R 3 .
Another difference in R 3 is that the nematic polarization γ is no longer defined as a velocity vector ω in S 2 is defined using two angles instead of one. However, we have provided an alternative quantity denoted J in S1 Text. We also use a smaller domain for our simulation in order to keep the same density as in R 2 (maintaining a similar ratio "volume occupied/volume domain"); thus we reduce the size of the box to L = 70μm. Otherwise, the other parameters remain those of the previous simulations (see Table 3), in particular we use α = 100 and β = .1 to be in the regime more susceptible of stream formation (at least in R 2 ).
First, we investigate the model with N = 1000 cells (low density). We plot the final configuration at t = 1000 unit times starting from two different initial conditions in Fig 14. As in Fig 7, we color code the cells depending on their orientation to help visualize cells moving in opposite directions. Notice that cells do not necessarily move parallel to each other (they can move orthogonal to each other). But after a transient period, only one or two directions remain as cells form either a flock or a stream. Indeed, we observe the formation of a flock (Fig 14-top  left) and of a stream (Fig 14-top right). Note that even when the flock develops (top left), few isolated cells (red) are still moving in opposing direction to to the main flow (blue). Thus, flock and stream configurations can emerge when the cell density is low.
The situation is different when we increase the density to N = 1500 and N = 2000. In this case we only observe the formation of streams (see Fig 14-bottom). Similar to the situation in R 2 , increasing the density reduces the possibility for the cells to rotate and therefore streams are more likely to occur. Plotting the time evolution of both the polarization ψ and nematic

Discussion
Whether self-organization exists in brain tumors is incompletely understood. The identification of multicellular structures in brain tumors suggest that single cells are able to coalesce into multicellular patterns, possibly as a result of self-organization. We have recently described such multicellular structures within gliomas [34]. As these structures are reminiscent of streams described in other systems, we have labeled these structures oncostreams. Oncostreams extend over 100-500μm long and 50-200μm wide, and contain elongated cells. In this work we aimed to understand the role of the elongated morphology of cells within such oncostreams in the formation of the multicellular patterns, and whether these patterns are the result of self-organization operating within gliomas. We propose an agent-based model that describes the motion of cancer cells and the emergence of pattern formation within gliomas. In our model, the morphology of the cells plays a key role in glioma pattern formation since cell eccentricity allows the cells to align (indirectly) PLOS COMPUTATIONAL BIOLOGY Self-organization in brain tumors: Influences of cell morphology and cell density to each other and eventually coalesce to form a flock or a stream. In the special case where cells are circular, cells cannot align and thus no flock or stream can formed. The emergence of such multicellular patterns is also governed by additional parameters, i.e. α (repulsion), β (steering), and cell density. Several phase diagrams summarizing the effects of these parameters have

PLOS COMPUTATIONAL BIOLOGY
Self-organization in brain tumors: Influences of cell morphology and cell density been estimated for various densities. In contrast to mean-field type models, the density drastically changed the dynamics. Flocking configurations became more sparse and streams' density increased as cellular density increased. This has important biological implications as the density of gliomas is constantly changing. We have discovered that glioma oncostreams indeed display characteristics of self-organization. Namely, our data strongly suggest that the emergence of the large-scale structures within brain tumors, flocks and streams, result from intercellular interactions. As the density of flocks and streams correlates with glioma aggressiveness, we can link the macroscopic behavior of brain tumors to the intercellular interactions of individual glioma cells. The formation of large scale patterns that result from intercellular interactions, and the fact that these structures determine glioma behavior, i.e., growth, invasion, aggressiveness, allow us to conclude that gliomas display clear evidence of self-organization, and that tumor self-organization plays an important role in tumor malignity.
Most of our experimental work was performed using mouse tissues and genetically engineered mouse models of brain tumors. However, the recognition of oncostreams in human malignant gliomas [35] strongly suggests that human brain tumors can also self-organize into structures that influence tumor malignity. We suggest that a novel approach to the treatment of brain tumors was to disassemble the oncostreams. Ongoing molecular studies of oncostreams indicate that this is feasible. Our data make important biological predictions.
Our model achieves flocking through cells of an eccentricity larger than 1, and parameters which determine cell repulsion, and cell steering. No function of cell adhesion is included into the model. This strongly suggests that the molecular intercellular interactions leading to stream and flock formation may regulate the degree of intercellular adhesion. It will be important to investigate the molecular basis of stream and flock formation in gliomas. The cytoskeleton is also likely to play a central role in stream and flock formation as the eccentricity has an optimal value to form oncostreams, and that, once exceeded, the capacity to form oncostreams decreases. Molecular mechanisms that might optimize cell eccentricity are currently not understood, and might yield important results concerning the molecular mechanisms that are necessary to form oncostreams and flocks. Our data also suggest that intra-glioma dynamics are cell density dependent, as the directionality of the cells changes significantly upon increases of cell density. As intraglioma dynamics are important to tumor growth and invasion, changes achieved by therapeutic cytotoxic agents, may not eliminate glioma growth, but rather alter its organization and directionality. As methods become available to affect the overall organization of gliomas, we suggest that the effects on intratumoral dynamics be considered in terms of drugs' mechanisms of action. Increasing density is also likely to affect the overall macroscopic tumor growth, as flocks (which move in one direction), are able to convert into streams (which move in two directions). As a consequence, the quality and distribution of tumor growth may change in response to partial cytotoxicity to a diffuse tumor capable of growing in more directions as it invades surrounding normal brain.
The phase diagrams indicate the potential existence of critical points and phase transitions at which non-aligned cells can become aligned and form a flock or a stream, a relationship which is also dependent on glioma cell density. The existence of critical points will be important as they might regulate the sudden scattering of cells, or their organization into patterns of collective motion that are likely to determine tumor growth.
Our work proposes a number of extensions which will be pursued in the future. For instance, it will be important to mix cells with different shapes (i.e., with different values for α and β) since not all the cells are identical (see for instance [36]). We will also study how cell eccentricity varies over time and even whether it influences cell mitosis, and/or the birth/death cycles. Increasing the density is also challenging numerically as the dynamics become singular when two cells overlap which is more likely with a birth process. To avoid this complication, one could explore a continuous description of the dynamics through a Partial Differential Equation (PDE) [30][31][32]. Such PDE description might provide some hindsight about the emergence of flock or stream in certain regimes (e.g. α � 1, β � 1).
Adding cell divisions will also raise new questions such as how fast do cell spread depending on their distance to the center of the tumor. It is yet to be determined whether cells will move faster close to the center (large density) or at the periphery of the tumor (low density).
Another extension of the model would consist in adding a "contact inhibition of locomotion" (CIL) to the cells [37]. The main idea is that cells would reduce their self-propulsion as they experience contact. As a result, we would expect that the perturbation due to the free transport component (see Fig 9) would be reduced and therefore flocks and streams would be more likely to occur at larger density. Such behavior would be consistent with experimental observations where 2D-cell layer clusters have been observed at large density [38].
In summary, through a detailed investigation of patterns of glioma growth and agent-based mathematical modeling we explain the importance of cell shape during glioma growth, and its consequences for glioma self-organization, aggressiveness and invasion. The long term consequences of glioma self-organization will impact our understanding of glioma biology, and suggest novel treatments.