Fine-grained simulations of the microenvironment of vascularized tumours

TL;DR

This paper introduces a computational tool for simulating the microenvironment of vascularized tumours, focusing on ecological and evolutionary dynamics that contribute to tumour heterogeneity and drug resistance.

Contribution

The paper presents a novel computational simulation method to study tumour microenvironment ecology and its evolutionary processes.

Findings

Tumour microenvironment acts as an ecological niche promoting heterogeneity.

Evolutionary dynamics within the microenvironment influence drug resistance.

Simulation results reveal how tumour clones adapt to ecological niches.

Abstract

One of many important features of the tumour microenvironment is that it is a place of active Darwinian selection where different tumour clones become adapted to the variety of ecological niches that make up the microenvironment. These evolutionary processes turn the microenvironment into a powerful source of tumour heterogeneity and contribute to the development of drug resistance in cancer. Here, we describe a computational tool to study the ecology of the microenvironment and report results about the ecology of the tumour microenvironment and its evolutionary dynamics.

Full text

Introduction

The tumour microenvironment is characterized by large chemical gradients and contains a mixture of normal and tumour cells. The presence of high gradients favours the formation of different ecological niches which can become an important source of tumour heterogeneity [1, 25, 4, 39, 3, 13].

The actual size of the niches is quite small, as it is mostly determined by the structure of the scaffold of capillary vessels that envelope and feed the tumour mass. This structure is highly irregular, with a chaotic blood vessel network, a missing lymphatic network, elevated acidity, poor oxygenation, and high interstitial fluid pressure. The distance between blood vessels may be as large as a few hundred $\mu$m, and the regions in between them can become highly hypoxic [39].

The fine-graininess and the large spatial variability are essential features in the formation of the ecological niches and they must belong to any mathematical model that aims to explain tumour heterogeneity as the result of Darwinian selection driven by the local environment. The importance of the microenvironment is widely recognized, and several researchers have tackled the problem of its description and understanding from the computational point of view, e.g. [24, 28, 15, 37, 30, 40, 41] , and the reviews [21, 10].

The studies, however, mainly focus on the biochemical, cellular and biophysical properties of the tumour microenvironment and not on its active role in providing varied evolutionary paths to genetic variants of the tumour cells. As it has already been pointed out [14, 12], the adaptive evolution of tumour clones (central concept of Darwinian dynamics) is driven by the formation of new environmental niches. Many practical difficulties limit the experimental study of the adaptation process, while computer simulations can shed light — albeit in a limited way — on the dynamics of many steps like the convergent evolution of different genotypes to the same phenotype, and the selective loss of specific cell functions.

Simulations of avascular solid tumours show that the microenvironment of these small cell aggregates is formed by rather homogeneous niches with smooth gradients of oxygen, of other nutrients, waste molecules and cell viability [18]. After the angiogenic transition, however, the microenvironment differentiates in unpredictable ways. To identify underlying processes and provide reasoning, we have developed a computational model that is a tool to understand the fine-grained features of a simplified tumour microenvironment, i.e. of a microenvironment that contains only tumour cells and blood vessels. This is a rather strong simplifying assumption, since experimental observations show that in addition to tumour cells the microenvironment contains several other important cell types such as stromal and immune cells. We wish to stress that this assumption is a necessary compromise, because no computer model can presently describe the huge biological complexity of the tumour microenvironment and capture its growth at appropriate spatial and temporal resolutions.

Our computational model combines two different simulation programs, namely, a lattice-free simulation of small avascular solid tumours and a lattice-based simulation of the blood vessels dynamics in solid tumours. As we discuss below, both models have been individually validated with experimental data, and here we provide a first validation of the combined model by correctly reproducing the relevant observed gradients of tumour oxygenation and other biochemical and biological features. Next, we use the combined model to explore how the tumour microenvironment takes shape in small solid tumours at the transition from the avascular to the vascular growth phases, the so-called angiogenic switch [5]. We follow the complex dynamics of this process in detail and in real time and show that, already at these early stages of tumour development, the diversity of environmental niches is high enough to promote adaptive evolution of tumour variants.

Brief outline of “Virtual Biology Lab” (VBL) and Tumorcode

In this section we give a short description of the two programs that have been merged. The present account is very brief as both programs have been described in detail in past papers (for VBL, see [18, 7, 8, 20, 26], for Tumorcode, see [22, 11, 33]).

VBL

VBL is a lattice-free, cell-based program that can simulate both: the growth and proliferation of disperse cells and of more complex, avascular cell clusters (tumour spheroids). It contains a detailed simulation of the metabolic activity of each cell including discrete events at the individual cell level (for instance: mitochondrial partitioning at mitosis), and has the potential to activate phenotypic differentiation. More specifically, the model has the following features:

• •

The model of human tumor cells includes both- the internal biochemical processes and a phenomenological description of the biomechanics of cells.

• •

Spatially, each cell is a two-compartment structure, the inside of the cell and the adjacent extracellular environment or intercellular space; this is important to handle simple diffusion and facilitated diffusion across the cell membrane at the same time.

• •

Discrete events are simulated and interleaved between successive deterministic steps; the random nature of some of them contributes to the stochasticity of the simulation as a whole.

• •

Each cell is characterized by its own phenotype which means that a specific set of parameters is linked to an individual cell.

• •

Enzyme activity is modulated both by pO2 concentration and by pH.

• •

Weighting of molecular paths and cellular mechanisms is possible.

• •

The program simulates dispersed tumor cells and cell clusters (tumor spheroids) including the intercellular space.

• •

The surrounding environment is an integral part of the simulation.

The numerical methods used in VBL have been described elsewhere [6], here we only mention that although the simulation is lattice-free, continuum processes are actually discretized on the continually-variable irregular lattice defined by the Delaunay triangulation [9] based on the cells’ centres which also defines the proximity relations among cells that are necessary to compute the cell-cell forces. We find that the computational geometry library CGAL (https://www.cgal.org) [29] is an easy tool that offers necessary triangulation features amended by the possibility to calculate alpha shapes. This is necessary to define the surface or contact zone of the spheroid with the environment. It is important to note that the simulation spans many orders of magnitude in time which makes the large system of differential equations describing the deterministic part of the program extremely stiff. Therefore we solve the equations by means of implicit integration methods, that allow for comparatively large time steps (of the order of 1-10 s in terms of simulated time) [19].

Figure 1 shows an example of a tumour spheroid simulated by VBL. As we have shown in the past, the main metabolic, morphologic and kinetic parameters are correctly reproduced both for isolated cells and for cells grown as three- dimensional clusters [6, 8].

Tumorcode

While cell-cell interactions like in VBL naturally produce nearly spherical structures, modelling the growth of blood vessels moves the complexity of the simulation one step higher to the level of tissue morphogenesis. The list of models is long. For an appropriate discussion of the existing models for tumor vasculature see [23].

Solid tumours grow in originally healthy tissue featuring a normal vasculature build up by a hierarchically organized arterio-venous blood vessel network that is then dynamically modified by the growing tumor. An integrative modeling approach thus has to address two issues: first, it has to find an appropriate representation of the original vasculature of the host tissue; and second, the dynamics of the given blood vessel network has to be defined, which includes the insertion of new vessels via angiogenesis as well as the removal of existing vessels (old and new) via vessel regression, and the modification of existing vessels via dilation, constriction or occlusion. Moreover, the network carries a blood flow and transports and emits oxygen, nutrients and drugs, which has to be represented, too. Various modeling approaches have been presented in the past ( for a review see e.g. [23, 36, 22] ), here we use the simulation framework Tumorcode described in [11], which does not only create artificial arterio-venous blood vessel networks representing the healthy initial vasculature, but includes vascular remodeling during solid tumour growth by angiogenesis, vessel dilation, regression and collapse. The created networks were tested against experimental data and successfully reproduced:

1. morphology and flow characteristics [31, 32],
2. interstitial fluid flow and pressure [35], and
3. tissue oxygenation [33] for both, healthy and tumour tissue.

Again, here is a short list of its main features of Tumorcode:

• •

Hemodynamics includes the phase separation effect (Fåhræus-Lindqvist effect) and different rheologies.

• •

Concentrations in the surrounding environment are computed according to a continuum model, i.e., according to a set of partial differential equations discretized by finite elements method solved with the Trilinos C++ library (https://trilinos.org) [17].

• •

Angiogenesis, i.e., the addition of new segments, is driven by the VGEF gradient and is partly stochastic [22].

• •

As the tumour grows, the vasculature and the blood flows change, and the modified shear stress leads to blood vessel dilation or to blood vessel collapse [22].

• •

A low local VGEF growth factor concentration produces blood vessel regression [22].

• •

The program computes the interaction with surrounding tissues to compute the extraction of oxygen and nutrients from blood vessels.

• •

The program returns blood flow, oxygen concentration, metabolite concentration, etc, in blood vessels and surrounding tissues.

• •

The blood vessel hematocrit acts as source for oxygen dissolved in tissue and the tissue simultaneously consumes the oxygen. In [34] we explain how the coupled set of non linear differential equations is solved to obtain the oxygen partial pressure for each blood vessel segment.

Figure 2 shows an example of a tumour vasculature simulated by Tumorcode. You can clearly distinguish the vasculature inside the tumour from the surrounding healthy tissue. Note the difference in length scale when comparing Figure 2 and Figure 1, and the fact that the tumour is a continuous structure ( governed by reaction-diffusion equations [11]) in Tumorcode. In contrast to VBL, it cannot describe the evolutionary features of the microenvironment at the single cell level.

The merging of the two programs

To merge the two programs into one, we had to match the spatial structures of the two programs. The cells are not confined to a lattice however the diffusion processes in VBL are calculated by means of a disordered Delaunay triangulation which is not permanent and has a variable size, as cells position change at each time step while the cells grow and proliferate. The situation is almost the reverse in Tumorcode where blood vessel segments belong to an FCC lattice and the continous space is discretized by a cubic lattice to solve the partial differential equations. We interpolated cell properties from the cells position to the surrounding nodes of the cubic lattice (compare Figure 3). In particular, the Tumorcode requires an oxygen consumption field and each cell is considered as point- like source of growth factor. Once the growth factor field exceeds the threshold, the vessel remodelling of Tumorcode starts.

We estimate that the accurate calculation of any substance concentration for a given cell requires $\mathcal{O}(N_{\mathrm{C}}N_{\mathrm{BV}})$ number of operations, where $N_{\mathrm{C}}$ is the number of cells and $N_{\mathrm{BV}}$ is the number of blood vessels segments. Since both $N_{\mathrm{C}}$ and $N_{\mathrm{BV}}$ easily exceed $10^{6}$ in the smallest-sized biologically relevant simulation runs, $\mathcal{O}(N_{\mathrm{C}}N_{\mathrm{BV}})\sim 10^{12}$ per time step, making the simulations utterly unmanageable. To reduce the computational complexity, we take into account that only the nearest blood vessels actually contribute to the inputs and outputs of the different molecules in the vicinity of a given cell trimming the large number of cell-blood vessel pairs to $\mathcal{O}(N_{\mathrm{C}})$. Moreover, the timesteps to follow cell biochemistry are of order 1-10 s while biological relevant timesteps for vascular remodelling are in the order of hours. Therefore we introduce two different time steps to avoid a huge number of useless sweeps over the blood vessel network (compare Figure 3).

Development of the microenvironment in small solid tumours

With the program described before, we studied selected aspects of the tumour microenvironment. To this end, we carried out simulation runs with three different tumour seeds:

1. at the center of a blood vessel network,
2. close to an arterial bifurcation,
3. close to a venous bifurcation. The details of each simulation are reported in the supplemental material. Here we display the results in a way that allows an easy comparison with measurements [38, 16].

As expected, the seeding is important in the very first steps of tumour growth: a seed close to blood vessels leads to fast growth and to a quick remodelling of the blood vessel network; for seeds far from blood vessels, the initial growth is slower and the subsequent initial development of the tumor microenvironment differs slightly.

Seed at the center of a blood vessel network

The computational approach provides an easy access to information that is hard to obtain with experiments and allows multiple ways of visualization (compare Figure 4). Here we take advantage of the information about the distance to the nearest blood vessel for every cell. We start with the distribution which is a defining feature of the microenvironment, as larger distances mean less oxygen and nutrients, and therefore starvation and death for many tumour cells. In Figure 5 we see how the distribution changes while the tumour grows. The distribution at the early stage ( 280 hours past seeding) ranges from 0 to 80 $\mu$m ensuring a optimal oxygenation to all cells and gets distorted until it covers a range of up to 160 $\mu$m at a later stage ( 580 hours past seeding) where cells no receive a sufficient input of oxygen and nutrients.

Intuitively the increase in the mean of cell-blood vessel distances leads to a decrease of available oxygen and nutrients, and to an increase of lactate ions and to a lower pH. Since the cell-blood vessel distance is known, we can display space-time data on the distributions of cell radii, pO2 and pH (Figure 6)

• the data confirms the intuition. We uploaded two videos where we visualize pO2 (fredrich:po2:2018) and pH (fredrich:pH:2018) during the the complete of the growth processes.

Finally, the changes in the microenvironment stretch the cell proliferation period and therefore modifies the cell cycle distribution as seen in Figure 7 where we can follow the formation of a necrotic core. At the beginning (280 hours past initial seed) the relative amount of cells in the different phases varies with the distance from the nearest vessel. Surprisingly this relative amount stays almost constant after the establishment of different niches. Compare with Figure 4 were we show the cell phase distribution in a more pictorial way. Even in such a small tumour we observe the initial development of necrotic regions which characterizes the development of the tumour microenvironment.

Actual tumours have irregular shapes, while the cells in VBL are not polarized, and when grown in uniform environments produce nearly spherical shapes. In the presented case the blood vessel network drives the growth in specific directions and we expect a marked deviation from sphericity already with small simulated tumours. Ideally, in a spherical tumour all cells on the surface have the same distance to the centroid. We find that the distribution of distances from the spheroid’s surface to the centroid a) broadens and b) deviates from symmetric form as the tumour grows which is a clear indication of the deformation of the cell cluster and of its deviation from sphericity.

Seed close to an arterial or venous bifurcation

Helminger et al. [16] measured the partial oxygen pressure and pH within different tumours. To compare our simulations with their experimental results, we seeded the initial tumour cell close to bifurcations, both arterial and venous, and sampled the cells along different lines (for details see supplemental material). Figure 9 is meant to illustrates the sampling procedure while we show data for a early and a late simulation state in Figure 10.

It is interesting to note that the left panel in first row in Figure 10 shows a dip in both pO2 and pH approximately midway between the two blood vessels. The right panel in the same row, taken at a later time, shows a maximum at roughly the same position because a new blood vessel sprouted into the intervening space. A detailed view of this dynamic change is shown in Figure 11, where the mean of the pO2 data is plotted for all intermediate time points as the new blood vessel grows.

Conclusions

As explained in the Introduction, the study of the ecology of the tumour microenvironment at the angiogenic switch is the primary motivation for the work described in this paper and for the development of the presented complex software tool. Indeed, the Figures 10 and 11 display a surprisingly large spatial and temporal variability already in a very small vascularized tumour. In particular, Figure 11 shows a continually changing and rugged landscape: this means that the niche diversity is large, and consequently that there is a high evolutionary pressure, and a high probability that the microenvironment selects different tumour clones even in small tumours. This is important in view of the long time required for tumour growth, as it has been estimated that a tumour requires about 10 years to reach a size of 1 cm [27]. According to the clonal selection hypothesis new and aggressive tumour clones develop through Darwinian selection during this extended growth time [27, 2]. The results discussed in this paper lend further credibility to the clonal selection hypothesis, and in our future studies we plan to demonstrate the spatial and temporal variations of different clone populations in this very complex tumour microenvironment. In order to highlight common and diverse evolutionary pathways in different cancers, we also plan to explore the adaptive evolution of cell clones in different microenvironments mimicking different organs or tissue types to reproduce the distinguishing features of selected tumour types.

Finally, we note that the molecular mechanisms that promote genotypic changes in tumour clones are well known and understood, but genotypic variability is only one feature of cancer’s evolutionary landscape. The other important feature is the variability of the environment that supervenes the genome and drives evolution itself, and with our computational tool we can start to explore its Darwinian dynamics [14].

Acknowledgments

The authors thank the Collaborative Research Center SFB 1027, the University of Trieste (FRA 2016), the University of Verona (RATs 2015), and the Deutscher Akademischer Austauschdienst (DAAD), for providing financial support to this project. Finally, we thank Benjamin Bogner for helping us with the computer graphics.

Author contributions statement

TF: conducted simulations, analyzed data, drafted results part of manuscript; EM: suggested scientific problem, added VBL to public domain, drafted other sections; TF and EM: designed an interface for VBL and Tumorcode, optimized parameters; RC: helped with biological aspects and parameter tuning; HR: developed Tumorcode program, supervised TF. All authors approve the content and reviewed the manuscript.

Additional information

**Supplementary information:

**Simulation details are described in a supplementary file.

**Source code availability:

**All source code containing the used parameters is distributed with this article as supplementary material and hosted at GitHub https://github.com/thierry3000/tumorcode.

**Data availability:

**We uploaded the full simulation raw data to zenodo with following digital object identifiers:

https://doi.org/10.5281/zenodo.2541688

https://doi.org/10.5281/zenodo.2541655

https://doi.org/10.5281/zenodo.2541678

https://doi.org/10.5281/zenodo.2541667

**Competing interests:

**The authors declare no competing interests.