Oxygen Transport in a Three-Dimensional Microvascular Network Incorporated with Early Tumour Growth and Preexisting Vessel Cooption: Numerical Simulation Study

Abstract

We propose a dynamic mathematical model of tissue oxygen transport by a preexisting three-dimensional microvascular network which provides nutrients for an in situ cancer at the very early stage of primary microtumour growth. The expanding tumour consumes oxygen during its invasion to the surrounding tissues and cooption of host vessels. The preexisting vessel cooption, remodelling and collapse are modelled by the changes of haemodynamic conditions due to the growing tumour. A detailed computational model of oxygen transport in tumour tissue is developed by considering (a) the time-varying oxygen advection diffusion equation within the microvessel segments, (b) the oxygen flux across the vessel walls, and (c) the oxygen diffusion and consumption within the tumour and surrounding healthy tissue. The results show the oxygen concentration distribution at different time points of early tumour growth. In addition, the influence of preexisting vessel density on the oxygen transport has been discussed. The proposed model not only provides a quantitative approach for investigating the interactions between tumour growth and oxygen delivery, but also is extendable to model other molecules or chemotherapeutic drug transport in the future study.

Yan Cai 1 and Jie Zhang 2 and Jie Wu 3 and Zhi-yong Li 1

Academic Editor:Weibo Cai

1, School of Biological Science and Medical Engineering, Southeast University, Nanjing 210096, China
2, Shanghai Pulmonary Hospital, Tongji University, Shanghai 200433, China
3, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiaotong University, Shanghai 200240, China

Received 5 October 2014; Revised 17 November 2014; Accepted 17 November 2014; 28 January 2015

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.

1. Introduction

There are a lot of controversies about exactly how a tumour is initiated, but it is generally known that the interactions between tumour cells with the host microenvironment have played a predominant role in the pathophysiological mechanisms of tumourigenesis and progression. Hypoxia is believed to be one of the important hallmarks of the abnormal metabolic microenvironments in malignant tumours [1]. Since the diffusion limit of oxygen in tissue is 100-200 μ m [2], the regions far from blood vessels become hypoxic. However, many cancers growing in a blood-supply-sufficient environment, such as breast cancer and glioma, were also found to form the hypoxic region in foci. This suggests that tumour blood vessels fail to provide adequate levels of oxygen for the proliferating tumour cells. A key reason for this phenomenon is the abnormal architecture and function of tumour microvasculature, including the heterogeneous spatial distribution, the uneven vessel diameters, the leaky vessel wall, the poor blood perfusion, and the low level of red blood cell (RBC) velocity [3]. The heterogeneity of tumour blood flow due to the disorganized, dysfunctional microvasculature causes hypoxic region in tumour tissue despite the presence of blood flow.

Hypoxic tumour cells can express high levels of angiogenic regulators including vascular endothelial growth factor (VEGF). The new vasculature in response to the upregulated angiogenic factors provides essential nutrients for rapid neoplastic expansion of tumours [4]. Except this generally known mechanism of neovasculature, the "cooption" of preexisting vessel network also plays a significant role in early tumour progression [5]. Holash et al. [5] studied early angiogenic events using the rat C6 glioma model. The results showed that even the smallest C6 gliomas at just 1 week after implantation was well vascularized by coopting of preexisting blood vessels. Further experiments [6-8] revealed that when a small amount of tumour cells was implanted into healthy tissue, they managed to coopt and migrate along host vessels, as well as producing many chemical substances, such as VEGF, Ang-1, and Ang-2, to change the microenvironment around the host vessels, which can induce immature changes in the host tissue vasculature, including vessel dilation, increased capillary permeability, and tortuosity [9]. The local microenvironment especially the oxygen concentration distribution is important in the switch from preexisting vessel cooption to angiogenic vascularization. Detailed mathematical models are needed to provide an analysis framework not only to quantify physiological parameters such as tissue oxygen diffusion and vessel wall oxygen permeability that are not directly measurable, but also to handle the dynamic oxygen delivery processes.

Mathematical models for oxygen transport in capillary-perfused tissue have been studied extensively from last century. As a famous pioneering study, Krogh [2] presented a simple model to describe oxygen delivery through uniformly spaced parallel capillaries and considered a single capillary and its surrounding tissue as a cylinder region. Since it is difficult to measure the tissue oxygen concentration in three-dimensions with a micron resolution experimentally, more and more mathematical modellings and simulations have been proposed to obtain a better quantitative understanding in the microvascular oxygen transport. Goldman et al. published a series of work [10-12] on oxygen delivery in skeletal muscle, which have shown the importance of many features neglected in the Krogh model, including the capillary tortuosity and anastomoses, nonlinear oxygen consumption in tissue, myoglobin-facilitated diffusion and nonlinear oxyhemoglobin dissociation in the RBCs and plasma. Their models include a detailed description of intravascular resistance to oxygen transport and are capable of incorporating realistic three-dimensional microvascular network geometries. Secomb et al. [13] developed a theoretical model based on a Green's function approach, for realistic three-dimensional network geometries derived from observations of skeletal muscle, brain, and tumour tissues. However, the influences of haemodynamics conditions on the oxygen delivery are poorly considered in the above models.

For the modelling of oxygen delivery in tumours, oxygen concentration field is generally introduced to simulate the dynamic interactions of tumour cells with the surrounding tissues, by considering one or more pathophysiological characteristics during the process of malignant tumours growth. The reaction-diffusion equation is often used to describe the distribution of oxygen concentration. However, the source term of the equation is difficult to handle without the inclusion of preexisting vessels. Anderson [14] assumed that oxygen production is proportional to the density of extracellular matrix (ECM). Some other studies [15, 16] included the production of oxygen only from the simulation region boundary as constant oxygen source, while the oxygen can only be used by tumour cells and transported by diffusion from the boundary. In our previous coupled model [17], the oxygen supply by microvessels was assumed to be a linear function of local blood haematocrit, which is a rough way to deal with the oxygen transport through vessel network. These models may be useful for the study of tumour cell growth in a given nutrient environment; they cannot provide the dynamic oxygen supply variation caused by the preexisting vessel network adaptation to tumour growth and its influence on tumour growth.

The main aim of this study is to establish a dynamic mathematical model of tissue oxygen transport by a preexisting three-dimensional microvascular network which provides nutrients for an in situ cancer at the very early stage of primary microtumour growth. At the same time, the expanding tumour will consume oxygen during its invasion to the surrounding tissues and cooption of host vessels. In addition, the preexisting vessel cooption, remodelling, and collapse are modelled by the changes of haemodynamic conditions due to the growing tumour. Based on our previous haemodynamical calculation in solid tumour [18], we develop a detailed computational model of oxygen transport in tumour tissue by considering (a) the time-varying oxygen advection diffusion equation within the microvessel segments, (b) the oxygen flux across the vessel walls, and (c) the oxygen diffusion and consumption within the tumour and surrounding healthy tissue. The current model will provide a framework for analyzing the dynamic responses of local oxygen environment to the early avascular tumour growth and host vessel cooption, and will enable us to study quantitatively the influence of oxygen transport on the initiation of primary microtumours by preexisting vessel cooption.

2. Method

In this section, we first describe the generation of initial preexisting vessel network in the model (Section 2.1). The haemodynamic calculation through the 3D microvasculature is explained in Section 2.2. The oxygen transport as well as other chemicals is modelled based on the obtained haemodynamic information (Section 2.3). The tumour growth and microvessels are coupled by the changes of haemodynamic and chemical environments. It is based on the facts as follows: (a) the dynamic changes of tumour microvessels influence the chemical environments especially the oxygen supply, which determine the behavior of tumour cells (Section 2.4) and (b) the changes of haemodynamic and chemical environments due to the growing tumour induce tumour vessel dilation, cooption, remodelling, and collapse (Section 2.5). The simulation algorithm is summarized in Section 2.6.

2.1. Preexisting Vessel Network

For the morphological analysis we consider vessel segments within a cube simulation domain Ω of 1 mm³ . A basic grid of 100 × 100 × 100 is generated uniformly in the cube with a distance of 10 μ m between the neighbouring nodes (Figure 1). The preexisting vasculature we designed in the basic model has a typical pattern of a normal arteriolar network. Parallel arterioles are distributed along [figure omitted; refer to PDF] axis with uniform distance of 10 grids, that is, 100 μ m. Three orders with different vessel diameters are defined to reflect the heterogeneous characteristic of these arterioles. The vessel diameter and proportion of each order are shown in Table 1. In addition, capillaries are designed stochastically as cross links between the arterioles, having diameters of 8 μ m.

Table 1: Vessel diameter and proportion of three orders of arterioles and capillaries in the model.

	Vessel diameter (μ m)	Proportion
Arterioles
Order 1	50	10%
Order 2	30	20%
Order 3	10	70%
Capillaries	8	NaN

(a) 3D preexisting microvascular network with typical pattern of a normal arteriolar network, including parallel distributed vessels with varying vessel diameter and capillaries for cross link. (b) The enlarged view of local microvascular network. Different colours represent different orders of arterioles and capillaries (arterioles order 1: black; arterioles order 2: blue; arterioles order 3: green; capillaries: red).

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

2.2. Haemodynamic Calculation

The haemodynamic model in this study is based on our previous work on the coupled modelling of intravascular blood flow with interstitial fluid flow [18, 22]. Briefly, the basic equation for the intravascular blood flow is the flux concentration and incompressible flow at each node. Flow resistance is assumed to follow Poiseuille's law in each vessel segment. The interstitial fluid flow is controlled by Darcy's law. The intravascular and interstitial flow is coupled by the transvascular flow, which is described by Starling's law. Blood viscosity is considered according to function of vessel diameter, local haematocrit, and plasma viscosity. In addition, vessel compliance and wall shear stress are correlated to vessel remodelling and vessel collapse, respectively, which will be explained in details in Section 2.5.

The main equations for blood flow calculation are as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the flow rate of each vessel segment, which amounts to zero at each node of the vessel network due to the assumption of flux conservation and incompressible flow. [figure omitted; refer to PDF] is the vascular flow rate without fluid leakage; [figure omitted; refer to PDF] is the transvascular flow rate. [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the length and radius of vessel segment. [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the intravascular pressure and the interstitial pressure of the vascular element. The total difference of [figure omitted; refer to PDF] from plane [figure omitted; refer to PDF] to [figure omitted; refer to PDF] is set to be 3.5 mmHg as the driving force of blood in the network (or the boundary condition). [figure omitted; refer to PDF] is the hydraulic permeability of the vessel wall. [figure omitted; refer to PDF] is the average osmotic reflection coefficient for plasma proteins; [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the colloid osmotic pressure of plasma and interstitial fluid.

The velocity of intravascular and interstitial flow satisfies [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the intravascular flow velocity and the interstitial flow velocity, respectively; [figure omitted; refer to PDF] is the hydraulic conductivity coefficient of the interstitium; [figure omitted; refer to PDF] is the surface area per unit volume for transport in the interstitium.

The distribution of red blood cells (RBCs) at microvascular bifurcation is calculated based on the approach proposed by Pries and Secomb [23]. The details of blood rheology simulation were described in Wu et al. [22].

From the haemodynamic simulation, we are able to obtain the intravascular flow velocity [figure omitted; refer to PDF] and the hematocritic [figure omitted; refer to PDF] in the microvessel network which are used in oxygen concentration calculation and vessel diameter which is used to estimate vessel remodelling and collapse.

2.3. Oxygen Concentration Calculation

The glioma cell and endothelial cell behaviours are coupled by the changes of the chemicals in the extracellular matrix (ECM), such as oxygen and matrix degradation enzymes (MDEs). The transport of these chemicals (oxygen and MDEs) is modelled by quasisteady reaction-diffusion equations. The ECM is treated as a continuum substance and can be degraded by MDEs, while the MDEs are governed by diffusion, produced by TCs and ECs, and decay of itself. The equations describing the interactions of TCs and ECs with ECM and MDE are [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the ECM and MDE concentration, separately. The [figure omitted; refer to PDF] and [figure omitted; refer to PDF] terms represent a tumour cell and an endothelial cell located at a node position [figure omitted; refer to PDF] . Their values are either 1 if a cell is present or 0 if it is not. [figure omitted; refer to PDF] is MDE diffusion coefficient, and, [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] are positive constants.

To obtain more realistic oxygen concentration field, the advection and diffusion of oxygen in the vessel network are introduced [20]. The computational space is separated into three domains to characterize three distinct physiological processes, which are the oxygen advection equation inside the vessel, the oxygen flux across the vessel wall and the free oxygen diffusion in the tissue.

Specifically, the oxygen transport inside the vessel is represented by the advection equation [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the free and bounded oxygen concentrations, respectively. [figure omitted; refer to PDF] denotes the intravascular blood velocity.

The equation between the hemoglobin bound oxygen and free oxygen is [figure omitted; refer to PDF] where [figure omitted; refer to PDF] denotes haematocrit obtained from the haemodynamic calculation; [figure omitted; refer to PDF] is the haemoglobin concentration within a red blood cell; [figure omitted; refer to PDF] is the haemoglobin oxygen saturation.

The free oxygen flux across the vessel wall satisfies Fick's law: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the tissue volume and the associated vessel wall area. [figure omitted; refer to PDF] is the oxygen flux which is obtained by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the free oxygen concentrations at the exterior and interior surfaces of the vessel, respectively; [figure omitted; refer to PDF] is the Bunsen solubility coefficient; [figure omitted; refer to PDF] is the vessel wall thickness. [figure omitted; refer to PDF] is the vessel wall permeability which is varied in different maturity level of vessel segments.

The interstitial fluid velocity is very slow due to the low interstitial pressure gradient in the tumour region. In fact, [figure omitted; refer to PDF] is almost 100 times smaller than [figure omitted; refer to PDF] in value according to the simulation results in our previous model. Therefore, we assume the free oxygen transported through the tissue space is governed by the oxygen diffusion equation which is not influenced by the interstitial fluid velocity [figure omitted; refer to PDF] .

Consider [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the tissue oxygen diffusion coefficient and [figure omitted; refer to PDF] is the consumption coefficient.

The initial condition of ECM density is set to be 1 and other chemicals' concentrations (oxygen and MDEs) are 0. No-flux boundary conditions are used in the simulation field. Since chemicals are transported much faster than the characteristic time for cell proliferation and migration, the chemicals' concentrations are solved to steady state at each time step of the simulation with an inner iteration step of 5 s.

2.4. Tumour Cell Phenotype

The probabilistic hybrid model for tumour cell growth is based on the previous work [17]. The 3D model is defined on a 100 × 100 × 100 grid to cover a 1 mm × 1 mm × 1 mm volume, so the grid length corresponds approximately to the size of a tumour cell, that is, 10 μ m.

We assumed three different phenotypes of glioma cells: the proliferating cells ( [figure omitted; refer to PDF] ), the quiescent cells ( [figure omitted; refer to PDF] ), and the necrotic cells ( [figure omitted; refer to PDF] ). Initially, we put 20 proliferating cells in the central area. Two thresholds of oxygen concentration for cell proliferation ( [figure omitted; refer to PDF] ) and cell survival ( [figure omitted; refer to PDF] ) are introduced to describe the effects of oxygen field on the tumour cell actions. To be specific, if there is enough oxygen ( [figure omitted; refer to PDF] ) and a space available, a tumour cell will proliferate to two daughter cells with a probability, defined as [figure omitted; refer to PDF] . [figure omitted; refer to PDF] is the tumour cell age range from 1 to [figure omitted; refer to PDF] and with an incremental 1 in each simulation time step. [figure omitted; refer to PDF] is the tumour cell proliferation time (set to be 9 hour, that is, 6 time steps). This assumption means an older tumour cell would be more likely to proliferate. One of the two daughter cells will replace the parent cell and the other cell will move to the neighbouring element. When the local oxygen concentration at a tumour site is below the cell survival threshold [figure omitted; refer to PDF] , the tumour cell is marked as necrotic cells and will not be revisited at the next time step. The necrotic cells have a probability of 20% to disappear to form available space for glioma cells and endothelial cells if they stay necrotic for more than 45 hours (30 time steps). When tumour cells satisfy the survival condition but there is no neighbouring space for them to proliferate, they will go to quiescent.

Each phenotype of tumour cell has different coefficients of the consumption rate of oxygen and the production rate of MDEs [17] (Table 2).

Table 2: Parameters of different phenotypes of glioma cells.

Phenotypes	MDE production	Oxygen consumption
Proliferating cells (P)	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
Quiescent cells (Q)	[figure omitted; refer to PDF] /5	[figure omitted; refer to PDF]
Necrotic cells (N)	[figure omitted; refer to PDF] /10	[figure omitted; refer to PDF]

2.5. Vessel Cooption, Remodelling and Collapse

Experimental and clinical studies both revealed that microvessel diameter increases in response to growth factors. Döme et al. [9] found that even a single tumour cell can induce radical changes in the host tissue vasculature in a mouse model of glomeruloid angiogenesis. In our model, we consider vessel dilation as the first signal of a preexisting vessel becomes an immature vessel. A vessel segment surrounded by active tumour cells will increase its radius [figure omitted; refer to PDF] with the rate 0.4 μ m/h if [figure omitted; refer to PDF] . At the same time, the permeability of vessel wall [figure omitted; refer to PDF] is increasing in a dilation vessel and satisfies [figure omitted; refer to PDF] where [figure omitted; refer to PDF] represents the vessel permeability value in tumour tissue.

For a preexist vessel, once vessel dilation occurs, the vessel segment was treated as an immature vessel and was allowed to vary in radius in response to the difference between intravascular pressure and interstitial pressure. The capillary compliance satisfies the empirical equation of Netti et al. [21]: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the origin radius of the capillary; [figure omitted; refer to PDF] is the compliance exponent; [figure omitted; refer to PDF] is the compliance coefficient. [figure omitted; refer to PDF] is the collapse pressure and represents the ability of a vessel segment remaining structure intact under the pressure difference.

Based on the above equation, when the immaturation level of the vessel segment become more serious, the [figure omitted; refer to PDF] increases due to the higher [figure omitted; refer to PDF] , which can induce vessel compressing. When vessel radius is compressed to below zero, the vessel segment collapses. Another factor leading to vessel collapse in tumours is the apoptosis of endothelial cells due to the disturbed blood flow or severe reduction of blood flow rate in the vessel. Wall shear stress (WSS) is used to estimate this kind of vessel regression, similar to our previous work [17]. The WSS of a vascular segment [figure omitted; refer to PDF] can be calculated as [figure omitted; refer to PDF] We assume that a circulated vessel, which is surrounded by the TCs, will collapse with a predefined probability if the WSS value in the vessel is less than a threshold [figure omitted; refer to PDF] ( [figure omitted; refer to PDF] is the mean WSS value of the vessels of Arterioles order 3). The probability is assumed to be proportional to the duration of low WSS in the vessel; that is, the longer the vessel experiences the low WSS, the more likely the vessel is to collapse if the criterion is satisfied.

2.6. Simulation Algorithm

Step 0. Initialize.

Step 0.1 . Create the 3D preexisting vessel network and initialize the model parameters.

Step 0.2. Put 20 proliferating cells with random ages from [figure omitted; refer to PDF] to [figure omitted; refer to PDF] in the central area.

Step 0.3. Set up fluid flow boundary conditions.

Step 1 . Haemodynamic calculation (see (1)-(2)), and obtain the flow information.

Step 2 . Calculate the chemical's concentration field.

Step 2.1 . Solve (3) to obtain the concentration distribution of ECM and MDE.

Step 2.2 . Solve (4)-(5) to obtain the oxygen transport inside the vessel using the intravascular blood velocity [figure omitted; refer to PDF] and haematocrit [figure omitted; refer to PDF] from Step 1.

Step 2.3 . Solve (6) and (7) to obtain the free oxygen flux across the vessel wall.

Step 2.4 . Solve (8) to obtain the free oxygen transported through the tissue space.

Step 2.5 . Update the oxygen concentration field through the simulation domain.

Step 3 . Determine the behavior of tumour cells according to the local oxygen concentration, the available space, and the cell age and update the tumour cell distribution.

Step 4 . Update the vessel network.

Step 4.1 . Vessel cooption and vessel radius and permeability changes (see (9)) in immature vessels.

Step 4.2 . Vessel remodelling according to (10) update the radius of vessel segments.

Step 4.3 . Certain vessel collapse based on the [figure omitted; refer to PDF] changing and WSS criterion (see (11)).

Step 5 . Go to Step 1.

Each time step increment ( [figure omitted; refer to PDF] ) corresponds to 1.5 hour. We will use nondimensional time unit instead of hours in the following results. All parameters used in the simulation are summarized in Table 3.

Table 3: Parameter values used in the simulation.

Parameter	Value	Description	Reference
[figure omitted; refer to PDF]	10 μ m	Lattice constant
[figure omitted; refer to PDF]	0.82	Average osmotic reflection coefficient for plasma proteins	Baxter and Jain (1989) [19]
[figure omitted; refer to PDF]	20 mmHg	Colloid osmotic pressure of plasma	Baxter and Jain (1989) [19]
[figure omitted; refer to PDF]	15 mmHg	Colloid osmotic pressure of interstitial fluid	Baxter and Jain (1989) [19]
[figure omitted; refer to PDF]	4.13 × 10-8 cm2 /mmHg s	Hydraulic conductivity coefficient of the interstitium	Baxter and Jain (1989) [19]
[figure omitted; refer to PDF] / [figure omitted; refer to PDF]	200 cm-1	Surface area per unit volume for transport in the interstitium	Baxter and Jain (1989) [19]
[figure omitted; refer to PDF]	10-9 cm2 s-1	MDE diffusion coefficient	Anderson (2005) [14]
[figure omitted; refer to PDF]	1.3 × 102 cm3 M-1 s-1	ECM degradation coefficient	Cai et al. (2011) [17]
[figure omitted; refer to PDF]	1.7 × 10-18 Mcells-1 s-1	MDE production by TC	Cai et al. (2011) [17]
[figure omitted; refer to PDF]	0.3 × 10-18 Mcells-1 s-1	MDE production by EC	Cai et al. (2011) [17]
[figure omitted; refer to PDF]	1.7 × 10-8 s-1	MDE decay coefficient	Anderson (2005) [14]
[figure omitted; refer to PDF]	1.27 × 10-15 μ mol/(μ m3 mmHg)	Bunsen solubility coefficient	Fang et al. (2008) [20]
[figure omitted; refer to PDF]	10-5 cm2 s-1	Oxygen diffusion coefficient	Anderson (2005) [14]
[figure omitted; refer to PDF]	6.25 × 10-17 Mcells-1 s-1	Oxygen consumption coefficient	Anderson (2005) [14]
[figure omitted; refer to PDF]	2.8 × 10-7 cm/mmHg s	Vessel permeability in tumour tissue	Baxter and Jain (1989) [19]
[figure omitted; refer to PDF]	6.5 mmHg	Vessel compliance coefficient	Netti et al. (1996) [21]
[figure omitted; refer to PDF]	0.1	Vessel compliance index	Netti et al. (1996) [21]

3. Results

3.1. Oxygen Concentration Distribution with Early Tumour Growth

Figure 2 shows the dynamic processes of oxygen transport during the tumour growth and invasion. The blue regions shown in the right column represent the tumour invasive areas. It is found that the tumour cells generally grow along the parallel arterioles of the preexisting vessel network, which can provide sufficient oxygen for tumour growth. At the same time, the growing tumour undergoes central necrosis, due to the oxygen consumption of proliferating tumour cells, resulting in the hypoxia in the tumour invasion region (shown in the left column in Figure 2). In the current model, the vessel cooption and remodelling are also considered as key factors to control the dynamic interactions between preexisting vessels and growing tumour. With the cooption and remodelling of preexisting vessels induced by the invasion of tumour cells, some microvessels become immature and surrounded by tumour cells to be intratumour vessels. Due to the compression of tumour cells on the vessel wall and the continuous reduction of WSS, these immature vessels (especially capillaries between arterioles) will undergo collapse. As a consequence, the surrounded tumour cells become necrotic and will invade to healthy tissue further, increasing the tumour invasive area and forming some clusters of surviving tumour cells (shown in the right column in Figure 2).

The distributions of oxygen concentration (left column) at plane [figure omitted; refer to PDF] and the tumour invasion region (right column) from [figure omitted; refer to PDF] view, at [figure omitted; refer to PDF] (a), [figure omitted; refer to PDF] (b), [figure omitted; refer to PDF] (c), and [figure omitted; refer to PDF] (d), respectively.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

3.2. Sensitivity to Preexisting Vessel Network

To access the sensitivities to the preexisting vascular network, we changed the microvessel density (MVD) to values of 0.5 time (Case [figure omitted; refer to PDF] ) and 1.5 time (Case [figure omitted; refer to PDF] ) of that in the baseline case (Case [figure omitted; refer to PDF] ). The distribution of integrated oxygen concentration [figure omitted; refer to PDF] along [figure omitted; refer to PDF] -axis at plane [figure omitted; refer to PDF] in three cases is shown in Figure 3. The dotted lines show the average oxygen concentrations through the plane. In the three cases, the average oxygen concentrations increase with the initial MVD and decrease with tumour growing. In addition, the reduction of the average [figure omitted; refer to PDF] is more significant in the less MVD case (Case [figure omitted; refer to PDF] ) than that in the baseline (Case [figure omitted; refer to PDF] ) and more MVD (Case [figure omitted; refer to PDF] ) cases, which suggests that the emergence of hypoxia-induced tumour angiogenesis will be targeted earlier in a blood-supply-deficient microenvironment.

Distribution of integrated oxygen concentration [figure omitted; refer to PDF] along [figure omitted; refer to PDF] -axis at plane [figure omitted; refer to PDF] in three cases of different initial MVD. The dotted lines represent the average oxygen concentrations. Case [figure omitted; refer to PDF] : the baseline model ((d) [figure omitted; refer to PDF] ; (e) [figure omitted; refer to PDF] ; (f) [figure omitted; refer to PDF] ); Case [figure omitted; refer to PDF] : a MVD value of 0.5 time of that in case [figure omitted; refer to PDF] ((a) [figure omitted; refer to PDF] ; (b) [figure omitted; refer to PDF] ; (c) [figure omitted; refer to PDF] ); Case [figure omitted; refer to PDF] : a MVD value 1.5 times larger than that in case [figure omitted; refer to PDF] ((g) [figure omitted; refer to PDF] ; (h) [figure omitted; refer to PDF] ; (i) [figure omitted; refer to PDF] ).

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

(e) [figure omitted; refer to PDF]

(f) [figure omitted; refer to PDF]

(g) [figure omitted; refer to PDF]

(h) [figure omitted; refer to PDF]

(i) [figure omitted; refer to PDF]

The oxygen concentration at the central region around [figure omitted; refer to PDF] is always below the average [figure omitted; refer to PDF] , which corresponds to the hypoxia in the tumour center. The hypoxia region will gradually expand with tumour growth, while the concentration gradient of oxygen will reduce at the same time. This indicates the oxygen concentration distribution in the modelling region towards equilibrium in the dynamics of oxygen consumption of tumour growth and oxygen production of preexisting vasculature. It is noteworthy that the oxygen distribution remains substantially unchanged after [figure omitted; refer to PDF] in Case [figure omitted; refer to PDF] (Figures 3(h) and 3(i)), which suggests the balance between oxygen consumption and production emerged earlier in a blood-supply-sufficient microenvironment.

3.3. Flow-Dependent and Flow-Independent Oxygen Transport

In the present study, flow-dependent oxygen transport was used instead of the simple treatment of oxygen in our previous studies [24] and many others in which the vessel was treated as point source of oxygen. A test case was designed in which the vessel segment was set as a point source of oxygen. The statistical results in Figure 4 show the proportion of oxygen supply to the tumour tissue by every arterioles order at the end of simulation ( [figure omitted; refer to PDF] ). Since the point source of oxygen is only related to the quantity of vessel segments in the test case, the smallest capillaries with largest number provide the highest percentage of oxygen (up to 35%). In the basic case, although the number of larger vessels (arterioles order 1) is the smallest in the four different sized vessels, the abundant blood perfusion allows them to be the primary provider of oxygen, but at limited perfusion locations. There are no noticeable differences between the two cases in the final results of tumour growth. However, flow-dependent oxygen transport can offer more realistic oxygen concentration field since the blood perfusion is known to be heterogeneous in the tumour tissue.

Figure 4: The proportion of oxygen supply to the tumour tissue by every arterioles order and capillaries in the basic case and the control case at [figure omitted; refer to PDF] .

[figure omitted; refer to PDF]

4. Discussion and Conclusion

In this work, we have proposed a dynamic mathematical modelling system to investigate the oxygen transport in a three-dimensional preexisting vessel network during the early tumour growth process. A 3D tree-like architecture network with different arterioles orders for vessel diameter is generated as preexisting vasculature in host tissue. To obtain more realistic oxygen concentration field, the advection and diffusion of oxygen in the vessel network are introduced. The computational space is separated into three domains to characterize three distinct physiological processes, which are the oxygen advection equation inside the vessel, the oxygen flux across the vessel wall, and the free oxygen diffusion in the tissue. The oxygen advection inside the vessel and the oxygen flux across the vessel wall are calculated based on the coupling haemodynamic environment including the intravascular blood flow and the interstitial fluid flow. In addition, the dynamic changes of vessel diameter and vessel wall permeability are introduced to reflect a series of pathological characteristics of abnormal blood perfusion in tumours such as vessel dilation, leakage, cooption, remodelling, and collapse.

The simulation focuses on the avascular phase of tumour development and stops before the emergence of angiogenesis phase. The results show the oxygen concentration distribution at different time points of tumour growth. In addition, the influence of preexisting vessel density on the oxygen transport has been discussed. In a case of blood-supply-deficient microenvironment, that is, preexisting vessel network with low MVD, the significant oxygen consumption will lead to the uneven distribution of oxygen concentration through the tumour tissue and eventually upregulate the hypoxia-induced growth factors such as VEGF to activate angiogenesis. Compared with the simple treatment of oxygen transport in which the vessel was treated as point source of oxygen, the modelling of flow-dependent oxygen transport can offer more realistic oxygen concentration field since the blood perfusion is known to be heterogeneous in the tumour tissue. The proposed model not only provides a quantitative approach for investigating the interactions between tumour growth and oxygen delivery, but also is extendable to model other molecules or chemotherapeutic drug transport in the future study.

Acknowledgments

This research is supported by the National Basic Research Program of China (973 Program) (no. 2013CB733800), the National Natural Science Foundation of China (no. 11302050 and no. 11272091), the Nature Science Foundation of Jiangsu Province (no. BK20130593), the Fundamental Research Funds for the Central Universities (Yan Cai, SEU), and the Priority Academic Program Development of Jiangsu Higher Education Institutions (Yan Cai, SEU).

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] A. L. Harris, "Hypoxia: a key regulatory factor in tumour growth," Nature Reviews Cancer , vol. 2, no. 1, pp. 38-47, 2002.

[2] A. Krogh The Anantomy and Physiology of Capillaries , Yale University Press, New York, NY, USA, 1922.

[3] D. Fukumura, R. K. Jain, "Tumor microvasculature and microenvironment: targets for anti-angiogenesis and normalization," Microvascular Research , vol. 74, no. 2-3, pp. 72-84, 2007.

[4] P. Y. Wen, S. Kesari, "Malignant gliomas in adults," The New England Journal of Medicine , vol. 359, no. 5, pp. 492-507, 2008.

[5] J. Holash, P. C. Maisonpierre, D. Compton, P. Boland, C. R. Alexander, D. Zagzag, G. D. Yancopoulos, S. J. Wiegand, "Vessel cooption, regression, and growth in tumors mediated by angiopoietins and VEGF," Science , vol. 284, no. 5422, pp. 1994-1998, 1999.

[6] P. Vajkoczy, M. Farhadi, A. Gaumann, R. Heidenreich, R. Erber, A. Wunder, J. C. Tonn, M. D. Mengerand, G. Breier, "Microtumor growth initiates angiogenic sprouting with simultaneous expression of VEGF, VEGF receptor-2, and angiopoietin-2," The Journal of Clinical Investigation , vol. 109, no. 6, pp. 777-785, 2002.

[7] M. Lorger, B. Felding-Habermann, "Capturing changes in the brain microenvironment during initial steps of breast cancer brain metastasis," The American Journal of Pathology , vol. 176, no. 6, pp. 2958-2971, 2010.

[8] C. Zhao, H. Yang, H. Shi, X. Wang, X. Chen, Y. Yuan, S. Lin, Y. Wei, "Distinct contributions of angiogenesis and vascular co-option during the initiation of primary microtumors and micrometastases," Carcinogenesis , vol. 32, no. 8, pp. 1143-1150, 2011.

[9] B. Döme, M. J. C. Hendrix, S. Paku, J. Tovari, J. Timar, "Alternative vascularization mechanisms in cancer: pathology and therapeutic implications," The American Journal of Pathology , vol. 170, no. 1, pp. 1-15, 2007.

[10] D. Goldman, A. S. Popel, "A computational study of the effect of capillary network anastomoses and tortuosity on oxygen transport," Journal of Theoretical Biology , vol. 206, no. 2, pp. 181-194, 2000.

[11] D. Goldman, A. S. Popel, "A computational study of the effect of vasomotion on oxygen transport from capillary networks," Journal of Theoretical Biology , vol. 209, no. 2, pp. 189-199, 2001.

[12] N. M. Tsoukias, D. Goldman, A. Vadapalli, R. N. Pittman, A. S. Popel, "A computational model of oxygen delivery by hemoglobin-based oxygen carriers in three-dimensional microvascular networks," Journal of Theoretical Biology , vol. 248, no. 4, pp. 657-674, 2007.

[13] T. W. Secomb, R. Hsu, E. Y. H. Park, M. W. Dewhirst, "Green's function methods for analysis of oxygen delivery to tissue by microvascular networks," Annals of Biomedical Engineering , vol. 32, no. 11, pp. 1519-1529, 2004.

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

[15] A. A. Patel, E. T. Gawlinski, S. . Lemieux, R. A. Gatenby, "A cellular automaton model of early tumor growth and invasion: the effects of native tissue vascularity and increased anaerobic tumor metabolism," Journal of Theoretical Biology , vol. 213, no. 3, pp. 315-331, 2001.

[16] S. C. Ferreira, M. L. Martins, M. J. Vilela, "Reaction-diffusion model for the growth of avascular tumor," Physical Review E: Statistical, Nonlinear, and Soft Matter Physics , vol. 65, no. 2, 2002.

[17] Y. Cai, S. Xu, J. Wu, Q. Long, "Coupled modelling of tumour angiogenesis, tumour growth and blood perfusion," Journal of Theoretical Biology , vol. 279, no. 1, pp. 90-101, 2011.

[18] J. Wu, S. Xu, Q. Long, "Coupled modeling of blood perfusion in intravascular, interstitial spaces in tumor microvasculature," Journal of Biomechanics , vol. 41, no. 5, pp. 996-1004, 2008.

[19] L. T. Baxter, R. K. Jain, "Transport of fluid and macromolecules in tumors. I. Role of interstitial pressure and convection," Microvascular Research , vol. 37, no. 1, pp. 77-104, 1989.

[20] Q. Fang, S. Sakadzic, L. Ruvinskaya, A. Devor, A. M. Dale, D. A. Boas, "Oxygen advection and diffusion in a three-dimensional vascular anatomical network," Optics Express , vol. 16, no. 22, pp. 17530-17541, 2008.

[21] P. A. Netti, S. Roberge, Y. Boucher, L. T. Baxter, R. K. Jain, "Effect of transvascular fluid exchange on pressure-flow relationship in tumors: a proposed mechanism for tumor blood flow heterogeneity," Microvascular Research , vol. 52, no. 1, pp. 27-46, 1996.

[22] J. Wu, Q. Long, S. Xu, A. R. Padhani, "Study of tumor blood perfusion and its variation due to vascular normalization by anti-angiogenic therapy based on 3D angiogenic microvasculature," Journal of Biomechanics , vol. 42, no. 6, pp. 712-721, 2009.

[23] A. R. Pries, T. W. Secomb, "Microvascular blood viscosity in vivo and the endothelial surface layer," American Journal of Physiology-Heart and Circulatory Physiology , vol. 289, no. 6, pp. H2657-H2664, 2005.

[24] Y. Cai, J. Wu, S. Xu, Z. Li, "A hybrid cellular automata model of multicellular tumour spheroid growth in hypoxic microenvironment," Journal of Applied Mathematics , vol. 2013, 2013.