BioMed Research International
Volume 2017 (2017), Article ID 8569328, 12 pages
https://doi.org/10.1155/2017/8569328

Research Article

In Vitro/In Silico Study on the Role of Doubling Time Heterogeneity among Primary Glioblastoma Cell Lines

M.-E. Oraiopoulou,1,2 E. Tzamali,² G. Tzedakis,² A. Vakis,1,3 J. Papamatheakis,4,5 and V. Sakkalis²

¹Department of Medicine, University of Crete, Heraklion, Greece
²Computational Bio-Medicine Laboratory, Institute of Computer Science, Foundation for Research and Technology-Hellas, Heraklion, Greece
³Neurosurgery Clinic, University General Hospital of Heraklion, Heraklion, Greece
⁴Gene Expression Laboratory, Institute of Molecular Biology and Biotechnology, Foundation for Research and Technology-Hellas, Heraklion, Greece
⁵Department of Biology, University of Crete, Heraklion, Greece

Correspondence should be addressed to V. Sakkalis; sakkalis@ics.forth.gr

Received 5 May 2017; Accepted 18 September 2017; Published 31 October 2017

Academic Editor: Sara Piccirillo

Copyright © 2017 M.-E. Oraiopoulou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The application of accurate cancer predictive algorithms validated with experimental data is a field concerning both basic researchers and clinicians, especially regarding a highly aggressive form of cancer, such as Glioblastoma. In an aim to enhance prediction accuracy in realistic patient-specific environments, accounting for both inter- and intratumoral heterogeneity, we use patient-derived Glioblastoma cells from different patients. We focus on cell proliferation using in vitro experiments to estimate cell doubling times and sizes for established primary Glioblastoma cell lines. A preclinically driven mathematical model parametrization is accomplished by taking into account the experimental measurements. As a control cell line we use the well-studied U87MG cells. Both in vitro and in silico results presented support that the variance between tumor staging can be attributed to the differential proliferative capacity of the different Glioblastoma cells. More specifically, the intratumoral heterogeneity together with the overall proliferation reflected in both the proliferation rate and the mechanical cell contact inhibition can predict the in vitro evolution of different Glioblastoma cell lines growing under the same conditions. Undoubtedly, additional imaging techniques capable of providing spatial information of tumor cell physiology and microenvironment will enhance our understanding regarding Glioblastoma nature and verify and further improve our predictability.

1. Introduction

Glioblastoma (GB), a grade IV glioma as categorized by the World Health Organization (WHO) [1], is one of the most aggressive brain cancer types [2] with a poor prognosis for the patient [3], despite the rapid advances in technology and novel therapeutics. One of the most characteristic features of GB that limits therapeutic potential is heterogeneity [4]; both different molecular GB subtypes [5, 6] and subclonal cell populations coexist within the same tumor [7–9]. Hence, the importance of individualized GB treatment and understanding of patient-specific GB pathophysiology is evident and research plans towards this aim are of great interest.

The use of the widely scientifically studied common GB cell lines passaged in lab conditions for decades [10] is nowadays questionable with respect to their clinical relevance in therapeutic outcome prediction and to their ability of representing the extensive heterogeneity observed among patients [11]. To this front, a common GB trend is the use of patient-derived GB cells to enable preclinical physiologic estimations and personalize therapeutic strategy. Basic researchers cooperate with clinicians in order to isolate GB cells and promote the establishment of short-term primary GB cell cultures [12–15], which provide additional results back to the patient. Established methods for biological research and early drug discovery utilize cell lines grown on plastic culture flasks. Over the years, the ability of these in vitro systems to provide biologically relevant answers and describe drug effects is limited due to the fact that they are too simplistic and do not include key players of the phenomenon. Hence, researchers seem to mobilize more realistic experimental approaches such as 3-dimensional (3D) cell cultures [16–20] and/or ex/in vivo implantations [14, 21–23] to better imitate cancer in a mechanistic and conditional way. Biological 3D models comprise an important step to describe the early phases of tumor progression before going to the complexity of in vivo systems.

Biological experiments are strongly linked with computational and mathematical (in silico) models. In silico models offer a systematic framework of understanding the underlying biological processes integrating knowledge and information from multiple biological experiments and/or clinical examinations [24]. By predicting the behavior of the system, new targeted experiments can be designed. In that way, the process of mathematical modeling validation is an iterative refinement procedure [25], which terminates when a valid and biologically plausible and concrete description of the system that reproduces the observed cellular behaviors and growth patterns is found. Several mathematical approaches have been proposed to describe the complex, multiscale spatiotemporal tumor evolution. According to their mathematical perspective, these approaches can be classified into continuum and discrete models. Continuous mathematical models are commonly used to describe tumors at tissue level focusing more on the collective, averaged behavior of tumor cells [26–28]. On the other hand, individual-cell-based models using discrete and hybrid discrete-continuous (HDC) mathematics can describe the behavior of each cancer cell individually as it interacts with its microenvironment. Individual-cell-based models are in general more suitable to describe in vitro experiments, animal models, and small-sized tumors [29–34].

In general, such mathematical models attempt to translate tumor physiology hallmarks [35] into computational parameters and the predicted output is subsequently validated using as ground truth either the experimental [36, 37] or the clinical results [38, 39]. As it is well-understood, both cell division and local spreading are responsible for cancer expansion [40, 41] comprising the most important aspects for cancer progress [30, 42]. Doubling time is defined as the average duration of cell growth and division as reflected by the cell cycle “clock” [43]. GB tumors have a remarkable rapid growth that has a critical role regarding the space-occupation and the development of intracranial pressure, usually the main reason of the GB symptomatology [44]. In previous studies, the significance of the proliferative rate has been shown. More specifically, in [45], the proliferation rates of different breast cancer patients are estimated from subsequent Magnetic Resonance (MR) images in conjunction with a simple logistic tumor growth model and show that the proliferation rate estimates could discriminate patient’s survival and response to therapy. In another study [46], the role of experimental and simulated diffusion gradients in 3D tumors affecting nutrient, oxygen, and drug availability within the tumor and subsequently controlling cell proliferative rate is examined. A mathematical model parameterized from monolayer experiments is used to quantify the diffusion barrier in 3D experiments. In the recent study [40], acquisition of physiologic parameters from multicellular tumor spheroids including proliferation and death spatial profiles is used to constrain and parametrize a mathematical agent-based model that addresses several cell growth mechanisms necessary to explain the experimental observations and reductively translates them to tumor progress over time.

This work utilizes primary tumor cells collected from GB patients and subsequently cultivated in vitro as 3D tumor spheroids. As an initial step towards understanding the GB heterogeneity among patients, we focus on proliferation. The aim of this work is first to mathematically study the important components affecting the growth dynamics of tumor spheroids when motility is inhibited, mainly including the inter- and intratumoral heterogeneity with respect to cell proliferation and, second, to parametrize the mathematical model based on experimentally estimated parameter values of primary GB cell lines in order to increase clinical relevance. Doubling times and the average cell sizes of in-house-established primary GB cell lines from three different patients are used. The well-known U87MG GB cell line is also used as control in the experiments. All the biological experiments are performed simultaneously under the same initial and growth conditions. A hybrid, individual, cell-based mathematical model is used to predict the growth curves of the tumor spheroids and parametrized based on the experimental data. Variations in several mathematical model parameters are explored in order to quantify their effect on tumor growth expansion. The simulated results are compared to the experimental data from the relevant 3D cell cultures and show that, in combination with the proliferation rate, additional factors like the mechanical cell contact inhibition are necessary to predict the in vitro evolution of the different GB cell lines under study.

2. Methods

2.1. Sampling Procedure

Brain tissue sample is collected from the lesions during biopsy or gross resection of patients with indications of GB based on symptoms and MR images, while still naïve from treatment and later histologically proved to be GB cases. For the purposes of this study, we used the primary cells of three different patients. The first is a 70-year-old male patient with de novo GB close to the left brain motor area, also called GBP03 cells. The second, called GBP06 cell line, was collected from a 47-years-old female patient with a tumor in the medulla proven to be a secondary GB, which was gradually evolved to grade IV from lower grades within a time period of approximately 20 years. The third sample, called GBP08, was provided by a 53-year-old male patient with also primary GB in the temporal-occipital left hemisphere. All samples are anonymously provided with the informed patients’ consent by the Neurosurgical Clinic of the General University Hospital of Heraklion, Crete, Greece, while the protocol has been approved by the Institutional Ethical Committees. Because of the relatively low success rate of the primary cell culture establishment, we are limited to these three GB cases for this work.

2.2. Primary Cell Cultures

Later to tissue sampling in saline solution, the specimens are immediately transferred to the lab where they are mechanically dissociated into smaller parts and supplemented culture medium is added (Dulbecco’s modified Eagle medium (DMEM) with 10% fetal bovine serum (FBS) and 1% gentamycin). After gradually removing all cell debris and dead tissue parts, cancer GB cells are cultured as monolayers in standard lab conditions.

As explained before, there is much heterogeneity between GB cases and the protocol of tissue handling is slightly modified per case. An ectopic, subcutaneous implantation to immunodeficient mice is a procedural step that takes place whether the conditional stability cannot be preserved in vitro so that it cannot be assured that the isolated GB cells will survive and proliferate in flask. Therefore, lab animals serve as “living incubators” and usually, after the first implantation, the cells are collected and recultured until the cell culture is successfully established. In this work, GBP03 cells are passaged once, while GBP06 and GBP08 cells are directly used. All possible steps are taken to avoid animal suffering at each stage of the experiments.

2.3. Doubling Time Assay

We use the GBP03, GBP06, and GBP08 primary cell lines as well as the U87MG cells (ATCC® HTB-14™, USA) as control line. In order to measure the doubling time intervals of the different cell types used we apply a simple protocol in adherent cultures. In a 24-well plate, 20000 cells/ml of supplemented DMEM are seeded per cell type at day zero. The plate is incubated in standard lab conditions for approximately a week. Whenever needed, cell culture medium is carefully renewed avoiding the adherent (active) cell population to be disturbed.

Every 24 hours after seeding, the culture medium of one well per cell type is removed and trypsin-EDTA (Sigma-Aldrich, Germany) 1x solution is added for 1-2 minutes. After another 1 minute of trituration in order to produce a single cell solution, all the context is removed from the well and is transferred to a 2 ml Eppendorf tube. As a final step, 4% formaldehyde is added to permanently fix the cells within the tube which is stored to the refrigerator for further use. The procedure is repeated up to the point that 100% cell confluence is achieved. The cell concentration for each cell type is measured with a 24-hour interval by using a hemocytometer.

2.4. Cell Size Estimation

A divided Petri dish is plated with a single cell solution of ~2000 cells/ml and is incubated in standard lab conditions overnight to let the cells adhere in the surface of the dish. Accordingly, brightfield images of attached single cells are captured in 40x magnification and known acquisition parameters to an inverted light microscope (Leica, Germany). To check size and shape homogeneity between each cell population so that to assure that the estimated average cell size will be representative, we capture a photograph of a single cell solution within the fixed grid dimensions of the hemocytometer.

2.5. D Spheroid Generation

We use the hanging-drop technique in order to produce spheroids from each cell type, as recommended in [16, 17, 47]. A single cell solution of 625 cells/50 ul of supplemented double-filtered DMEM is initially seeded per well in a 96-well hanging-drop plate (3D Biomatrix, USA). Two rows of wells per cell type are plated so that approximately 24 spheroids are produced. Agarose solution of 1% w/v is added to plate’s reservoirs to prevent evaporation of the droplets. After 2–4 days of cells aggregating at the bottom of each droplet, we can consider that the spheroids are finally formed. The growth progress of the spheroids is monitored over time via photographs taken under set acquisition parameters to an inverted light microscope (Leica, Germany) for predecided critical time points (2-day interval).

2.6. Data Analysis

The average doubling time of each cell line is estimated using exponential linear regression on the doubling time data. The average cell size of each cell line is estimated by segmenting the area of approximately 10 randomly selected cells in brightfield images to ImageJ [48] and averaging. The tumor expansion of the 3D spheroids is again estimated based on the area shown in their brightfield images. The growth curve is estimated by the mean area value ± standard deviation over time. All the above measurements are evaluated per cell type and many experiments are performed for each cell type.

2.7. Computational Model Implementation of Tumor Spheroids

A simplistic HDC mathematical model is used to describe the observed tumor growth of the 3D in vitro experiments. In the context of the HDC model, each individual cell is described by a discrete cellular automaton, while the local microenvironment is approximated by partial differential equations (PDE). In the following, a concise description of the HDC model is provided, while more thorough description can be found in [49].

2.7.1. Computational Domain

To simulate a central slice of the 3D in vitro tumor spheroids, we set up a 2D regular lattice of size = 5 mm. We assume that each square lattice site fits a single cell; thus the lattice site defines the cell size as well. The same lattice is used by both the discrete and the continuous compartments.

2.7.2. Continuous Compartment

For simplicity, we assume that oxygen is the only limiting molecule required by the cells in order to proliferate. The spatiotemporal evolution of oxygen is described by the partial differential equation (PDE) shown in (1). Oxygen is assumed to diffuse through the domain with diffusion coefficient , decays naturally at a rate , and is consumed by the tumor cells at a rate . The term is 1 if there is a tumor cell at the location or 0 otherwise.

2.7.3. Discrete Compartment

Each tumor cell is an individual entity with its own traits. Sets of these traits are assumed to represent a cellular phenotype. A more detailed description of the cell life cycle can be found in [49, 50].

In this work, two mechanisms of tumor cells are mainly considered: proliferation and death. Cellular movement has been neglected considering that the protocol of the in vitro experiments does not conditionally allow cell motility. Cells die if the local oxygen concentration drops below a defined threshold . When a cell dies, its location is immediately treated as empty space. On the other hand, the live cells incrementally prepare for proliferation at every time step, until the cell age reaches their doubling time. At that moment, the cell searches for a nearby empty space at the 1-Moore neighborhood. If no empty space is available, the search is expanded to the 2-Moore neighborhood (see Figure 1) and the process is repeated up to -Moore neighborhood, where is defined as the proliferation depth and determines the maximum neighborhood size. Examples of Moore neighborhood can be seen in Figure 1. If more than one empty space is found in the same neighborhood, one of them is randomly chosen.

(a)

(b)

(c)

Figure 1: Example of a cell (shown in black) attempting to proliferate. Firstly, the cell searches the 1-Moore neighborhood highlighted by the gray squares in (a). Being unable to find an empty space, it searches the 2-Moore neighborhood indicated by the gray squares in (b) and (c). As an empty space is found, the orange cell is pushed towards the empty space as shown in (b). The latter movement frees the empty space on the 1-Moore neighborhood and allows the proliferating cell to place an identical cell (also shown in black) to the adjacent empty space (c).

As shown in Figure 1, when an empty space is found on a neighborhood other than the 1-Moore, cells are pushed away from the location of the proliferating cell towards the empty space in order to create an empty space to the 1-Moore neighborhood. Then the cell resets its cell age and places a copy of itself at the adjacent empty space. If no empty space has been found, the cell enters a quiescent state at which it constantly searches for empty space, without further increasing its age. The extended proliferating rim describes the maximum distance over which a cell is capable of pushing other cells away in order to create space for its proliferation and reflects the mechanical growth inhibition processes observed in growing cell populations [40].

3. Results

In this work, the in vitro estimated doubling times and cell sizes of three in-house-established primary GB cell lines, as long as of the U87MG cells, are used to initialize the individual-cell-based mathematical model in an attempt to predict their different growth patterns. A sensitivity study is performed where the effects of important factors affecting tumor spheroid expansion such as the doubling time, the cell size, and the depth of the proliferative rim and the coexistence of multiple clones with different proliferative capacities within the tumor are computationally explored. We argue that, as expected, proliferation is one of the most defining characteristics regarding tumor expansion and that tumor predictive computational models should prioritize these remarkable variances between individuals and not just based on theoretically defined values.

3.1. In Vitro 2D Cultures

3.1.1. Cell Size Estimation

A usual answer of what a common human (cancer) cell diameter could be is about 10 to 100 microns [51, 52], and actually, most computational approaches assume cell size within 10–30 microns [29]. In 2D cultures of low confluence, the cell size and shape are in resting state and not crucially influenced by neighboring cells. As depicted in Figure 2, there is much homogeneity in U87MG culture with the cells conforming a rather prolonged typically observed shape, with a soma cell size varying between 19 and 24 microns in diameter (see also Table 1). On the contrast, all primary cells used in our study are smaller and typically round with not many cellular protrusions compared to U87MG cells, yet cells of the same cell line appear to differ within the same population. In case of U87MG cells, it is expected that after all these years in lab conditions there is not much morphological diversity within the cell population and that the cell soma size adequately represents the cell line. On the other hand, regarding primary cells, the cell size is only an average of all possible phenotypes within each cell line. More specifically as denoted in Table 1, GBP03 cells have an average cell diameter of 19 microns, while GBP06 are approximately 16 microns and GBP08 are close to 15 microns in diameter. Also, U87MG cells, when growing in adherent cultures, intrinsically form aggregates when much confluent. On the contrary, the primary cells studied here seem to continue as monolayers no matter the level of confluence. Obviously, the average cell size of a certain cell population, no matter how well represented in 2D, it is not maintained when growing in 3D culturing since other physiological parameters that will be discussed next also affect the cell surface-to-volume ratio altering both size and shape.

Table 1: Mean cell sizes and doubling times (±standard deviation) as estimated from the in vitro experiments for the respective cell lines (first column). The in silico values used to initialize the HDC model regarding the doubling time are also shown.

Cell type
In vitro estimations
In silico values

Cell diameter (μm)
Doubling time (h)
Doubling time (h)

U87MG
21.5
30.8 ± 2.5
33

GBP03
19
25.4 ± 0.5
25

GBP06
16
23.5 ± 0.7
23

GBP08
15
23.0 ± 1.5
22

(a)

(b)

(c)

Figure 2: U87MG cells along with primary GB cells growing as monolayers ((a) 40x magnification) and as hanging-drop spheroids (initial day in (b) and final day in (c), 4x magnification). Scale bars are 50 and 100 microns, respectively. The initial day is set to be the first day of cell aggregation in spheroidal shape after seeding, meaning Days 2–4. Accordingly, the final day is the time point where spheroids start to deform and decompose, usually approaching well’s borders. This day is Day 14 for most primary spheroids.

3.1.2. Doubling Time Estimation

Based on literature, glioma cells usual doubling time ranges from 24 h to a couple of days [53], but more often established primary GB cell lines are recorded to vary few days [12, 54, 55]. Particularly for U87MG cells, they are supposed to have a population doubling time approximating 34 hours, according to their product sheet (ATCC HTB-14, USA). Our measurements presented in Table 1 are in line with the bibliographic records. Specifically, U87MG cells have a mean doubling time of 30.8 ± 2.5 h, which is the slowest division between the cell types we use. Among the primary cell lines, GBP03 cells divide approximately every 25.4 ± 0.5 h, while GBP06 and GBP08 have similar doubling times estimated at 23.5 ± 0.7 h and 23.0 ± 1.5 h, respectively.

3.2. In Vitro GB Spheroids

The hanging-drop technique used here to generate the 3D spheroids is a method conditionally approaching the real avascular tumoral state in vivo [17]. The spheroid size is determined with optical microscopy and monitored over time. It should be noted that, the imaging approach used here cannot give any quantitative estimate of the compactness of the cells or any other spatial information including the number of the cells, the cell size, shape, and polarity, which are definitely different between 2D and 3D structures.

In general, we observe that both primary and U87MG cells need approximately 4 days from single cell solutions to aggregate into spheroidal structures, while during this starting period, they seem to suppress proliferation capacity. However, most often, primary cells aggregate sooner than U87MG ones after seeding.

Figure 2 illustrates the growth area of the in vitro spheroidal domains as imaged in 2D brightfield images at the initial and final day. The growth curves of each cell line are shown in Figure 3. An apparent difference between patients, but also between primary and conventional cell lines, can be observed. To be more specific, all primary spheroids grow larger than the U87MG cells. GBP06 and GBP08 primary spheroids follow an initial fast growing, exponential phase that slows down after approximately 6 days. U87MG spheroids have an almost linear growth pattern. It has to be clarified that the spheroids reach the well’s borders before the plateau and decay phases are observed. The patients GBP06 and GBP08 adopt a high growth pattern, while the patient GBP03 follows an intermediate growth rate closer to the U87MG cell line. As already mentioned, especially for the primary cell lines, the initial distribution of the subclones, when plating the cells (Day zero), is random. This eventually leads to a multifactorial subclonal spheroid growth integrated to average estimations.

Figure 3: Growth of the tumor spheroid area over time for the in vitro experiments of each cell line.

3.3. Computational Parameter Study

Prior to parametrizing and predicting the growth pattern of the multicellular spheroids, a simple parameter study is performed to determine the extent at which the doubling time and cell size affect the 3D growth simulation, as well as explore the effect of additional parameters that could play a significant role in tumor expansion including the depth of the proliferative rim and intratumoral heterogeneity.

The discrete and the continuous part of the computational model are parametrized accordingly to meet the experimental setup as shown in Table 2. The length of the computational domain equals 5 mm to resemble approximately the size of the hanging-drop plate. Both the oxygen decay rate and the cell’s oxygen consumption rate were adopted from [29]. To numerically solve the PDE (1), its parameters have been nondimensionalized by using , , and , which correspond to the maximum oxygen concentration, the computational iteration time, and the domain length, respectively. Dirichlet boundary conditions are used to lock the boundaries to the maximum oxygen concentration to simulate the so-assumed adequate and stable nutrients’ availability, since the culture medium during the experiment is periodically refreshed. Also, the alternating directions implicit method is used to numerically solve the PDE [56, 57].

Table 2: The computational parameters used to initialize the HDC model.

Parameter
Value

Domain length,
5 mm (methods-computational domain)

Cell (& lattice) size,
14–20 μm (methods-computational domain)

Iteration time,
8 h (methods-computational domain [49])

Oxygen consumption,
6.25 10−17 M cell⁻¹ s⁻¹ (methods-computational domain [29])

Maximum oxygen,
6.7 10⁻⁶ M O2 cm⁻³ (methods-continuous compartment [29])

Oxygen decay rate,
0.0125 (ND) (methods-continuous compartment [29, 30, 58])

At first, we explore the effect of doubling time on tumor expansion keeping the rest modeling parameters constant. Specifically, we assume a tumor cell of size equal to 18 μm and consider a depth of proliferative rim equal to 2 cells, while varying the doubling time from 15.5 h to 35.5 h. Figure 4(a) shows the growth curves of the tumors with different doubling times. As expected, increased proliferative capacity results in increased tumor expansion. If a reference time point is picked at 10 days, we can calculate the absolute increase of area yielded by the decrease of the doubling time. When the doubling time is reduced from 35.5 h to 30.5 h, the area increases by approximately 24.46%; while comparing the respective areas between the doubling times 20.5 h and 15.5 h, the area is increased by 54.87%. We can thus conclude that the expansion area is affected more, when the doubling times are lower. As expected, the effect is accumulative; thus if a later/earlier time point was picked the differences would increase/decrease, respectively.

(a)

(b)

(c)

Figure 4: Growth of the tumor spheroid area over time as predicted from the computational model related to altering doubling time from 15.5 h to 35.5 h (a), the cell size from 14 to 20 microns (b), and the proliferation depth from 1 to 5 (c).

We also explore the effect of cell size on the observable tumor expansion. It should be noted that if counting of the tumor cell population was possible on the in vitro experiments, then this parameter would make no difference. We vary the cell size from 14 to 20 μm, while keeping the doubling time constant and equal to 25.5 h and the proliferation depth equal to 2 cells. Figure 4(b) shows that, by increasing the cell size, the tumor expansion increases as well, as expected. Indicatively, by comparing the values at simulation time 10 days, the area relatively increases by 21.5%, 29.8%, and 31.1% as the cell size increases from 14 μm, 16 μm, and 18 μm to 16 μm, 18 μm, and 20 μm, respectively.

The depth of the proliferative rim significantly affects the tumor expansion as it increases the number of proliferative cells. Figure 4(c) illustrates the effect that different proliferation depths have on the tumor area over time. The proliferation time was set to 25.5 h and the cell size to 18 μm. At the reference point of 10 days, as the proliferation depth increases from 1 to 5 cells with a step of 1 cell, the area increases relatively to its previous value by 94.7%, 58.4%, 38.9%, and 31.3%. In other words, a considerable higher expansion of the tumor area (94.7%) is observed when the proliferation depth is increased from 1 to 2, as compared to a change from depth 4 to 5. As the proliferation depth increases, less cells enter the quiescent state and proliferate instead; this is why the growth area is increased.

To further investigate the role of heterogeneity between our cases, we proceed by performing simulations which contain multiple phenotypes identical in all traits except for their respective doubling time. All phenotypes have their cell size set to 18 μm and proliferation depth () equal to 2 cells. The proliferation time is randomly selected for each phenotype at the beginning of the simulation from a uniform distribution in the interval (15.5, 35.5) hours. As shown in Figure 5, to illustrate the impact of the phenotypic multitude, two scenarios are considered inspired by [29]: one at which the number of phenotypes is 100 (shown in green line) and another where 10 phenotypes are randomly selected (shown in purple line). Additionally, given the randomness of the phenotypic initialization, each experimental scenario is repeated 50 times. Figure 5 also shows the area expansion over time for three monoclonal examples with doubling times 15.5 h (red dashed line), 25.5 h (blue dashed line), and 35.5 h (yellow dashed line). Figure 6 illustrates the doubling time of the populations that survive over time. As it can be seen, the mean minimum and the mean maximum values of the doubling time are constant for a long period of time indicating the presence of both the fastest and the slowest populations within the tumor, yet the frequency of these populations becomes progressively unequal with the fastest population to actually overpopulate within the tumor. Thus, a decline to minimum values of the mean doubling time is observed.

Figure 5: Monoclonal and polyclonal tumor area expansion. For the polyclonal case two scenarios are considered: one at which the number of phenotypes is 100 (green line) and another where 10 phenotypes are randomly selected (purple line). Each experiment is repeated 50 times and the corresponding standard deviation is also shown. The mean area of three monoclonal examples with doubling times 15.5 h (red dashed line), 25.5 h (blue dashed line), and 35.5 h (yellow dashed line) is also illustrated.

Figure 6: Doubling time of the populations that survive over time in a polyclonal tumor. Two scenarios are considered: one at which the number of phenotypes is 100 and another where 10 phenotypes are randomly selected. Each experiment is repeated 50 times. The minimum, maximum, and average doubling times for both scenarios are shown, as well as their corresponding standard deviations.

3.4. Comparison of Biological and Computational Results

In the following, we assume monoclonal populations and parametrize the mathematical model based on the estimated experimental values for the doubling time and cell size for the different GB cell lines. We also parameterize the model without taking into account the in vitro estimates of cell sizes and keep the cell size and all the other parameters constant in all the experiments. Parameters within the range of the experimental biological observations are chosen to achieve the best-fitting growth curves. It has to be noted that both the simulated and the biological experiments have an initial seeding population of approximately 625 cells per spheroid per cell type. The simulations show that the in vitro estimates of cell sizes do not improve the model predictability and that accounting only for differences in doubling time among GB lines results in very similar growth curves.

Table 1 shows the parameters used by the in silico model regarding the doubling time. Figure 7 shows the in vitro growth curves and the in silico predicted ones for all the GB cell lines. Based on the selected doubling time values and keeping the proliferation depth equal to 2, the growth curves of U87MG and GBP03 cell lines are closely approximated by the in silico model. However, the GBP06 and GBP08 cell lines diverge significantly from the in vitro results indicating that proliferation alone is necessary, but not sufficient to explain the tumor expansion of different GB cell lines growing under the same initial conditions. Hence, additional phenomena should be taken into account. For example, increasing the proliferative depth and/or consider the possibility that multiple phenotypes with various proliferative capacities coexist within such tumors, then the in vitro and in silico growth curves would come in line as our parameter study analysis previously revealed. Alternatively one could advocate that GBP06 and GBP08 contain phenotypes with higher proliferation depth than U87MG (and GBP03) which are expected to thrive in compact environments such as a solid spheroid. It should be noted that the proliferative depth could also be affected by the development of extracellular matrix (ECM) substrate in 3D cultures, even in the conditional absence of a relevant substrate [17], as in our biological experiments. This, along with antagonistic and synergetic relationships of subclones within the growing spheroid, could alter the mechanical responses of dividing cells, reflected in terms of proliferation depth to our mathematical model. However, our biological approach did not take into account a priori this parameter, but it was the computational approach that indicates such possible behavior suggesting that ECM production and distribution might also be different in different cell lines.

Figure 7: In vitro spheroidal growth as opposed to in silico for all four cell types with the final chosen sets of doubling times as shown in Table 1 and fixed proliferation depth equal to 2. Two additional simulated growth curves are depicted with different proliferation depth values for the GBP06 (, yellow dash dotted line) and the GBP08 (, purple dash dotted line) spheroids.

Figure 7 also shows the simulated growth curves for the GBP06 and GBP08 after changing their proliferation depth values from 2 to 4 and 3, respectively. The in vitro data better correlate the relevant in silico data. Also notice that setting the proliferation depth of GBP06 higher than the GBP08 is important to achieve their corresponding growth patterns, where GBP06 grows faster than GBP08, given that the doubling time of the former is higher than the latter and that small differences in their cell sizes are not adequate to reverse their growth patterns. Another point that should be marked is that the subsequent decline observed after Day 8 in the in vitro growth curves of these two cell types cannot be predicted by the computational model. This is because the computational model we use does not account for inhibitory stimuli that are probably developed in real growing tumors, since this was beyond the scope of this study.

4. Discussion

This work utilizes primary tumor cells collected from GB patients and subsequently cultivated in vitro as 3D tumor spheroids and computational approaches to study, experimentally parametrize, and predict the growth dynamics of tumor spheroids focusing on proliferation. At first, a parameter study is performed in order to evaluate the extent to which important factors such as the doubling time, the cell size, and the depth of the proliferative rim, as well as the coexistence of multiple clones with different proliferative capacities within the tumor, affect tumor spheroid expansion when motility is inhibited. The experimentally estimated doubling times and cell sizes of three in-house-established primary GB cell lines, as long as of the U87MG cells, are then used to parametrize the computational individual-cell-based model.

Overall the parameter study verifies the significant effect of proliferation (depicted in both the cellular doubling time and the depth of the proliferative rim) on tumor expansion [40] and underlines additional factors that could play an important role on tumor growth curves including the intratumoral heterogeneity that has been widely observed in GB. We also observe that a multiclonal population with the same mean proliferation exhibits a greater tumor expansion than the corresponding monoclonal population because fitter clones survive over time driving tumor expansion at higher rates. Furthermore, the clonal heterogeneity within the tumor mass allows different clones to be selected every time an experiment is performed. Thus, a variation is observed in the growth curves. The variance is cumulative, increases over time, and can reach a difference of 100 μm in radius after 14 days of growth (Figure 5). Furthermore, the simulations also show that although the mean growth curves are quite similar, the variance highly depends on the initial number of different clones coexisting within the tumor mass such that fewer initial clones in the population produce higher variability (Figure 6).

Comparing the in vitro experiments with the in silico predictions, we observe that although the proliferation rate is necessary, yet it is not sufficient, to describe the growth curves we observe experimentally. The simulations show that additional factors including the intratumoral heterogeneity together with the overall proliferative capacity reflected in both the proliferation rate and the mechanical cell contact inhibition can predict the evolution of different GB cell lines. Nevertheless, further investigation of the underlying mechanisms is critical.

In general, compactness of the spheroids can be assigned to two factors in mesoscopic terms: (a) the cellularity, in means of cells’ size and shape given the space, and (b) the levels of stress tolerance, reflecting their response against internal forces within the spheroid which vary between division and entering quiescence state, also known as “contact inhibition.” As smaller in size and quicker regarding divisions, GBP06 and GBP08 cells appear to grow larger in 3D over time than the other two cell types mainly because of their promoted proliferative capacity reflected by the higher proliferation depth in the respective simulated growth curves (see Figure 7). However, this is only an assumption for our in silico trials since there is no indication of spheroids cell density and proliferation depth to our experimental protocol and this is a limitation of our method needed to be taken into account in future work.

The migratory capability of our cells is conditionally blocked to our experiments so that it can be assumed to play a minor role in the proliferative characteristics studied here. However, when the different cell populations grow in 3D, both ECM can be produced, and the cell shape and polarity could also be affected, such that cell-to-cell and cell-to-matrix adhesion properties could be further explain the divergence observed over time in growth patterns between the in vitro and in silico experiments.

We suggest that, instead of using bibliographic values usually referenced by common GB cell lines, cell doubling time was found to critically enhance the in silico predictability but is insufficient to holistically describe differences in tumor growth over time among the different GB cell lines. The mechanical cell responses to internal forces obtained during the growth of a compact tumor should be further investigated experimentally, as well as the important role of intratumoral heterogeneity. The importance of quantitative methods to provide spatial information of proliferative, quiescent, and necrotic cells as well as additional features including the remodeling of ECM and phenotypic distribution regarding intratumoral heterogeneity affecting tumor expansion becomes evident.

5. Conclusions

The massive proliferation is a defining characteristic of the tumor nature, essential for its progress. When focusing on such a greed form of cancer, such as GB, constantly growing intra-axially and aggressively disturbing brain functionality, proliferation underlying processes become incompatible in cancer progress. In GB, heterogeneity is another typical hallmark, not only among patients with differences between GB molecular subtypes, but more unexpectedly, between different regions of the same tumor with the presence of intratumoral subclonal dormancies. We claim that future research should be based on primary cells directly collected from patients and that common cell lines should only serve as landmarks to unite studies of different groups. For every primary established cell line, not only molecular but also physiological parameters should be estimated to enable a more precise future clustering of different GB cases. Estimations starting with the typical doubling time as shown here and evolving to more delicate features such as delineation of necrotic and hypoxic regions or invasive capability or others are highly important. To this front, computational models may serve as predictor tools not only for estimating cancer progress [59], but also for designing targeted biological experiments. Simulations of cancer progress, either in vitro or in silico, should not anymore be based on theoretical values, especially if clinical translation is of interest. If we target the holistic description of tumor evolution, we should follow a stepwise approach, where computational tools can definitely help in identifying the most important parameters affecting the final outcome.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This study was supported by grants from the European Commission under the 7th Framework Programme Project Computational Horizons in Cancer (CHIC).

References

D. N. Louis, A. Perry, G. Reifenberger et al., “The 2016 World Health Organization Classification of Tumors of the Central Nervous System: a summary,” Acta Neuropathologica, vol. 131, no. 6, pp. 803–820, 2016.

D. Schiffer, Brain tumor pathology: current diagnostic hotspots and pitfalls, Springer, Dordrecht, Netherlands, 2006.

D. Sturm et al., Hotspot Mutations in H3F3A And IDH1 Define Distinct Epigenetic And Biological Subgroups of Glioblastoma, vol. 22, no. 4, pp. 425–437, 2012.

M.-D. Inda, R. Bonavia, and J. Seoane, “Glioblastoma multiforme:A look inside its heterogeneous nature,” Cancers, vol. 6, no. 1, pp. 226–239, 2014.

M. L. Goodenberger and R. B. Jenkins, “Genetics of adult glioma,” Cancer Genetics, vol. 205, no. 12, pp. 613–621, 2012.

F. Ali-Osman, “Brain tumors,” in Contemporary Cancer Research, p. 393, Humana Press, Totowa , NJ, USA, 2005.

A. Sottoriva, I. Spiteri, S. G. M. Piccirillo et al., “Intratumor heterogeneity in human glioblastoma reflects cancer evolutionary dynamics,” Proceedings of the National Acadamy of Sciences of the United States of America, vol. 110, no. 10, pp. 4009–4014, 2013.

Y. Atlasi, L. Looijenga, and R. Fodde, “Cancer Stem Cells, Pluripotency, and Cellular Heterogeneity. A WNTer Perspective.,” Current Topics in Developmental Biology, vol. 107, pp. 373–404, 2014.

J. Chen, R. M. McKay, and L. F. Parada, “Malignant glioma: Lessons from genomics, mouse models, and stem cells,” Cell, vol. 149, no. 1, pp. 36–47, 2012.

M. J. Clark, N. Homer, B. D. O'Connor et al., “U87MG decoded: the genomic sequence of a cytogenetically aberrant human cancer cell line,” PLoS Genetics, vol. 6, no. 1, Article ID e1000832, 2010.

M. Allen, M. Bjerke, H. Edlund, S. Nelander, and B. Westermark, “Origin of the U87MG glioma cell line: Good news and bad news,” Science Translational Medicine, vol. 8, no. 354, Article ID 354re3, 2016.

C. S. Mullins, B. Schneider, F. Stockhammer, M. Krohn, C. F. Classen, and M. Linnebacher, “Establishment and Characterization of Primary Glioblastoma Cell Lines from Fresh and Frozen Material: A Detailed Comparison,” PLoS ONE, vol. 8, no. 8, Article ID e71070, 2013.

L. F. Pavon, L. C. Marti, T. T. Sibov et al., “In vitro analysis of neurospheres derived from glioblastoma primary culture: A novel methodology paradigm,” Frontiers in Neurology, vol. 4, Article ID Article 214, 2014.

B. L. Carlson, J. L. Pokorny, M. A. Schroeder, and J. N. Sarkaria, “Establishment, maintenance and in vitro and in vivo applications of primary human glioblastoma multiforme (GBM) xenograft models for translational biology studies and drug discovery,” in Curr Protoc Pharmacol, vol. 14, Chapter 14, unit 14.16, 2011.

K. M. Joo, J. Kim, J. Jin et al., “Patient-Specific Orthotopic Glioblastoma Xenograft Models Recapitulate the Histopathology and Biology of Human Glioblastomas In Situ,” Cell Reports, vol. 3, no. 1, pp. 260–273, 2013.

R. Foty, “A simple hanging drop cell culture protocol for generation of 3D spheroids,” Journal of Visualized Experiments, no. 51, Article ID e2720, 2011.

White paper, 5 Reasons Cancer Researchers Adopt 3D Cell Culture: A Review of Recent Literature, 3D Biomatrix Inc., Michigan, USA, 2013.

J. Friedrich, C. Seidel, R. Ebner, and L. A. Kunz-Schughart, “Spheroid-based drug screen: considerations and practical approach,” Nature Protocols, vol. 4, no. 3, pp. 309–324, 2009.

M. Zanoni, F. Piccinini, C. Arienti et al., “3D tumor spheroid models for in vitro therapeutic screening: A systematic approach to enhance the biological relevance of data obtained,” Scientific Reports, vol. 6, Article ID 19103, 2016.

M. Vinci, S. Gowan, F. Boxall et al., “Advances in establishment and analysis of three-dimensional tumor spheroid-based functional assays for target validation and drug evaluation,” BMC Biology, vol. 10, article no. 29, 2012.

M. A. Grotzer, A. Neve, and M. Baumgartner, “Dissecting brain tumor growth and metastasis in vitro and ex vivo,” Journal of Cancer Metastasis and Treatment, vol. 2, no. 5, pp. 149–162, 2016.

R. A. Morgan, “Human tumor xenografts: The good, the bad, and the ugly,” Molecular Therapy, vol. 20, no. 5, pp. 882–884, 2012.

T. Voskoglou-Nomikos, J. L. Pater, and L. Seymour, “Clinical predictive value of the in vitro cell line, human xenograft, and mouse allograft preclinical cancer models,” Clinical Cancer Research, vol. 9, no. 11, pp. 4227–4239, 2003.

V. Sakkalis, S. Sfakianakis, E. Tzamali et al., “Web-based workflow planning platform supporting the design and execution of complex multiscale cancer models,” IEEE Journal of Biomedical and Health Informatics, vol. 18, no. 3, pp. 824–831, 2014.

A. Roniotis, K. Marias, V. Sakkalis, G. D. Tsibidis, and M. Zervakis, “A complete mathematical study of a 3D model of heterogeneous and anisotropic glioma evolution,” in Proceedings of the 31st Annual International Conference of the IEEE Engineering in Medicine and Biology Society: Engineering the Future of Biomedicine, EMBC 2009, pp. 2807–2810, usa, September 2009.

A. Gerisch and M. A. Chaplain, “Mathematical modelling of cancer cell invasion of tissue: local and non-local models and the effect of adhesion,” Journal of Theoretical Biology, vol. 250, no. 4, pp. 684–704, 2008.

K. R. Swanson, R. C. Rockne, J. Claridge, M. A. Chaplain, E. C. Alvord Jr., and A. R. A. Anderson, “Quantifying the role of angiogenesis in malignant progression of gliomas: In Silico modeling integrates imaging and histology,” Cancer Research, vol. 71, no. 24, pp. 7366–7375, 2011.

S. M. Wise, J. S. Lowengrub, H. B. Frieboes, and V. Cristini, “Three-dimensional multispecies nonlinear tumor growth—I: Model and numerical method,” Journal of Theoretical Biology, vol. 253, no. 3, pp. 524–543, 2008.

A. R. A. Anderson, “A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion,” Mathematical Medicine and Biology, vol. 22, no. 2, pp. 163–186, 2005.

A. R. Anderson, K. A. Rejniak, P. Gerlee, and V. Quaranta, “Modelling of cancer growth, evolution and invasion: bridging scales and models,” Mathematical Modelling of Natural Phenomena, vol. 2, no. 3, pp. 1–29, 2007.

F. Michor and K. Beal, “Improving Cancer Treatment via Mathematical Modeling: Surmounting the Challenges Is Worth the Effort,” Cell, vol. 163, no. 5, pp. 1059–1063, 2015.

J. P. O’Connor et al., “Imaging biomarker roadmap for cancer studies,” Nature Reviews Clinical Oncology, 2016.

A. L. Baldock, R. C. Rockne, A. D. Boone et al., “From patient-specific mathematical neuro-oncology to precision medicine,” Front Oncol, vol. 3, p. 62, 2013.

J. S. Lowengrub, H. B. Frieboes, F. Jin et al., “Nonlinear modelling of cancer: bridging the gap between cells and tumours,” Nonlinearity, vol. 23, no. 1, pp. R1–R91, 2010.

M. G. McNamara, S. Sahebjam, and W. P. Mason, “Emerging biomarkers in glioblastoma,” Cancers, vol. 5, no. 3, pp. 1103–1119, 2013.

J. C. Anderson, C. W. Duarte, K. Welaya et al., “Kinomic exploration of temozolomide and radiation resistance in Glioblastoma multiforme xenolines,” Radiotherapy & Oncology, vol. 111, no. 3, pp. 468–474, 2014.

C. Athale, Y. Mansury, and T. S. Deisboeck, “Simulating the impact of a molecular 'decision-process' on cellular phenotype and multicellular patterns in brain tumors,” Journal of Theoretical Biology, vol. 233, no. 4, pp. 469–481, 2005.

T. E. Yankeelov, N. Atuegwu, D. Hormuth et al., “Clinically relevant modeling of tumor growth and treatment response,” Science Translational Medicine, vol. 5, no. 187, p. 187ps9, 2013.

M. De Jong, J. Essers, and W. M. Van Weerden, “Imaging preclinical tumour models: Improving translational power,” Nature Reviews Cancer, vol. 14, no. 7, pp. 481–493, 2014.

N. Jagiella, B. Müller, M. Müller, I. E. Vignon-Clementel, and D. Drasdo, “Inferring Growth Control Mechanisms in Growing Multi-cellular Spheroids of NSCLC Cells from Spatial-Temporal Image Data,” PLoS Computational Biology, vol. 12, no. 2, Article ID e1004412, 2016.

B. Hegedüs, A. Czirók, I. Fazekas, T. Bábel, E. Madarász, and T. Vicsek, “Locomotion and proliferation of glioblastoma cells in vitro: Statistical evaluation of videomicroscopic observations,” Journal of Neurosurgery, vol. 92, no. 3, pp. 428–434, 2000.

D. Hanahan and R. A. Weinberg, “Hallmarks of cancer: the next generation,” Cell, vol. 144, no. 5, pp. 646–674, 2011.

A. Bertuzzi, A. Gandolfi, C. Sinisgalli, G. Starace, and P. Ubezio, “Cell loss and the concept of potential doubling time,” Cytometry, vol. 29, no. 1, pp. 34–40, 1997.

P. Black, “Management of malignant glioma: role of surgery in relation to multimodality therapy,” Journal of NeuroVirology, vol. 4, no. 2, pp. 227–236, 1998.

N. C. Atuegwu, L. R. Arlinghaus, X. Li et al., “Parameterizing the logistic model of tumor growth by DW-MRI and DCE-MRI data to predict treatment response and changes in breast cancer cellularity during neoadjuvant chemotherapy,” Translational Oncology, vol. 6, no. 3, pp. 256–264, 2013.

H. B. Frieboes, M. E. Edgerton, J. P. Fruehauf et al., “Prediction of drug response in breast cancer using integrative experimental/computational modeling,” Cancer Research, vol. 69, no. 10, pp. 4484–4492, 2009.

Y.-C. Tung, A. Y. Hsiao, S. G. Allen, Y.-S. Torisawa, M. Ho, and S. Takayama, “High-throughput 3D spheroid culture and drug testing using a 384 hanging drop array,” Analyst, vol. 136, no. 3, pp. 473–478, 2011.

C. A. Schneider, W. S. Rasband, and K. W. Eliceiri, “NIH Image to ImageJ: 25 years of image analysis,” Nature Methods, vol. 9, no. 7, pp. 671–675, 2012.

G. Tzedakis, E. Tzamali, K. Marias, and V. Sakkalis, “The importance of neighborhood scheme selection in agent-based tumor growth modeling,” Cancer Informatics, vol. 14, pp. 67–81, 2015.

G. Tzedakis, E. Liapis, E. Tzamali, G. Zacharakis, and V. Sakkalis, “A hybrid discrete-continuous model of in vitro spheroid tumor growth and drug response,” in Proceedings of the 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBC 2016, pp. 6142–6145, usa, August 2016.

M. M. Melicow, “The three steps to cancer: a new concept of cancerigenesis,” Journal of Theoretical Biology, vol. 94, no. 2, pp. 471–511, 1982.

D. A. Guertin and D. M. Sabatini, Cell Size Control, in in eLS2001, Hoboken, NJ, United States, John Wiley & Sons, in eLS2001, 2001.

A. L. Stensjoen et al., “Growth dynamics of untreated glioblastomas in vivo,” Neuro-Oncology, vol. 17, no. 10, pp. 1402–1411, 2015.

E. Sandén, S. Eberstål, E. Visse, P. Siesjö, and A. Darabi, “A standardized and reproducible protocol for serum-free monolayer culturing of primary paediatric brain tumours to be utilized for therapeutic assays,” Scientific Reports, vol. 5, Article ID 12218, 2015.

C. M. L. Machado, A. Schenka, J. Vassallo et al., “Morphological characterization of a human glioma cell line,” Cancer Cell International, vol. 5, article no. 13, 2005.

J. Douglas, “Alternating direction methods for three space variables,” Numerische Mathematik, vol. 4, pp. 41–63, 1962.

E. L. Wachspress and G. J. Habetler, “An alternating-direction-implicit iteration technique,” Journal of the Society for Industrial and Applied Mathematics, vol. 8, pp. 403–423, 1960.

P. Hinow, P. Gerlee, L. J. McCawley et al., “A spatial model of tumor-host interaction: application of chemotherapy,” Mathematical Biosciences and Engineering, vol. 6, no. 3, pp. 521–546, 2009.

E. Skounakis et al., “DoctorEye: A clinically driven multifunctional platform, for accurate processing of tumors in medical images,” The Open Medical Informatics Journal, vol. 4, no. 3, pp. 105–115, 2010.