Skip to main content
Advertisement
  • Loading metrics

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

Abstract

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 . 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 with the noticeable difference that stream formation is more ubiquitous compared to flock formation.

Author summary

Self-organization is the formation of large-scale multicellular patterns that result exclusively from the interactions amongst constituent single cells. To establish the existence of self-organization in brain tumors we used agent-based modeling based on data extracted from static and dynamic genetically engineered and implantable mouse glioma models. Implementation of our model in identifies the dynamics that lead to formation of flocks (cells moving in a single direction), streams (cells moving in two directions), and cells moving as swarms or scattering. Increasing cellular density reduced formation of flocks and increased the formation of streams both in and in . These results demonstrate the detailed mechanism leading to self-organization in brain tumors. As increasing density of oncostreams correlates with tumor malignancy, we establish a pathophysiological link between self-organization of glioma tumors and glioma malignancy. We propose the dismantling of oncostreams as a new therapeutic approach to the treatment of brain tumors.

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 [14]. 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 [1217]). 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 [1820].

thumbnail
Fig 1. Representative areas of elongated cells versus rounded cells found in glioma tumors.

Hematoxilin and eosin staining of a genetic engineered mouse glioma tumor expressing the following genetic lesions: Nras overexpression, shp53 downregulation, shATRx downregulation). A) Outlines and arrows demarcate multicellular structures formed by elongated cells (*) from areas of rounded cells (+). Scale bar 50μm. B) Image magnification of (A) highlighting in black broken lines the morphological differences between elongated (red) and rounded cells (blue). Scale bar 20μm.

https://doi.org/10.1371/journal.pcbi.1007611.g001

We propose an agent-based model that utilizes two mechanisms: i) self-propulsion, ii) cell-cell avoidance due to non-overlapping constraints. Since the cells have an ellipsoid shape, cell-cell 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 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 since cells avoiding each other no longer move aligned or in opposite direction as in . However, our simulations show that flock and stream formation do still occur in 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 [2427]. 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 [3032] 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 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 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 with d the spatial dimension (d = 2 or 3), moving with velocity i where c > 0 is the speed (supposed constant) and 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 Fig 2-left). 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 Vi between cells that measures the tension exerted on cell i generated by the surrounding cells: (1)

thumbnail
Fig 2.

Left: a cell i is described by its position xi, orientation ωi and its elliptic shape determined by the two morphological components a and b. Right: function Φ relies spacing between cell rij (1) into tension that generates repulsion when two cells touch each other.

https://doi.org/10.1371/journal.pcbi.1007611.g002

The quantity rij is referred to as the normalized distance between the centers of the cells i and j. For instance, if a = b we recover that rij is simply the norm ∥xjxi∥/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 (i.e. d = 2) and can be generalized to by defining rij as follows: (2) where e ∈ (0, 1) is the eccentricity of an ellipse defined as .

To prevent overlapping, the function Φ has to be singular at the origin. We choose the following smooth function with compact support (see Fig 2-right): (3)

As rij decreases, Φ increases resulting into repulsion. We have now defined all the quantities required to define our agent-based model.

Definition 1 Consider for i = 1..N and the dimension space d = 2 or d = 3. The self-propelled dynamics are defined as: (4) (5) where α, β and c are positive constant, Vi is the tension defined in (1) and is 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 Vi (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. . Indeed, in that case, turning will have no effect on the reduction of the tension Vi. Thus, steering effects can only occur when if ab.

Remark 2 Notice that rij cannot be considered a distance between cells i and j as it is not symmetric (i.e. rijrji in general). Thus, the influence of the cell i on j is in general different from the cell j on i, i.e .

Results

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.

thumbnail
Table 1. Parameters used for the simulations of Figs 3 and 5.

https://doi.org/10.1371/journal.pcbi.1007611.t001

We perform the simulation in a square domain Ω = [0, L]×[0, L] with periodic boundary conditions. For the initial condition, the positions of the N particles {xi}{i = 1…N} are distributed uniformly Ω and their directions {ωi}{i = 1…N} are taken randomly on the unit circle . In Fig 3, we draw the final configuration of the dynamics after T = 1000 unit time for two types of cells: circular shape (a = b = 4μm, e = 0) and elliptic shape (a = 5.5μm, b = 3μm, e = .84). We observe that circular cells have no particular spatial organization (Fig 3-left) while elliptic cells have formed clusters moving in the same directions (Fig 3-right). The full simulation is also available (see S1 Video).

thumbnail
Fig 3. Snapshot of the simulation of the model starting from a uniform distribution.

After t = 1000 unit of time, circular cells (a = b = 4μm, e = 0) do not form any flocking pattern (left) whereas elliptic cells (a = 5.5μm, b = 3μm, e = .84) move in a common direction (right).

https://doi.org/10.1371/journal.pcbi.1007611.g003

Statistical characterization.

To further investigate the dynamics, we introduce several statistics to characterize the emergent behavior.

Definition 3 Consider a velocity distribution . We denote by ψ the polarization: (6)

Similarly, we define the nematic polarization [12]: (7) where θi is the angle between the direction ωi and the horizontal axis and 〈〉 denotes the averaging over the indices i (i.e. ).

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).

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 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. ij) with ij defined: (8)

Clusters are defined as the connected components for this relationship: two cells i0 and j0 belong to the same cluster if there exists particles i1, …, ik (a path) such that (9)

The cluster size denotes the number of cells belonging in the cluster k. The average cluster size is defined as the expected cluster size over all the positions: (10) where denotes the cluster containing the cell i.

We illustrate the three statistics used in Fig 4.

thumbnail
Fig 4. The statistics used to characterize the dynamics: The polarization ψ (6), the nematic polarization γ (7) and the clustering (8) and (9).

https://doi.org/10.1371/journal.pcbi.1007611.g004

In Fig 5left, 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.

thumbnail
Fig 5. Polarization ψ (6) over time while varying the eccentricity of the cell e.

Ellipsoid cells that align will lead to an increase of ψ close to its maximum 1. For the circular cells (blue curve), the polarization remains very low as no streams emerge from the dynamics (left). The polarization ψ over eccentricity e during the final time t = 1000 of the left figure will form a parabola. By increasing the eccentricity, there is no fundamental impact on the polarization coefficient of the cells (right).

https://doi.org/10.1371/journal.pcbi.1007611.g005

To further investigate the relationship between polarization and eccentricity, we plot in Fig 5right 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 (xi, ω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 xjxi (see Fig 6). Thus, one can write ωi = (cos θij, sin θij)T using the basis . In particular,

thumbnail
Fig 6. Indirect alignment of two cells i and j.

Both cell i and j will rotate to align with the orthogonal vector to xjxi.

https://doi.org/10.1371/journal.pcbi.1007611.g006

We deduce that: where . Moreover, elementary geometry shows that xij = |xjxi| cos θij and yij = |xjxi| sim θij and thus: (11)

Notice that C ≤ 0 since . As a consequence, if i and j stay close enough (e.g. ), 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 (xjxi). By a similar argument, ωj will be orthogonal to (xjxi) as well.

Therefore, instead of a direct alignment between ωi and ωj, we have an indirect alignment: both vectors will align with (xjxi). 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, (xjxi) is an hyperplane, thus even if ωi and ωj become orthogonal to (xjxi), 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.

Emergence of streams.

To illustrate the formation of streams (see Eqs (4) and (5)), we perform simulations within the parameter constraints of: α = 100, β = .1 (strong repulsion, low 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 Ωnem (see S1 Text) and then color each cell i blue if 〈ωi, Ωnem〉>0 or in red if 〈ωi, Ωnem〉<0.

thumbnail
Fig 7. Snapshots of the simulation of the dynamics at t = 1000 unit time for various cell densities (N = 1000, 1500 and 2000).

Red cells are moving in opposite direction to the blue cells. Flocking appears when the density is low (A) but then the dynamics start to converge to stream formation as we increase the density (B-C).

https://doi.org/10.1371/journal.pcbi.1007611.g007

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).

thumbnail
Fig 8. Polarization ψ and nematic polarization γ for the simulations of Fig 7.

A flock occurs at low density (i.e. N = 1000) where ψ and γ converge approximately to 1, whereas streams emerge at larger density (i.e. N = 1500 and N = 2000).

https://doi.org/10.1371/journal.pcbi.1007611.g008

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 perfect flock with no overlapping (i.e. ωi = Ω for all i and rij > 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): (12) with rij given by (1). The dynamics (4) and (5) is a gradient descent of the potential V (i.e. ) combined with a free transport component (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.

thumbnail
Fig 9. Sketch representation of the potential energy V (12) over the configuration space {(xi, ωi)}i = 1,…,N.

Stream can be seen as local equilibrium whereas flock are global equilibrium. When the parameter α (repulsion) is increased, the stream configuration become more stable.

https://doi.org/10.1371/journal.pcbi.1007611.g009

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 α ∈ [0, 200] and β ∈ [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 〈ψ〉 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.

thumbnail
Fig 10.

Top-left: for each pair α, β, we estimate the polarization ψ (scatter point) at the end of 5 simulations. The average is then computed to construct a surface plot. Top-right: we use a local regression (’loess’) to estimate ψ as a function of α and β reducing the fluctuation. Center-left: heat-map representation of the polarization ψ as a function of a, b using the smooth estimation of ψ. The contour ψ = .8 will be used to determine the region when the dynamics generate flock. Center-right: we perform a similar analysis as the left figure using the nematic polarization γ. Bottom: the average size cluster (10) for various values of α and β. The estimation has been smoothed using local regression (’loess’). We then deduce the region when the cluster size is below a certain threshold.

https://doi.org/10.1371/journal.pcbi.1007611.g010

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 (10). We use the radius R = 10μm to define the clusters (i.e. two cells are connected if their distance ∥xjxi∥ is less than 10μm). The average size cluster 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:

  1. i). flocking: {a, b such that ψ >.8},
  2. ii). streaming: {a, b such that γ >.7 and ψ <.8},
  3. iii). scattering: {a, b such that }.

The results are given in Fig 11. For most of the parameters α and β, the dynamics converge to a flock.

thumbnail
Fig 11. Combining the results of Fig 10, we create a phase diagram consisting of three regions for the configuration: Flocking (ψ >.8), streams (γ >.7 and ψ <.8) and scattering ().

https://doi.org/10.1371/journal.pcbi.1007611.g011

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 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.

thumbnail
Fig 12.

A) Average polarization ψ for various parameters α and β with N = 1500 and N = 2000. Notice that the polarization is significantly smaller compared to the case N = 1000 (Fig 10). B) Nematic polarization γ for N = 1500 and N = 2000. γ remains close to 1 for any values of α when β is small. C) The average cluster size behave similarly as in the case N = 1000 with smaller clustering for small value of α.

https://doi.org/10.1371/journal.pcbi.1007611.g012

thumbnail
Fig 13. Phase diagram when the total number of cells N is 1500 (left) and 2000 (right).

As we increase the density, the regions for flocking configurations drastically reduce. However, streams are still form when β is small.

https://doi.org/10.1371/journal.pcbi.1007611.g013

Dynamics in 3D

Finally, we would like to study the dynamics (4) and (5) in . There are several key differences between and 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 (xjxi). In , we concluded that nearby cells must be aligned at equilibrium, i.e. ωi = ωj or ωi = −ωj. This is no longer the case in : ωi and ωj could be orthogonal to (xjxi) without being aligned. Therefore, it is uncertain if one should observe the emergence of flock or stream formations in .

Another difference in is that the nematic polarization γ is no longer defined as a velocity vector ω in 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 (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 ).

thumbnail
Table 3. Parameters used for the simulations in (Figs 14 and 15).

https://doi.org/10.1371/journal.pcbi.1007611.t003

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 , 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 polarization J in Fig 15 confirms our observations. The quantity J always converges to 1 whereas the nematic polarization ψ stays low except when N = 1000.

thumbnail
Fig 14. Snapshots of simulation in at t = 1000 unit time with number of particles N equal to 1000 (top), 1500 and 2000 (bottom).

Flocking and streaming appear when the number of particles is low depending on the initial condition (top left and top right respectively). Whereas only stream emerges when the number of particles is higher, N equal 1500 and 2000 (bottom—left and right respectively). We color code the cells in blue or red depending on the direction in comparison to the nematic average (see S1 Text).

https://doi.org/10.1371/journal.pcbi.1007611.g014

thumbnail
Fig 15. The evolution of the polarization ψ and the nematic polarization J for the simulations presented in Fig 14.

When the number of particle is low (N = 1000), flock and stream emerge depending on the initial condition leading to an average value for ψ is around 0.6 and J close to 1. However, for larger density (N = 1500 and N = 2000), only streams emerge since the low polarization ψ is low and the nematic polarization J is close to 1.

https://doi.org/10.1371/journal.pcbi.1007611.g015

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) 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 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) [3032]. 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.

Supporting information

S1 Text. Explicit expression of the model and definition of nematic average.

https://doi.org/10.1371/journal.pcbi.1007611.s001

(PDF)

S1 Video. Circles VS ellipses.

Numerical simulation comparing the dynamics with circles (i.e. a = b = 4μm) and ellipses (a = 5.5μm, b = 3μm). Fig 3 corresponds to a screenshot of this video taken at t = 1000 unit time.

https://doi.org/10.1371/journal.pcbi.1007611.s002

(AVI)

S2 Video. Stream formation.

At large density, stream formations are more likely to occur (see also Fig 7). After a transient time, the dynamics generate streams that move in opposite direction. Cells have been colored depending on their direction: blue cells are the one moving to the right, red cells move left and white cells move up or down.

https://doi.org/10.1371/journal.pcbi.1007611.s003

(AVI)

References

  1. 1. Deisboeck T. and Couzin I. D. Collective behavior in cancer cell populations. Bioessays, 31(1):190–197, 2009. pmid:19204991
  2. 2. Ribba B., Alarcón T., Marron K., Maini P., and Agur Z. The use of hybrid cellular automaton models for improving cancer therapy. Cellular Automata, p.444–453. Springer, 2004.
  3. 3. Anderson A. and Quaranta V. Integrative mathematical oncology. Nat Rev Cancer, 8(3):227–234, March 2008. pmid:18273038
  4. 4. Harpold H., Alvord E., and Swanson K. The Evolution of Mathematical Modeling of Glioma Proliferation and Invasion. Journal of Neuropathology and Experimental Neurology, 66(1):1–9, January 2007. pmid:17204931
  5. 5. Moussaïd M., Guillot E., Moreau M., Fehrenbach J., Chabiron O., Lemercier S., Pettré J., Appert-Rolland C., Degond P., and Theraulaz G. Traffic Instabilities in Self-Organized Pedestrian Crowds. PLoS Comput Biol, 8(3):e1002442, March 2012. pmid:22457615
  6. 6. Ballerini M., Cabibbo N., Candelier R., Cavagna A., Cisbani E., Giardina I., Lecomte V., Orlandi A., Parisi G., and Procaccini A. Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study. Proceedings of the National Academy of Sciences, 105(4):1232, 2008.
  7. 7. Couzin I.D., Krause J., Franks N.R., and Levin S.A. Effective leadership and decision-making in animal groups on the move. Nature, 433(7025):513–516, 2005. pmid:15690039
  8. 8. Rabani A., Ariel G., and Be’er A.. Collective Motion of Spherical Bacteria. PLos One, 8(12):e83760, December 2013. pmid:24376741
  9. 9. Ariel G., Sidortsov M., Ryan S., Heidenreich S., Be’er M., and Bär A. Collective dynamics of two-dimensional swimming bacteria: Experiments and models. Physical Review E, 98(3):032415, 2018.
  10. 10. Palsson E. and Othmer H. A model for individual and collective cell movement in Dictyostelium discoideum. Proceedings of the National Academy of Sciences, 97(19):10448–10453, 2000.
  11. 11. Menzel A. and Ohta T. Soft deformable self-propelled particles. EPL (Europhysics Letters), 99(5):58001, 2012.
  12. 12. Li X., Balagam R., He T-F., Lee P., Igoshin O., and Levine H. On the mechanism of long-range orientational order of fibroblasts. Proceedings of the National Academy of Sciences, 114(34):8974–8979, 2017.
  13. 13. Chen X., Dong X., Be’er A., Swinney H., and Zhang H. P. Scale-invariant correlations in dynamic bacterial clusters. Physical review letters, 108(14):148101, 2012. pmid:22540824
  14. 14. Pletjushkina O., Ivanova O., Kaverina I., and Vasiliev J. Taxol-treated fibroblasts acquire an epithelioid shape and a circular pattern of actin bundles. Experimental cell research, 212(2):201–208, 1994. pmid:7910561
  15. 15. Xiong Y., Rangamani P., Fardin M-A., Lipshtat A., Dubin-Thaler B., Rossier O., Sheetz M., and Iyengar R. Mechanisms controlling cell size and shape during isotropic cell spreading. Biophysical journal, 98(10):2136–2146, 2010. pmid:20483321
  16. 16. Balagam R., Litwin D., Czerwinski F., Sun M., Kaplan H., Shaevitz J., and Igoshin O. Myxococcus xanthus gliding motors are elastically coupled to the substrate as predicted by the focal adhesion model of gliding motility. PLoS computational biology, 10(5):e1003619, 2014. pmid:24810164
  17. 17. Ariel G., Shklarsh A., Kalisman O., Ingham C., and Ben-Jacob E.. From organized internal traffic to collective navigation of bacterial swarms. New J. Phys., 15(12):125019, December 2013.
  18. 18. Börzsönyi T. and Stannarius R. Granular materials composed of shape-anisotropic grains. Soft Matter, 9(31):7401–7418, 2013.
  19. 19. Wegner S., Börzsönyi T., Bien T., Rose G., and Stannarius R. Alignment and dynamics of elongated cylinders under shear. Soft Matter, 8(42):10950–10958, 2012.
  20. 20. Baule A. and Makse H. Fundamental challenges in packing problems: from spherical to non-spherical particles. Soft Matter, 10(25):4423–4429, 2014. pmid:24898797
  21. 21. Cucker F. and Smale S. Emergent Behavior in Flocks. IEEE Transactions on automatic control, 52(5):852, 2007.
  22. 22. Vicsek T., Czirók A, Ben-Jacob E., Cohen I., and Shochet O. Novel type of phase transition in a system of self-driven particles. Physical Review Letters, 75(6):1226–1229, 1995. pmid:10060237
  23. 23. Cho H., Jönsson H., Campbell K., Melke P., Williams J., Jedynak B., Stevens A., Groisman A., and Levchenko A. Self-organization in high-density bacterial colonies: efficient crowd control. PLoS biology, 5(11):e302, 2007. pmid:18044986
  24. 24. Lowengrub J., Frieboes H., Jin F., Chuang Y-L., Li X., Macklin P., Wise S., and Cristini V. Nonlinear modelling of cancer: bridging the gap between cells and tumours. Nonlinearity, 23(1):R1–R91, January 2010. pmid:20808719
  25. 25. Saut O., Lagaert J-B., Colin T., and Fathallah-Shaykh H. A multilayer grow-or-go model for GBM: effects of invasive cells and anti-angiogenesis on growth. Bulletin of mathematical biology, 76(9):2306–2333, 2014. pmid:25149139
  26. 26. Deisboeck T., Wang Z., Macklin P., and Cristini V. Multiscale Cancer Modeling. Annu Rev Biomed Eng, 13, August 2011. pmid:21529163
  27. 27. Degond P., Peurichard D. Modelling Tissue Self-Organization: From Micro to Macro Models. Multiscale Models in Mechano and Tumor Biology, p.93–108, Springer, 2017. https://doi.org/10.1007/978-3-319-73371-5_5
  28. 28. Swanson K., Bridge C.J. Murray D., and Alvord E. Jr. Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion. Journal of the Neurological Sciences, 216(1):1–10, 2003. pmid:14607296
  29. 29. Hawkins-Daarud A., Prudhomme S., Zee K., and Oden J. Bayesian calibration, validation, and uncertainty quantification of diffuse interface models of tumor growth. J. Math. Biol., 67(6-7):1457–1485, December 2013. pmid:23053536
  30. 30. Morale D., Capasso V., and Oelschläger K.. An interacting particle system modelling aggregation behavior: from individuals to populations. J. Math. Biol., 50(1):49–66, January 2005. pmid:15692840
  31. 31. Filbet F., Laurençot P., and Perthame B. Derivation of hyperbolic models for chemosensitive movement. Journal of Mathematical Biology, 50(2):189–207, 2005. pmid:15480673
  32. 32. Peruani F., Deutsch A., and Bär M.. A mean-field theory for self-propelled particles interacting by velocity alignment mechanisms. The European Physical Journal Special Topics, 157(1):111–122, 2008.
  33. 33. Baker G., Yadav V., Motsch S., Koschmann C., Calinescu A., Mineharu Y., Camelo-Piragua S., Orringer D., Bannykh S., Nichols W., deCarvalho A., Mikkelsen T., Castro M., and Lowenstein P. Mechanisms of glioma formation: iterative perivascular glioma growth and invasion leads to tumor progression, VEGF-independent vascularization, and resistance to antiangiogenic therapy. Neoplasia, 16(7):543–561, 2014. pmid:25117977
  34. 34. Lowenstein P., Motsch S., Yadav V., Zamler D., Koschmann C., F. Nunez-Aguilera, Calinescu A., Kamran N., Dzaman N., MG. Castro. Streams, swirls, and neurospheres: in vivo self-organization of brain tumors revealed by mathematical and biological modeling. American Society for Cell Biology (ASCB) Annual Meeting San Diego, California, 12-16, 2015.
  35. 35. Comba A., Dunn P., Argento A.E., Kadiyala P., Motsch S., Kish P., Kahana A., Zamler D., Muraszko K., Castro M.G., and Lowenstein P.R. Oncostreams: novel dynamics pathological multicellular structures involved in glioblatoma growth and invasion. Journal of Clinical and Translational Science 3(s1):111, 2019.
  36. 36. Kansal A. R., Torquato S., Chiocca E. A., and Deisboeck T. S. Emergence of a subpopulation in a computational model of tumor growth. Journal of Theoretical Biology, 207(3):431–441, 2000. pmid:11082311
  37. 37. Davis J., Luchici A., Mosis F., Thackery J., Salazar J., Mao Y., Dunn G., Betz T., Miodownik M., and Stramer B. Inter-cellular forces orchestrate contact inhibition of locomotion. Cell, 161(2), 361–373, 2015. pmid:25799385
  38. 38. Deisboeck T., Berens M., Kansal A., Torquato S., Stemmer-Rachamimov A., and Chiocca E. Pattern of self-organization in tumour systems: complex growth dynamics in a novel brain tumour spheroid model. Cell proliferation, 34(2), 115–134, 2001. pmid:11348426