Fine-grained simulations of the microenvironment of vascularized tumours

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.

fine-grained simulations of the microenvironment of vascularized tumours thierry fredrich 1 , Heiko Rieger 1 , Roberto chignola 2 & edoardo Milotti 3 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.
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][2][3][4][5][6] .
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 μm, and the regions in between them can become highly hypoxic 4 .
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. [7][8][9][10][11][12][13] , and the reviews 14,15 .
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 16,17 , 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 19 . 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.  . 500 μm thick slice trough tumour vasculature of 1 cm lateral size. For 3D illustration we extent vessels with radius bigger than 20 μm for another 750 μm above and below the slice. According to the lower panels: arteries are displayed in varying reddish color and veins in varying blueish color (unit is also μm). Capillaries and young vessel are shown in green.

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,[20][21][22][23] , for Tumorcode, see [24][25][26] ).
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  www.nature.com/scientificreports www.nature.com/scientificreports/ 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 pO 2 concentration and by pH.
• Weighting of molecular paths and cellular mechanisms is possible. The numerical methods used in VBL have been described elsewhere 27 , 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 28 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) 30 . 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 21,27 . 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 31 .
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. 24,31,32 ), here www.nature.com/scientificreports www.nature.com/scientificreports/ we use the simulation framework Tumorcode described in 25 , 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 www.nature.com/scientificreports www.nature.com/scientificreports/ experimental data and successfully reproduced: 1) morphology and flow characteristics 33,34 , 2) interstitial fluid flow and pressure 35 , and 3) tissue oxygenation 26 for both, healthy and tumour tissue.
Again, here is a short list of its main features of Tumorcode: www.nature.com/scientificreports www.nature.com/scientificreports/ • Hemodynamics includes the phase separation effect (Fåhraeus-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) 36 . • Angiogenesis, i.e., the addition of new segments, is driven by the VGEF gradient and is partly stochastic 24 .
• 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 24 . • A low local VGEF growth factor concentration produces blood vessel regression 24 .
• 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.  www.nature.com/scientificreports www.nature.com/scientificreports/ • The blood vessel hematocrit acts as source for oxygen dissolved in tissue and the tissue simultaneously consumes the oxygen. In 26 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 Figs 1 and 2, and the fact that the tumour is a continuous structure (governed by reaction-diffusion equations 25 ) 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 www.nature.com/scientificreports www.nature.com/scientificreports/ equations. We interpolated cell properties from the cells position to the surrounding nodes of the cubic lattice (compare Fig. 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 N N ( ) C BV  number of operations, where N C is the number of cells and N BV is the number of blood vessels segments. Since both N C and N BV easily exceed 10 6 in the smallest-sized biologically relevant simulation runs,  ∼ N N ( ) 10 C BV 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 N ( ) 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 Fig. 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 37,38 .
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 Fig. 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 Fig. 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 μm ensuring a optimal oxygenation to all cells and gets distorted until it covers a range of up to 160 μm at a later stage (580 hours past seeding) where cells not 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, pO 2 and pH (Fig. 6) -the data confirms the intuition. We uploaded two videos where we visualize pO 2 (fredrich:po2:2018) and pH (fredrich:pH:2018) during 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 Fig. 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 Fig. 4 were we show the cell phase distribution in a more pictorial way. Even in such a www.nature.com/scientificreports www.nature.com/scientificreports/ 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 (Fig. 8) 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. 38 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 Fig. 10.
It is interesting to note that the left panel in first row in Fig. 10 shows a dip in both pO 2 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 Fig. 11, where the mean of the pO 2 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 Figs 10 and 11 display a surprisingly large spatial and temporal variability already in a very small vascularized tumour. In particular, Fig. 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 39 . According to the clonal selection hypothesis new and aggressive tumour clones develop through Darwinian selection during this extended growth time 39,40 . 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 16 .

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.