article: Oligodendrocyte precursor cells (OPCs) have remarkable properties: they represent the most abundant cycling cell population in the adult normal brain and they manage to achieve a uniform and constant density throughout the adult brain. This equilibrium is obtained by the interplay of four processes: division, differentiation or death, migration and active self-repulsion. They are also strongly suspected to be at the origin of gliomas, when their equilibrium is disrupted. In this article, we present a model of the dynamics of OPCs, first in a normal tissue. This model is based on a cellular automaton and its rules are mimicking the ones that regulate the dynamics of real OPCs. The model is able to reproduce the homeostasis of the cell population, with the maintenance of a constant and uniform cell density and the healing of a lesion. We show that there exists a fair quantitative agreement between the simulated and experimental parameters, such as the cell velocity, the time taken to close a lesion, and the duration of the cell cycle. We present three possible scenarios of disruption of the equilibrium: the appearance of an over-proliferating cell, of a deadless/non-differentiating cell, or of a cell that lost any contact-inhibition. We show that the appearance of an over-proliferating cell is sufficient to trigger the growth of a tumor that has low-grade glioma features: an invasive behaviour, a linear radial growth of the tumor with a corresponding growth velocity of less than 2 mm per year, as well a cell density at the center which exceeds the one in normal tissue by a factor of less than two. The loss of contact inhibition leads to a more high-grade-like glioma. The results of our model contribute to the body of evidence that identify OPCs as possible cells of origin of gliomas.
  Contrary to a long-lasting belief, there exists a population of proliferating cells everywhere in the mammalian adult brain [1,2]. These cells are the oligodendrocyte precursor cells (OPCs) and they differentiate into oligodendrocytes [3], but also in some circumstances to astrocytes. They can be found everywhere in the central nervous system, in the gray and white matter [1], contrary to neurogenic areas that are mostly situated in the hippocampus, the olfactory bulb and all along the lateral ventricule [4]. OPCs represent the most important cycling population in the adult human normal brain [1]. In the mouse, recent studies have revealed that their density is strikingly uniform throughout the brain and that this is achieved through self-repulsion mediated by contact-inhibition [5]. OPCs constantly survey their surroundings by actively extending and retracting filopodia. Growing filopodia retract when they contact another filopodium, causing the cells to be equally spaced from one another and to occupy non-overlapping domains. Their density is also remarkably constant in time: the cells achieve the homeostasis of the normal brain tissue by balancing the differentiation and death process by proliferation, which allows them to keep their density constant on average. Through this dynamical interplay between proliferation, differentiation and death, the precursors can be very rapidly mobilised after a brain injury that has killed mature oligodendrocytes. The surrounding OPCs proliferate and migrate to the lesion area, where they differentiate and replenish the injured area by new mature oligodendrocytes. This is indeed what happens in the case of an acute (such as in traumatic or vascular accidents) demyelinisation episode. The case of multiple sclerosis is more complex, since an active remyelinisation can be performed in early lesions, but is completely inhibited in chronic ones. This disease is thus characterized by the existence of large demyelinated areas that do not heal, leading to serious disability for the patients [6]. The oligodendrocyte precursor cells are suspected to be at the origin of some gliomas. Gliomas account for almost 80% of primary malignant brain tumors, and they result in more years of life lost than any other tumor [7]. Standard treatment only confers a modest improvement in progression and overall survival [8], underscoring the pressing need for the development of new therapies. Gliomas are classified with respect to their malignancy from grade I to IV by the World Health Organization [9]. From grade II, gliomas are invasive, and isolated glioma cells can be found beyond the MRI limits of signal abnormalities, including T2-weighted and FLAIR sequences [10,11]. This characteristic explains the systematic recurrences that are observed after oncological treatments. Diffuse low-grade gliomas, i. e. grade II gliomas, are primary brain tumours that affect young adults, contrary to high-grade glioma that is more frequent in elderly people. At the beginning of their evolution, the low-grade gliomas grow slowly and continuously without angiogenesis [12]. During this phase, we have shown that the cell density in the tumor center is higher than in the normal tissue, but stays moderate. The number of proliferating cells in the tissue also remains low [11]. In [13] a diffusion-proliferation model has been used to estimate the date of birth of low-grade gliomas, and a side result was an estimate of the proliferation coefficient around 1 yr−1. During this slow-growing phase, the mean radius of these tumours measured on MR images increases linearly with time [14,15]. After several years of growth, and despite treatments (surgery, chemotherapy and radiotherapy), angiogenesis is suddenly triggered and low-grade gliomas inexorably evolve into more aggressive forms (glioblastoma multiforme), impeding the social and professional life of the patients [16]. Clearly, new therapeutic strategies are urgently needed. Identifying the cell of origin of gliomas (that could be specific to each type of glioma) could constitute one of these new approaches, by inducing the development of targeted drugs against a specific type of cells. Targeted drugs could be expected to be more efficient against the tumor cells, but to have less side effects on normal tissue than non-specific drugs. Many gliomas express markers characteristic of oligodendrocyte progenitors and the fact that OPCs remain widespread in the adult brain and still actively proliferate in the adult central nervous system make them good candidates for the accumulation of mutations and thus gliomagenesis [17–19]. The high ability of these progenitors to regulate their proliferation according to exogeneous signals may play a role in their susceptibility to transformation. Indeed, in vitro, it is possible to break the homeostasis by adding a growth factor (PDGF) and trigger uncontrolled proliferation: in vivo, in rodents, overexpression of PDGF is often used to induce gliomas [20]. This hypothesis is supported by several experiments with murine animals that show that the OPCs can become highly proliferative and form a malignant tumour after reactivation by inducing mutations [18,21,22]. In the domain of cancer research, mathematical models have been developed for decades to help address the essential questions on tumor growth (for a review, see [23]). At the genetic scale, in the pre-cancerous step, some models describe the dynamics of mutation accumulation that lead to cancer [24,25]. At the cellular scale, other models describe how mutations (or phenotypic changes) can disrupt the homeostasis of hierarchically organized tissues [26–29]. At the tumor level, yet other models focus on the tumor growth dynamics after the appearance of the first tumor cells, their interaction with the micro-environment and the influence of treatments, with discrete [30,31], continuous [32–34] and hybrid approaches [35–37]. The model we introduce here belongs to the class of models at the cellular level, and it addresses the problem of the cell-of-origin. We focus on a very specific tumor (gliomas) and a specific type of precursor cells (OPCs). The stem cells and progenitors density in hierarchically organized tissues such as hematopoietic system, epithelium or colonic crypts is regulated by a more differentiated downstream cell population [38], but this is not the case for OPCs that regulate their density by themselves, through contact inhibition and a density-dependent lifetime. To our knowledge, there is no mathematical model of the dynamics of normal OPCs and the disruption of their homeostasis that leads to gliomas. We therefore designed a specific and original model of the dynamics of these progenitor cells, based on a cellular automaton including the main processes that the OPCs undergo, i. e. self-repulsion, proliferation, migration, differentiation/death. We show that our model is able to reproduce the homeostasis and the repair process that characterize OPCs. In a second part of the paper, we describe the possible genesis of a glioma, by modeling the appearance of an over-proliferating cell, a non-proliferating cell or a cell that has lost its contact inhibition. When a tumor develops, we compare its characteristics to histological samples of human low-grade gliomas.
  Our automaton model consists in a collection of cells which evolve according to a given set of rules. Each cell is modelled by a sphere. The center of the sphere represents the cell soma (modelled as a point), and the sphere models all the filopodia and extensions around it, see Fig 1 (a). The cell radius is the same for all cells in a given simulation and its value, typically R = 50 μm, stays constant during the simulation. The three-dimensional space is continuous, i. e. the position of a cell in space is a triplet of real numbers, the coordinates of the center of the sphere. In order to avoid border issues, we choose periodic boundary conditions. At each iteration, the cells are updated, one after the other, in a random order in order to avoid undesirable correlations. At each iteration, each cell can participate in three different processes: proliferation, migration, differentiation/death (described in what follows). To be able to track individual cells easily, a number is attributed to each cell and stays the same for the whole life of the cell. The data on human biopsies we are using have already been published elsewhere, and we show in this article a new presentation of the data. The clinical data have been presented in details in [11]. Briefly, we searched the database of a previously published study performed in the Sainte-Anne Hospital in Paris (France) and focused on patients in whom a low-grade glioma was newly diagnosed by MRI-based serial stereotactic biopsies according to the Talairach stereotactic method between January 1992 and December 2001. These patients gave their informed consent for storage of the surgical samples for further analyses. We retrospectively selected nine cases of untreated adult patients. We analyzed 44 biopsy samples. After surgery, the biopsy samples were fixed in formalin-zinc and individually embedded in paraffin. Serial sections were cut at 6 μm. Sections used for immunohistochemistry were microwaved in citrate buffer (pH = 6) for antigen retrieval (Micro MedMicro MEDT, Hacker instruments, Winnsboro, SC) at 98deg C for 30 minutes. Ki-67 (MIB-1) immunostaining revealed MIB-1 positive cells (i. e. , cycling cells).
  When starting with a small number of cells, after a transient state characterized by the appearance of damped oscillations, the cell number stabilises, see Fig 2 (b). As in experiments, in our model the cell density reaches a constant value. The equilibrium is a dynamical one since there is still proliferation and differentiation. When a cell disappears, a neighboring cell will move and proliferate in order to fill the gap and to maintain a constant density. The cell density at equilibrium depends on the value of the lifetime threshold and on the value of the proliferation coefficient as will be explained later with a simplified model. The characteristic time of the sharp increase in density at the beginning of each simulation (due to proliferation) is 1/λ. As the density increases, the overlapping between cells increases. This has two consequences: first, the rate of proliferation decreases and, second, many lifetime clocks are triggered. After being triggered, the clocks are increased by a unit at each iteration, and when a clock reaches the lifetime threshold, its cell disappears, see Fig 2 (a). The interplay between proliferation and delayed differentiation leads to oscillations before the cell density reaches a steady state. The period of the oscillations and the coefficient of attenuation depend on the lifetime threshold and on the proliferation coefficient. In order to study these oscillations, we introduce a simplified version of the model, aiming at keeping only the characteristic features that lead to oscillations in the model, but also in order to facilitate the derivation of the mean field approximation. In this case, the three dimensional space is discretised into sites. Each cell occupies only one site (the cells are point-wise in this version). At each iteration, each cell can undergo either proliferation or differentiation/death. The proliferation rule is the same as in the full model. Since cells are point-wise, overlaps are not possible anymore. Therefore, the density-dependent differentiation rule becomes: if one or more neighboring sites of the cell are occupied, the lifetime clock is triggered. We can calculate the mean field approximation of this simplified model. We first define w (t) the mean cell density that corresponds to the density of occupied sites in the automaton. We also define c (t) the density of sites full of cells that haven’t triggered their clock yet. There is no cell motion in this simplified model. We show that the temporal evolution of the two cell densities c (t) and w (t) satisfies respectively Eqs (1) and (2), where D is the lifetime threshold, λ is the proliferation coefficient, and τ is the time calibration parameter, see S1 Appendix: d c (t) d t = λ w (t) [ (1 - w (t) ] - τ w (t) c (t) (1) d w (t) d t = λ w (t) [ (1 - w (t) ] - τ w (t - D) c (t - D) (2) The study of the steady state is detailed in S1 Appendix. To summarise, we show that in the approximation D ≫ 1 λ ≫ 1 the steady state is: w ∞ ≈ 1 - 1 λ D - τ λ D 2 (3) and c ∞ ≈ τ D - τ λ D 2 (4) The study of stability is also developed in S1 Appendix. We show that the period T of the damped oscillations is given by the following expression: T = D 1 - 1 λ D (5) When λD ≫ 1, we can further reduce the expression of the period to: T ≈ D. In S1 Fig, the period T either calculated with the analytical expression (plain line), measured in simulations with the simplified cellular automaton with point-wise cells or measured in simulations with the cellular automaton with spheres, is plotted as a function of lifetime threshold D and as a function of the proliferation coefficient λ. In [45,46], the authors report that after a lesion, the cell density in the lesion increases, displays an overshoot and then reaches its equilibrium value, without any sustained oscillations. In our model, as we showed before, if the cells have a fixed value of the lifetime duration, the system that has been brought far form equilibrium by a lesion for exemple, comes back to equilibrium, with several oscillations. One possible simple way to increase the damping coefficient of the oscillations is to broaden the distribution of lifetime thresholds. We thus choose the lifetime threshold of each cell D in a shifted exponential distribution of lifetimes: D = D 0 - D 1 ln (1 - r) where r is a random number uniformly between 0 and 1 and D0 and D1 are two parameters that characterize respectively the minimum lifetime threshold and the width of the distribution. Since there is no underlying spatial network in the sphere model, the spatial length that will serve as a reference is the radius of a cell, which, experimentally [5], has been found to be around 50 μm. The total volume we consider in the simulations is the volume of a cube of 20-30 cell radii. In the automaton, time is measured as a number of iterations. In order to compare the simulation results with the experimental ones, a calibration of the time scale is necessary by defining the duration of one iteration in a physical time unit. We have chosen one iteration of the automaton to represent 1 h. Two parameters have to be fixed in the simulations: the proliferation coefficient and the lifetime threshold parameter. For the choice of these parameters, we have a first constraint on the lifetime threshold. The latter must be larger than the duration of the cell cycle so that the cells have time to proliferate several times before dying (λD ≫ 1). This condition allows to reach the cell density observed in experiments. In order to fine-tune the values of the time-related parameters we review what is experimentally known, in rodents. Both [45] and [46] present results of the order of magnitude of the time necessary to close a lesion. In these articles, toxin-induced focal lesions are caused in rat brain with ethidium bromide. This drug is cytotoxic and kills all nucleated cell types, creating a lesion of around 500 μm diameter. Right after the lesion, the cells at the border proliferate more and migrate. The cell density inside the lesion increases sharply, overshoots the normal density, reaches a maximum after ten days and finally goes back to normality after 15-30 days. As the time to return to equilibrium is close to the value of the mean lifetime threshold (in the case λD ≫ 1), it means that the value of D0 should be close to 20 days. The second temporal order of magnitude can be found in [5], where the cycle duration of the OPCs was measured in vivo in mice at around 20 days. We found that the parameters λ = 0. 05 h−1, D0 = 25 days and D1 = D0/2 are a very good compromise between the different constraints. We will now present the results of the simulations that can be compared with the experimental data. The simulation presented at first is the no-lesion case. In the simulations, we start with one cell so we first have to wait for the system to reach its equilibrium point. In Fig 3 (a), the total cell density in space is represented. The first remark is that the oscillations are highly damped (as expected) when the total cell density converges towards its equilibrium value at the beginning of the simulation for times t < 80 days. At equilibrium, in the total volume of simulation (a cube of 1 mm3) there are around 2000 cells, meaning that the cell density is equal to 2000 cells per mm3, see Fig 3 (a). In experiments, 160 cells were counted in a volume of 0. 06 mm3, so the experimental density is approximately equal to 2700 cell/mm3 [5]. The two values agree well. Also, in the equilibrium state, other variables have the same order of magnitude as in experiments: We also simulated the repopulation after a lesion killed all the cells in an area. At time t = 125 days after the beginning of the simulation, a lesion of 500 μm diameter is induced see Fig 3 (c). In the total cell density (see Fig 3 (b) ), the lesion is materialized by a peak in the dying/differentiating cells (dashed blue curve), followed by another peak in the proliferating curve (plain blue curve) when the cells at the border of the lesion begin to replenish the empty space (see also S1 Video). The cell density inside the area of the initial lesion overshoots and goes back to its equilibrium value with further oscillations 40 days after the induction of the lesion, see Fig 3 (b). This duration is compatible with the experimental estimates. Fig 3 shows the system of cells at the time of the lesion (c) and after 6 days of evolution (d), when the empty space of the lesion has been filled by new cells (appearing in blue on the figure, meaning that their lifetime clock is low), created by the proliferation and the migration of the cells at the border of the lesion. This behaviour mimics very well what happens in the experiments, where the induction of proliferation and migration of cells around a lesion has been also observed [5]. In this section, we explored the possibility that in our simulations, the appearance of a cell with only one different property than the existing cells could lead to the formation of a glioma. We compared the results of the simulations in each cases with the properties of real low-grade gliomas: a) there is a limited increase in the total and in the proliferating cell density in the tumor; b) the radius of the tumor increases linearly, with a low velocity and c) the tumor is invasive with diffuse boundaries. Concerning the first point, we had the opportunity to have access to stereotactic biopsies from different spots of the same human tumor and we showed that in the center of the tumor, the cell density was higher than that of the surrounding normal tissue, but only by a factor close to 2 [11] (up to a factor of three larger in some cases, but on average 34 samples, a factor of 1. 7), see Fig 4 (a). We studied three possible scenarios for the genesis of a glioma: the appearance of either an immortal cell, or a contact-inhibition free or an over proliferating cell, that transmits its abnormal property to its daughter cells. In the simulations, this abnormal original cell corresponds to a given number fixed at the beginning of the simulation (for exemple the 2000th cell). We vary only one parameter at a time. The first scenario that we studied is the acquisition of a lifetime advantage by the abnormal original cell: its lifetime clock threshold is longer than the duration of simulations, so that the cell cannot die during the simulation. All its other properties are maintained the same (proliferation, migration, contact inhibition). All its daughter cells keep the same property. Therefore, normal dying cells are progressively replaced by these immortal cells. A tumor forms and the cell density inside the tumor reaches a new equilibrium at a larger value than in the normal case, see Fig 5 (a). This is in good agreement with Eq (3) that tells us that the cell density at equilibrium increases with the lifetime. However, since the new cells do not die, the density of proliferating cells decreases down to zero, see Fig 5 (b). This feature is not consistent with the histological samples and the experimental observations on low-grade gliomas (see Fig 4 (a) and 4 (b), and in [11]) where the proliferation inside the tumor is low but higher than in the normal tissue. The second scenario that we studied is the bypass of the saturation of the proliferation rule. This rule stipulates that a cell 1 can proliferate only if the distances between the newly created cell and all the cells that overlap cell 1 are larger than the cell radius R. The abnormal original cell and its daughters are thus not limited in their proliferation, so the cell density increases very fast, see Fig 5 (c) and 5 (d). A tumor forms again but in contrast with the previous scenario, no new equilibrium of the cell density inside the tumor is reached. The cell density explodes and the situation is be closer to what happens in high-grade gliomas (where cells at the center of the tumor are closely packed), than to low-grade gliomas. In order to keep a cell density reasonable, as in low-grade gliomas, the saturation-of-proliferation rule seems to be necessary. The last scenario we studied was the appearance of an abnormal cell with a higher proliferation coefficient. This abnormal original cell is associated with a proliferation coefficient five times higher (0. 25 h−1) than the proliferation coefficient of all the other cells around it (0. 05 h−1). All the progeny of this mutated cell has the same high proliferation coefficient. The holes left by normal cells that die are quickly filled by highly-proliferating ones. As the small highly-proliferating core progressively grows, glioma and normal cells compete for space. Fig 6 depicts the temporal evolution of such a glioma (see also S2 Video). Where both cells types cohabit, Fig 7 (a) reveals that normal cell density decreases with time (blue circles and lines), whereas the glioma cell density increases (red circles and lines). After a long time, the normal cells disappear and all the cells bear the mutation. In the tumor, the over proliferating cells reach a new equilibrium: the cell density is constant and higher than before, as predicted by Eq (3) and as it can be observed in Fig 7 (c). Since the third scenario (an over-proliferating cell) corresponds to the better agreement with the cell density and proliferating cell density observed in real low-grade gliomas, we kept only this scenario for further comparisons with clinical properties of gliomas. As real low-grade gliomas, the simulated tumour is invasive: the border of the tumor are not regular and the tumor cell density decreases progressively from the center to the normal tissue, see Figs 6 and 7 (a). Since the normal cells already have migration properties, there is no need in our model to introduce specific migration properties for tumor cells in order to obtain an invasive tumor. Moreover, the evolution of the mean radius of the core in the three cartesian directions is plotted and fitted, see Fig 7 (b). We find that the growth of the core is linear, and the velocity of growth of this core in the simulations is around 1. 1 mm/yr.
  In this article, we present a model of oligodendrocyte precursor cells dynamics. The model is based on a cellular automaton and the rules that we define are close to the biological reality: the virtual cells can divide, disappear as OPCs (either differentiate or die) after a given lifetime that depends on the local cell density, and move in the direction that will minimise their overlaps with other cells. Our model reproduces the main properties of the OPCs that have been described experimentally in a normal tissue: it allows to achieve the same homeostasis with a constant density and the same behavior after a lesion when the surrounding cells proliferate and migrate in order to fill the injured area. Our model also achieves a very good quantitative agreement with experiments in rodents, since the theoretical and experimental cell density, mean velocity, cycle duration, fraction of proliferating cells are very similar. There is a strong suspicion that OPCs and more generally glial progenitor cells are at the origin of gliomas [19,20,22]. As OPCs retain some of their characteristics from development to adulthood, including the ability to proliferate and migrate, they could easily transform to an uncontrolled growth state. We proposed three different scenarios for the genesis of gliomas: the appearance of a deadless cell, a cell that has lost its contact inhibition for proliferation and an over proliferating cell. The deadless cells could also correspond to non-differentiating cells, since in our simulations, we do not separate death from differentiation. Those three scenarios are not the only ones possible for gliomatogenesis, but are among the most probable [47]. we found that these three scenarios lead to the formation of gliomas, but with very different characteristics: Another feature of real gliomas that is very well reproduced by our in silico glioma is the linear growth of the tumor radius. Actually, this property of the tumor could have been predicted: as long as the rules of the automaton include linear proliferation and migration, the tumor is invasive and the model naturally reproduces the linear increase of the radius of the tumor at large time, observed clinically. This linear behaviour of the tumor radius is the same as the constant velocity observed by the propagating front obtained with the reaction-diffusion model used for gliomas [50]. Here, in the scenario of the appearance of an over-proliferating cell, we measured that the velocity of increase of the tumor radius is around 1. 1 mm/yr. , which is close to the mean velocity measured for human low-grade gliomas of 2 mm/yr [14]. However, this quantitative agreement has to be taken with a grain of salt since the model has been calibrated with mouse data. The proliferation coefficient and lifetime clock may not be the same in humans. Our model shows that the dynamics of OPCs is compatible with the fact that they can be at the origin of gliomas. We provide different scenarios that could lead to the formation of low-grade and high-grade gliomas. Challenging this result with clinical data from histology will constitue the objective of our future work.
summary: Gliomas are the most common brain tumors and result in more years of life lost than any other tumor. Standard treatments only confer a limited improvement in overall survival, underscoring the need for new therapies. Finding the type of cells at the origin of these tumors could lead to the development of new drugs, specifically targeted towards these cells. The oligodendrocyte precursor cells are suspected to be these cells of origin, because they continue to proliferate through all the adult life. In this article, we present a model of the dynamics of these cells, first in the normal brain, and then we extrapolate our model to the pathological situation. We study several scenarios where, from the normal situation, a cell appears with one property different from those of the normal cells. We show that the alteration of only one of the properties of these cells in the model can lead to the formation of gliomas with different aggressiveness and very similar to real gliomas, reinforcing the suspicion that the precursor cells are at the origin of gliomas.
keywords: cell motility
  medicine and health sciences
  cell cycle and cell division
  biological cultures
  cell processes
  cancers and neoplasms
  precursor cells
  cell differentiation
  oncology
  neurological tumors
  physiological processes
  developmental biology
  homeostasis
  cell cultures
  research and analysis methods
  invasive tumors
  animal cells
  glioma cells
  cell biology
  neurology
  physiology
  cell migration
  biology and life sciences
  cellular types
  cultured tumor cells
  glioma
title: Modeling the dynamics of oligodendrocyte precursor cells and the genesis of gliomas