## Abstract

Patients diagnosed with glioblastoma, an aggressive brain tumour, have a poor prognosis, with a median overall survival of less than 15 months. Vasculature within these tumours is typically abnormal, with increased tortuosity, dilation and disorganization, and they typically exhibit a disrupted blood–brain barrier (BBB). Although it has been hypothesized that the ‘normalization’ of the vasculature resulting from anti-angiogenic therapies could improve drug delivery through improved blood flow, there is also evidence that suggests that the restoration of BBB integrity might limit the delivery of therapeutic agents and hence their effectiveness. In this paper, we apply mathematical models of blood flow, vascular permeability and diffusion within the tumour microenvironment to investigate the effect of these competing factors on drug delivery. Preliminary results from the modelling indicate that all three physiological parameters investigated—flow rate, vessel permeability and tissue diffusion coefficient—interact nonlinearly to produce the observed average drug concentration in the microenvironment.

## 1. Background

With a median overall survival of less than 15 months [1], glioblastomas (GBMs) are challenging tumours with unsatisfactory responses to chemotherapy [2–5]. While it is known that many factors, from genetic variants [6] to microscopic tissue oxygenation [7], influence drug effectiveness, the drug delivery within the tumour is a fundamental consideration. Tumours often have abnormal, disorganized vasculature with increased tortuosity, shunting, poor perfusion and permeability, resulting in impaired blood flow. These abnormalities can result in elevated interstitial fluid pressure and hinder the delivery of therapeutic agents to tumours [8–11].

Poor or uneven regional distribution of therapeutic agents in the tumoural tissues of brain tumours, well known to be highly heterogeneous, is believed to contribute to disappointing clinical trial results. Further, complementary drugs such as anti-angiogenic agents, when used along with cytotoxic agents, modulate drug delivery in complex ways. Antiangiogenics have been hypothesized to improve drug delivery through ‘normalization’ of the vasculature [7,12,13]. However, it has also been suggested that the restoration of blood–brain barrier (BBB) integrity might limit the delivery of concomitant cytotoxic agents and hence their effectiveness [14,15].

Regional measurement of tissue drug delivery is thus a foundational measurement that has profound implications for the timing and dosing of drug combinations, yet one that remains elusive in clinical practice. The challenge lies in the complexity of the tumour microenvironment and the spatial and temporal heterogeneity of the tumours. Structural and functional changes in tumour vasculature as well as changes in the cellularity in the interstitial space can impair drug delivery [16]. Determining regional drug distribution in tumours in individual patients non-invasively is a crucial tool needed for improving our understanding of how best to utilize our current drug therapies and evaluate new treatment regimens in individual patients.

Therefore, in order to improve our understanding of how best to use drugs to target these (and other) aggressive tumours, the clinical community needs a tool to measure drug penetration and delivery to brain tumours in individual patients and longitudinally over time [17]. *In vivo* magnetic resonance imaging (MRI) is a non-invasive tool that has been demonstrated to interrogate microscopic tissue features sensitive to the microvascular space, tissue perfusion and interstitial delivery [7,12,18–21] but methods to use such MRI data to precisely predict therapeutic drug delivery distribution and kinetics directly do not currently exist. For example, figure 1 shows a post-contrast magnetic resonance (MR) image as well as a positron emission tomography (PET) image for the radiolabelled [11C] temozolomide for two patients before and after anti-angiogenic therapy. The post-contrast MR image is a map of the vascular permeability as the gadolinium agent extravasates into the interstitial areas leading to corresponding areas of enhancement. The PET image is similarly a spatial map of the distribution of the temozolomide, the chemotherapy agent. These images *suggest* that the drug delivery, at least in these patients, appears to be higher in areas of increased vascular permeability and then decreases with decreased permeability. However, it is difficult to test this further from the data available from the images alone. Nonetheless, we believe that MRI can indeed be the basis for such a predictive tool, with the advantages of its ubiquitous clinical availability and non-invasiveness. Further, we believe that MRI methods sensitive to the physical properties of the tumour microenvironment can be used to provide the needed inputs into a detailed, realistic mathematical model of the tumour microenvironment, designed and validated to accurately predict drug delivery. We refer to this approach as ‘Model + MRI’.

In this paper, we aim to provide a ‘proof of concept’ of our Model + MRI approach by formulating and validating a mathematical model of perfusion in a vascular network, the actual structure of which is empirically derived from high-resolution three-dimensional optical microscopic features of the vascular network in mouse models, both morphological and physiological. This will allow us to have highly realistic computational models of the essential features that govern local drug delivery and tissue distribution including the permeability of the drug to regions of disrupted and intact BBB, and the complexities of the interstitial space. The empirical data for these model inputs will come from advanced imaging methods that have been developed over the last several years [18,22–28].

## 2. Material and methods

We have previously developed a hybrid discrete-continuum mathematical model of flow through vascular networks allowing for dynamic adaptation of the vessels, with applications to tumour-induced angiogenesis and drug delivery, wound healing and retinal development [29–37]. Recent extensions to the model have incorporated permeability of the vasculature and diffusion of the extravasated agent through heterogeneous three-dimensional tissue.

Advances in optical imaging technology have enabled us to develop vascular anatomical network (VAN) models based on two-photon microscopy in mouse models [26,38]. These ‘virtual voxels’, representing realistic vascular networks, were used as input to the mathematical models to calculate the flow and distribution of agents within the ‘virtual voxel’.

We have preliminary imaging data for patients with glioblastoma enrolled in an ongoing clinical trial receiving an anti-angiogenic agent. Patients received MR as well as PET imaging before and after treatment with the anti-angiogenic therapy. In addition to standard anatomical MRI, we also performed dynamic perfusion and permeability imaging in order to measure blood flow, blood volume and vascular permeability at a voxel level. PET was performed using and radiolabelled [11C] temozolomide [39,40], thus allowing for the visualization of the spatial distribution of the drug (cf. figure 1). Parametric maps were generated for cerebral blood flow, volume and permeability [18] as well as for the standardized uptake value [41] of radiolabelled temozolomide and quantified within regions of interest identified by the clinician.

### 2.1. Modelling perfusion and tissue diffusion

#### 2.1.1. Model description

In this work, we model a small block of rodent brain tissue (0.6 × 0.6 × 0.6 mm) that hosts—and is perfused by—a highly tortuous interconnected network of flowing capillary elements (figure 2*a*). For simulation purposes, the network itself is partitioned into nodes (junctions) and edges (connections) and each element (node/edge) is assigned a range of intrinsic attributes (radius, length, vessel type, volume and conductivity). This model was created using structural images with fluorescein-labelled blood. A three-dimensional mask of the vasculature was obtained from the angiogram and the graph and a mesh of the vasculature was generated [38]. This provided an accurate three-dimensional vascular structure in which to carry out the flow simulations. Summary statistics for the size distribution of the vessel radii are given in table 1.

Several attributes can be extracted directly from optical imaging (two-photon and optical coherence tomography) of the rodent brain and these have been incorporated into our computational (i.e. *in silico*) vascular bed. We note that, in contrast with our earlier modelling studies, the network used here is completely lattice-free and inlet and outlet pressure boundary conditions can be assigned with great flexibility. While we have assumed a cylindrical geometry for capillary elements at present, it should be noted that this assumption can be easily relaxed if image data suggests otherwise.

The host tissue surrounding the vessels is modelled using a regular cubic mesh comprising 30 × 30 × 30 cubes, of equal edge length *hx, hy, hz*. Given the scale of the block of brain tissue (0.6 × 0.6 × 0.6 mm), this means that the spatial grid size is of order 20 µm (figure 2*b*). Beginning with a full mesh we then proceed to remove all blocks that lie completely within vessels and identify all vessel segments that provide capillary surface area to each remaining tissue element—this information is used later when calculating the diffusive flux from a vessel to an adjacent cube (and vice versa).

#### 2.1.2. Flow and diffusion

Computing the nodal pressure field during the full simulation is an essential component of the process and is used to update the tracer concentration within each capillary segment at each time step. We fix an injection rate for the blood flow (although constant pressure drop simulations can also be readily considered), which is effectively used as a boundary condition to determine the nodal pressure distribution and elementary flows within the system. We assume that for each vascular element of shape factor *G* (the ratio of cross-sectional area to the square of the perimeter), length *L* and cross-section *A*, the flow *Q* is given by a Poiseuille-type law
where *g* is the element conductance, *μ* is the fluid viscosity in the element and Δ*P* the pressure difference acting across the element.

By applying the appropriate pressure gradient across the network, the pressure field inside a network can be obtained by applying the mass conservation law at each node *i* (assuming that the flowing fluid is incompressible)
where *Q _{i}*

_{,j}is the flow between node

*i*and node

*j*. We use Cholesky factorization to solve this system of linear equations and determine the pressure value at each node. For this work, we assume the flow of a simple Newtonian fluid, although complex blood rheology can be incorporated straightforwardly (see [31] for full details).

In the absence of drug leakage across the vessel membrane, we can apply mass conservation at nodes to update the contrast agent concentrations within each capillary over time. The initial concentration is set to zero for all elements and we then inject a benign contrast agent into the inlet arterioles at dimensionless unit concentration.

If we consider the configuration depicted in figure 3, where *C _{i}* and

*Q*

_{i}_{,j}represent, respectively, the concentration and the flow in capillary element

*i*; the new tracer concentration after a time step Δ

*t*in capillary 1 is given by where

*C*

_{old}is the old tracer concentration and

*V*

_{1}the volume of capillary 1.

However, when tracer leakage in the surrounding tissue is considered, this approach needs to be amended to account for the diffusion of contrast agent into the host tissue. To satisfy mass conservation in this scenario, we need to consider (i) the convective flow of tracer in the vessels (as described above); (ii) the diffusive flux of agent from the vessels to the tissue (and vice versa); (iii) diffusive transport between tissue blocks; and (iv) tracer decay in the tissue.

These four main flow mechanisms are expressed as a coupled system of differential equations as described below
2.1
2.2where *S*V = vessel concentration; *V*_{v} = vessel volume; *Q*_{v} = flow in vessel V; MP_{Vi} = vessel membrane permeability (MP) between vessel *V* and tissue block *i*; *S*_{Ti} = tissue block *i* concentration, *A*_{Ti} = contact area between vessel *v* and tissue block *i*; *S*_{T} = tissue block *T* concentration; *V*_{T} = tissue block *T* volume; *D*_{T} = diffusion coefficient in tissue; *σ* = tissue intake coefficient.

The term describes the mass flow of drug carried to vessel *V* by its neighbouring vessels *V _{i}* via convection; refers to the mass flow of drug leaving vessel

*V*via convection; refers to the mass flow of drug leaving vessel

*V*via diffusion into the surrounding tissue;

*D*

_{T}Δ

*S*

_{T}is the expression of diffusive flow between tissue blocks; refers to the mass flow leaving the tissue block to its surrounding vessels via diffusion, and

*σS*

_{T}describes the mass flow of drug consumed by the tissue.

The elementary flows *Q*_{Vi} and *Q*_{V} are obtained after solving the pressure field in the vasculature network, and they are proportional to the blood flow rate *Q*, which is treated as an input. We also note that a homogeneous MP is assumed in the entire network.

Solving this coupled system using finite difference methods yields the contrast agent concentration in each vessel *v* and tissue block (*i*, *j*, *k*) at every time step Δ*t* as described in the following scheme

To ensure mass conservation of the drug in the vasculature system, the time step has to be chosen carefully to ensure that the dimensionless drug concentration in every vessel and tissue block remains within the interval [0, 1]. Therefore, we consider the case where every vessel and tissue block gets a drug inflow of concentration 1 from its surrounding vessels and tissue blocks whilst assuming a zero drug outflow. This yields the following expression for the time step:

In fact, for the computational results here, was always satisfied for the range of parameters we used.

## 3. Results

In order to determine the most important factors affecting tracer distribution within a tissue voxel, a range of simulations was undertaken by varying flow rate, transmural transport coefficient (i.e. vessel permeability), and extracellular matrix diffusion coefficient. A summary of the parameter values used in the simulations is given in table 2 (cf. [26,33,35]).

Images of tracer evolution for all parameter combinations explored are shown in figure 4, with each row corresponding to a given set of input parameters (see figure caption for parameter values).

We begin by investigating the combined effects of vessel permeability and ECM diffusion upon the *average vessel concentration* at a fixed flow rate of 10^{−8} m^{3} s^{−1} (this average is taken over all vessels in the voxel). Figure 5*a*, corresponding to slow diffusion through the ECM (*D*_{T} = 10^{−7} m^{2} s^{−1}), shows that contrast agent builds up rapidly within the capillary network following infusion. The build-up of agent in the bed is faster when the capillaries are less leaky, as expected: however, we observe that the vessel concentration asymptotes towards a steady-state value that is insensitive to vessel permeability. By contrast, under conditions of faster tissue diffusion (*D*_{T} = 10^{−6} m^{2} s^{−1}), we find that the steady-state vessel concentration increases with decreasing vessel permeability (figure 5*b*)—the rate of diffusive transport through the tissue block is now relatively high and so contrast agent is efficiently removed from the vasculature throughout its transit.

Corresponding graphs of the *average tissue concentrations* are shown in figures 6*a*,*b*. For slow tissue diffusion, the results exhibit the opposite trend to that observed for vessel concentration (i.e. as vessel permeability increases, average tissue concentration increases more rapidly). When tissue diffusion is fast, however, an interesting effect is observed: we find that the highest average level of contrast agent in the tissue corresponds to an intermediate value of vessel permeability. Moreover, it can be seen that all of the simulations yield tissue concentrations that are far smaller than those seen under slow diffusion conditions. It appears that this combination optimizes the vascular supply of contrast agent close to the inlet of the voxel, where it then diffuses rapidly through the tissue and is quickly re-captured via diffusive transport from tissue into low-concentration, downstream vessels; effectively removing it from the system.

Finally, we have also studied the effect of flow rate on tracer delivery to the tissue (figure 7). Here, we observe that average tissue concentrations asymptote to higher values at higher flow rates.

Hence, we can conclude that all three physiological parameters investigated thus far—flow rate, vessel permeability and tissue diffusion coefficient—interact nonlinearly to produce the observed average tracer concentration in the voxel. Preliminary results from the modelling indicate that delivery of the agent to the tumour bed is affected by flow, permeability and diffusion in the interstitial space. We find that reductions in the permeability of the vasculature, as would be typified through restoration of the BBB in response to anti-angiogenic therapy, could potentially reduce the delivery of the drug to the tumour, while increased flow could increase drug delivery.

## 4. Discussion and conclusion

Mathematical modelling combined with empirical data from *in vivo* models can provide insights into the role that the vascular architecture, permeability and microenvironment heterogeneity can play in drug distribution within a tumour bed during anti-angiogenic therapy. In humans, PET/MR with radiolabelled drugs (such as temozolomide) can be used to non-invasively validate these models and provide patient specific guidance for therapy. As described above, our modelling suggests that under certain conditions of vascular flow and diffusion through the interstitial space, a decrease in the permeability of the BBB, such as resulting from anti-angiogenic therapies, might actually result in reductions in the drug delivered to the tumours, despite normalization of the vasculature.

This proof-of-concept allows us to explore the use of mathematical models combined with the parameters measured using advanced MR images to predict drug delivery. This modelling will serve as an important link between our optical microscopic measurements of the tumour microenvironment (vascular network, permeability and accessible interstitial space) and the macroscopic MRI, and will ultimately validate our use of MRI to predict drug delivery. In the future, we hope to take our realistic computational models of the vascular and interstitial spaces (incorporating detailed empirical information on vascular perfusion, permeability and interstitial transport) [37–40], and then quantitatively test the ability of several advanced MRI approaches to assess these microenvironmental features from our non-invasive MR datasets [17]. We can then directly validate, in our animal models, the ability of the Model + MRI data to provide specific regional assessment of drug delivery, assessed at both microscopic levels statically using PET. Finally, this approach can be validated in our human patients with the ultimate goal to use patient specific MRI parameters of vascular permeability, flow and diffusion (all parameters easily obtain with non-invasive MRI) to inform the Model + MRI and predict the patient specific response to a therapeutic agent.

## Authors' contributions

A.B. and M.W. carried out the computational simulations, participated in data analysis and drafted the manuscript; Y.-F.Y. and C.C. helped with the data acquisition and analysis; patient data were obtained by T.T.B. and E.R.G. as part of a clinical trial; they also helped with the study design and manuscript; S.Mc., M.A.J.C., J.K.-C., T.D., D.B. and B.R. conceived of the study, designed the study, coordinated the study, participated in data analysis and helped draft the manuscript. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

S.M.c. and M.A.J.C. gratefully acknowledge support of EPSRC grant no. EP/N014642/1 (‘EPSRC Centre for Multiscale Soft Tissue Mechanics—With Application to Heart & Cancer’). M.A.J.C. would like to thank the Isaac Newton Institute for Mathematical Sciences for its hospitality during the programme ‘Coupling geometric PDEs with physics for cell morphology, motility and pattern formation’ supported by EPSRC grant no. EP/K032208/1.J.K.-C., B.R. and Y.-F.Y. gratefully acknowledge NIH/NCI grant no. U01CA154601, E.R.G. acknowledges NIH/NCI grant no. K23CA169021 and C.C. acknowledges NIH/NIBIB grant no. R01EB014894. T.T.B. was supported by NIH/NCI grant R01-R01CA129371-01.

## Footnotes

One contribution of 12 to a theme issue ‘Coupling geometric partial differential equations with physics for cell morphology, motility and pattern formation’.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.