Mathematical and computational models of drug transport in tumours.

Downloaded from rsif.royalsocietypublishing.org on March 27, 2014

Mathematical and computational models of drug transport in tumours C. M. Groh, M. E. Hubbard, P. F. Jones, P. M. Loadman, N. Periasamy, B. D. Sleeman, S. W. Smye, C. J. Twelves and R. M. Phillips J. R. Soc. Interface 2014 11, 20131173, published 12 March 2014

Supplementary data

"Data Supplement" http://rsif.royalsocietypublishing.org/content/suppl/2014/03/10/rsif.2013.1173.DC1.htm l

References

This article cites 33 articles, 6 of which can be accessed free

Subject collections

Articles on similar topics can be found in the following collections

http://rsif.royalsocietypublishing.org/content/11/94/20131173.full.html#ref-list-1

computational biology (268 articles)

Email alerting service

Receive free email alerts when new articles cite this article - sign up in the box at the top right-hand corner of the article or click here

To subscribe to J. R. Soc. Interface go to: http://rsif.royalsocietypublishing.org/subscriptions


Mathematical and computational models of drug transport in tumours rsif.royalsocietypublishing.org

C. M. Groh1,2, M. E. Hubbard3, P. F. Jones5, P. M. Loadman6, N. Periasamy6, B. D. Sleeman2,7, S. W. Smye4,8, C. J. Twelves9 and R. M. Phillips6 1

Research Cite this article: Groh CM, Hubbard ME, Jones PF, Loadman PM, Periasamy N, Sleeman BD, Smye SW, Twelves CJ, Phillips RM. 2014 Mathematical and computational models of drug transport in tumours. J. R. Soc. Interface 11: 20131173. http://dx.doi.org/10.1098/rsif.2013.1173

Received: 16 December 2013 Accepted: 20 February 2014

Subject Areas: computational biology Keywords: computational modelling, mathematical modelling, drug delivery, drug transport, drug binding, pharmacokinetic profiles

Author for correspondence: M. E. Hubbard e-mail: [email protected]

Electronic supplementary material is available at http://dx.doi.org/10.1098/rsif.2013.1173 or via http://rsif.royalsocietypublishing.org.

Klinik und Poliklinik fu¨r Strahlentherapie (Medizinische Physik), Universita¨tsklinikum Wu¨rzburg, Josef-Schneider-Strasse 11, 97080 Wu¨rzburg, Germany 2 School of Mathematics, 3School of Computing, and 4Division of Medical Physics, University of Leeds, Leeds LS2 9JT, UK 5 Leeds Institute of Biomedical and Clinical Sciences, University of Leeds, Wellcome Trust Brenner Building, St James’s University Hospital, Leeds LS9 7TF, UK 6 Institute of Cancer Therapeutics, University of Bradford, Bradford BD7 1DP, UK 7 Division of Mathematics, University of Dundee, Dundee DD1 4HN, UK 8 R&D Department, Leeds Teaching Hospitals NHS Trust, Leeds LS2 9LN, UK 9 Leeds Institute of Cancer and Pathology, University of Leeds, St James’s Hospital, Leeds LS9 7TF, UK The ability to predict how far a drug will penetrate into the tumour microenvironment within its pharmacokinetic (PK) lifespan would provide valuable information about therapeutic response. As the PK profile is directly related to the route and schedule of drug administration, an in silico tool that can predict the drug administration schedule that results in optimal drug delivery to tumours would streamline clinical trial design. This paper investigates the application of mathematical and computational modelling techniques to help improve our understanding of the fundamental mechanisms underlying drug delivery, and compares the performance of a simple model with more complex approaches. Three models of drug transport are developed, all based on the same drug binding model and parametrized by bespoke in vitro experiments. Their predictions, compared for a ‘tumour cord’ geometry, are qualitatively and quantitatively similar. We assess the effect of varying the PK profile of the supplied drug, and the binding affinity of the drug to tumour cells, on the concentration of drug reaching cells and the accumulated exposure of cells to drug at arbitrary distances from a supplying blood vessel. This is a contribution towards developing a useful drug transport modelling tool for informing strategies for the treatment of tumour cells which are ‘pharmacokinetically resistant’ to chemotherapeutic strategies.

1. Introduction A characteristic feature of solid tumours is the presence of a poorly organized and dysfunctional vasculature. The blood vessels that develop in response to angiogenic stimuli are structurally and functionally different from those of normal tissues, leading to poorly perfused areas of the tumour [1]. This results in the establishment of a microenvironment that can have profound effects on tumour biology and response to chemotherapy. The tumour microenvironment is characterized by gradients of oxygen tension, nutrient status, catabolite concentrations, extracellular pH, cell proliferation rates and a multitude of biochemical changes that enable cells to adapt to these hostile conditions [2], unless conditions become so extreme that tumour cells cannot survive and regions of necrosis occur. Where delivery of cancer drugs to cells within the tumour microenvironment is impaired, these cells are ‘pharmacokinetically resistant’; this form of resistance, distinct from cellular resistance, is increasingly recognized as a barrier to effective treatment [3]. In this paper, we use the generic term ‘chemotherapy’ to cover both conventional cytotoxic drugs and novel ‘targeted’ therapies, because the challenges of drug delivery apply to both. The variations in microenvironment in vivo are extremely complex, and depend on distance from a supporting blood vessel. Important insights can,

& 2014 The Author(s) Published by the Royal Society. All rights reserved.


(a)

2

(b)

BV

N N

Glut-1

(d)

rosis nec

BV B V

Figure 1. HCT116 human colorectal tumour xenografts. (a) Histological sections at low magnification showing blood vessels (BV) surrounded by cords of viable cells and necrotic areas (N) further distanced from the vessels. (b) Higher magnification of a blood vessel and surrounding viable area. (c) Section through a clinical sample of a transitional cell carcinoma of the bladder immunostained for Glut-1, a glucose transporter which is upregulated under hypoxic conditions. Glut-1positive staining, visible as darker, denser regions (brown membrane staining), can be seen distal from the supporting blood vessel (BV). (d ) Schematic of the tumour cord where the central blood vessel (BV) is surrounded by viable cells, but as distance from the BV increases (as indicated by the arrow) the microenvironment becomes more extreme, leading to regions of necrosis. Cells that reside some distance away from the BV are ‘pharmacokinetically resistant’ to chemotherapy as drugs cannot effectively penetrate through multi-cellular layers. The solid scale bars represent 500 mm (a,c) and 100 mm (b). (Online version in colour.)

however, be gained from a much simpler representation, that of a ‘tumour cord’, where a ‘collar’ of cells surrounds a supporting blood vessel, with cells at a distance from the vessel comprising regions of necrosis (figure 1). This model forms the biological platform for the studies described within this paper, which focuses on predicting the ability of drugs to penetrate through several layers of cells from a central blood vessel, and reach more distant cells. This ability of drug to reach cells some distance from the vessel is key to the effectiveness of chemotherapy, for both cytotoxics and targeted drugs, because these cells are typically resistant to current cytotoxic treatment, mainly owing to inadequate drug penetration. The interplay between cellular factors and drug transport is complex, but drug delivery can be broken down into three stages: drug delivery to a tumour is determined by the supply of drug via the blood vessels in tumours; the flux or movement of drug through the tumour mass; and sequestration, which includes binding to cellular or extracellular components and metabolism of the drug [3]. The impact of each stage on drug delivery will vary depending upon the pharmacology of individual drugs and the biological properties of individual cancers (e.g. differential expression of targets or drug-metabolizing enzymes). It is recognized from the outset that the models describing these processes represent considerable simplifications of complex biology and geometry [4], which may be better described by a more sophisticated computational framework [5]. However, we believe simpler mathematical models play an important role in developing more complex schemes by

providing what might be described as ‘semi-quantitative’ understanding of the biological transport mechanisms, and offering a simple approach to assessing the relative merits of different protocols. Changing drug administration protocols varies the concentration and exposure time of drugs within the central blood vessel (both in practice and in our tumour cord models), and these models allow us to investigate how these factors influence drug delivery. To test a range of scheduling options in purely experimental animal tumour models is expensive and conflicts with the aim of reducing, refining and replacing (the 3Rs) animal experiments that is currently promoted by institutions such as the Medical Research Council. In silico modelling offers the promise of being able to test multiple experimental scenarios and streamline the search for drug treatment regimens that optimize drug delivery to tumour cells throughout the tumour microenvironment. There exists a plethora of models describing the transport of drugs in tissue, ranging from compartmental models that account for exchange of drug within spatially distinct intracellular compartments [6–9] to continuum models describing the transport over macroscopic tissue scales [10–14]. If modelling is to have greater predictive impact on the development of new therapeutic agents, then it is important that the relative merits and limitations of these different descriptions are clearly understood. The key questions this paper seeks to address are ‘Does each of these models give similar results for the variation in drug concentration in the tumour cord?’ and ‘How do the administration schedule and cell response affect drug

J. R. Soc. Interface 11: 20131173

(c)

rsif.royalsocietypublishing.org

BV


Three distinct modelling approaches are considered, each of which represents drug delivery from a central blood vessel to a surrounding tumour cord. Each model assumes axial uniformity—the dependence of the drug concentration on the distance from the central vessel does not vary along the vessel, to reduce the complexity (as in reference [16]). A genuinely multi-dimensional model is developed in §2.2, followed by two simplified, one-dimensional models (§§2.3 and 2.4) in which radial symmetry is assumed. The interaction of the chemotherapeutic agent with the microenvironment is restricted to drug binding only. This binding model, which is common to all three approaches, is discussed in §1. The in vitro experiments used to parametrize the model were conducted over relatively short time-scales and contained no elimination mechanism for the drug: explicit modelling of decay or elimination (other than clearance via the supplying vessel) is left to future work.

2.1. Binding model The interaction of the chemotherapeutic agent with the microenvironment of cells is described by a three-compartment model, composed of extracellular space (volume V1) with

k1

C2

V2

V1

3 with C3 £ C0

Figure 2. A three-compartment model of drug distribution in tissue. C1 represents the extracellular drug concentration, C2 the free intracellular drug concentration and C3 the bound intracellular drug concentration. (Online version in colour.)

a concentration C1 of free drug and intracellular space (volume V2) with concentrations C2 and C3 corresponding to free and bound drug, respectively (where, in this model, the term bound includes both DNA-intercalated drug and drug bound to the cell in other ways). The binding is described by a simple kinetic model: drug binds reversibly to sites within the cell. Applying the principle of mass action leads to three coupled ordinary differential equations which describe the system, V1 V2 and

dC1 ¼ ak1 (C2 C1 ), dt

dC2 ¼ ak1 (C1 C2 ) V2 k2 C2 (C0 C3 ) þ V2 k2 C3 dt V2

dC3 ¼ V2 k2 C2 (C0 C3 ) V2 k2 C3 , dt

(2:1) (2:2) (2:3)

in which k1 is the rate constant for the transmembrane transport of drug, a is the area of the interface between the extracellular and intracellular spaces (the surface area of the cell), k2 and k22 are the drug association and dissociation rates, respectively, and C0 is the concentration of binding sites within the cell. It is assumed that mixing in each compartment is instantaneous— that is, intracellular and extracellular diffusion are assumed to be fast on the scale of an individual cell. This model is illustrated by the two-dimensional schematic in figure 2. Note that the model presented in (2.1)–(2.3) is an extension of those in [12] and [17], in which drug binding is non-saturable and nonreversible. It would be straightforward to modify the binding model in this way, or to account for Michaelis–Menten-type transport should this be required, as in [8,18,19], for example. Values for the kinetic rate constants for the binding process are derived from a bespoke experimental binding assay. Doxorubicin binding to DLD-1, a colorectal adenocarcinoma cell line, was studied by incubating 106 tumour cells suspended in tissue culture medium (378C) with a range of doxorubicin concentrations (0–100 mM). At evenly distributed times between 0 and 2 h, a fixed volume of the culture was removed, and cells pelleted by centrifugation. Free doxorubicin in the supernatant was extracted and measured by a sensitive and specific high-performance liquid chromatography technique [20]. Binding was calculated by comparison with cell-free doxorubicin solutions incubated in parallel. The data from these experiments are included in the electronic supplementary material. DLD-1 was also the cell line used to estimate the transport rate, k0, taken from [6]. Two sets of experiments were performed, providing both steady-state and time-dependent data. A chi-squared minimization, described in the electronic supplementary material, was then performed to provide estimates for the parameters k1, k2, k22 and C0, given in table 1.


2. Models

C1

extracellular space intracellular space k2 C3 k–2


delivery?’ Where differences do appear, we will seek to explain the reasons for them and the consequences for the choice of modelling approach. Three approaches to modelling the spatio-temporal evolution of drug concentrations in a tumour cord are compared, each of which is representative of a class of models: (i) a multi-dimensional cell-centre model that defines a network of nodes (each node corresponding to a computational cell which is identifiable with a biological cell), in which drug transport is defined locally between nodes and their nearest neighbours; (ii) a compartmental model, which makes use of the concentric-layer structure of tumour cords; and (iii) a continuum model that assumes Fickian diffusion in the cylindrical geometry of the cord. The first of these approaches is amenable to multi-scale modelling [5,15], because each node may be characterized by a bespoke microenvironment consisting of, for example, a cell cycle and molecular pathways. The remaining models are tailored to the tumour cord geometry, so are less flexible but much simpler (and faster) computationally. In §2, after outlining the underlying binding model, which is parametrized by experimental data for the cytotoxic drug doxorubicin, a description of each spatio-temporal model is given, emphasizing the relationship between the three discrete transport models. In §3, the model predictions are compared for two scenarios: in the first, each model is tested using a single set of model parameters (and hence a single homogeneous biological environment) estimated from bespoke in vitro experimental data, allowing us to investigate the influence that the choice of mathematical approach to drug transport has on the predictions (with the same binding model). The second scenario explores the effect on the model predictions of varying the pharmacokinetic (PK) profile and model parameters representing the drug’s binding affinity, a biological characteristic which can, to some degree, be controlled by the administration of other drugs. The results are then discussed in §4 in relation to the choice of modelling approach.


variable

value

l

1.6 1025 m

source of estimate

vessel radius

histology [21]

L r

2.0 10 m 1.0 1025 m

cord radius (vessel þ approx. nine cells) cell radius

histology [22,23] histology [24]

d a k0 k1 k2 k22 kv D C0

0.0625 1.94028 105 m21

extracellular – intracellular volume ratio parameter membrane surface– tissue volume ratio

histology [6,25] pffiffiffiffiffiffiffiffiffiffiffi a ¼ 2=(r 1 þ d)

2.5 1026 m s21 1.0 1026 m s21

permeability between cells permeability across cell membrane

[6] experiment

0.90 1026 mM21 s21

drug association rate

experiment

14.0 10 s 1.25 1027 m s21

drug dissociation rate permeability across vessel wall

experiment estimated

5.0 10211 m2 s21 2.6 103 mM

interstitial diffusion rate binding site concentration

D ¼ 2k0r experiment

25 21

2.2. Multi-dimensional cell-centre model The most general approach considered here is a multidimensional model in which each biological cell in the cord is represented as a computational node (cell centre) characterized by geometrical and biological information: position, cell radius, binding rates and concentration of binding sites. Each of these nodes in isolation acts in accordance with the threecompartment binding model illustrated in figure 2, providing the potential to insert a bespoke model of the microenvironment for each node/cell. The binding model is augmented by a description of the spatial transport of the drug, described locally by defining transport terms between nodes and their nearest neighbours. A representative geometry of a slice through the tumour cord is illustrated in figure 3 (cf. figure 1b). We assume uniformity in the axial direction—variations along the vessel are assumed negligible—so only two-dimensional cord cross sections are considered in this paper. The terms volume and area are used as though the model were fully three dimensional but, for the sake of simplicity, the factor of the cord length is omitted from all of the analysis, because it cancels. Each computational node (see figure 3) is assigned associated interstitial and intracellular volumes, d1Vi and d2Vi, respectively, where d1 and d2 are volume fractions, parametrized by d in table 1 (see the electronic supplementary material for full details) and three concentrations: interstitial, intracellular free and intracellular bound drug (i) (i) denoted by C(i) 1 , C2 and C3 , respectively. The transport between nodes is assumed to be proportional to the concentration difference and the interface area, Aij, between connected nodes. Under these assumptions, the spatial variation represented by introducing distinct computational nodes can

be included in the model by adapting equations (2.1)–(2.3) to include terms for transport between nodes and having a separate set of equations for each node, i.e.

d1 Vi

X X dC(i) (j) (i) 1 Aij k0 (C1 C(i) Aij kv (C(j) ¼ v C1 ) 1 )þ dt j[V j[N i

i

(i)

þa

d 2 Vi

k1 (C(i) 2

C(i) 1 ),

dC(i) (i) 2 ¼ a(i) k1 (C(i) 1 C2 ) dt (i) (i) d2 Vi {k2 C(i) 2 (C0 C3 ) k2 C3 }

and

d2 Vi

(2:4)

dC(i) 3 dt

(i) (i) ¼ d2 Vi {k2 C(i) 2 (C0 C3 ) k2 C3 }:

(2:5) (2:6)

Here, Aij is the area of the interface between nodes i and j (see the electronic supplementary material) and a (i) is the area of the interface between the intra- and extracellular space surrounding node i. N i and V i are sets of indices of, respectively, the nodes and blood vessels neighbouring node i: all simulations presented in this paper have been carried out with a single central vessel. Cv(t) is a predefined PK profile in the blood vessel and determines the flux of drug through the vessel wall. At the outer boundary, this model automatically implies zero flux of drug, because there are no connections to nodes outside the cord geometry. This models a situation in which the tumour cord is surrounded by similar cords, providing a symmetry boundary condition. It is assumed, in this model, that the transport of drug is limited to the interstitium—drug is only exchanged between neighbouring nodes via the compartments representing extracellular space. Internalization and binding terms in equations (2.4) –(2.6) are analogous to those for the binding model. The additional rate parameter in this model is the


24

description

4


Table 1. Summary of model parameter values for baseline studies. In the final column, ‘experiment’ refers to the fitting to experimental data described in §2.1 and the electronic supplementary material, and ‘histology’ indicates estimation from histological tissue images, such as those illustrated in figure 1 or at www. virtualpathology.leeds.ac.uk, or from the cited references. The parameter k0 has been estimated based on the value of r1 (transport rate between cell layers) in the multi-layer model of [6]. The parameter kv has been chosen to give slower transport across the vessel wall than across the cell membrane, though in the disordered, leaky tumour vasculature this is likely to be highly variable. Note that the volume fractions used in §2 are defined by d1 ¼ d/1 þ d and d2 ¼ 1/1 þ d, where they are combined with the relevant compartmental volumes, Vi. We chose these to match the volumes of biological cells, though this is not necessary.


Vi ¼ p((2i 1)d þ 2l)d :

(2:10)

Ai ¼ 2p(l þ id):

Figure 3. A representative tumour cord geometry for the two-dimensional cell-centre model in which the computational nodes (cell centres) surround a central vessel. The radius associated with each node is illustrated by the surrounding circle and the connectivity is designated by the connecting lines. (Online version in colour.)

(2:11)

Note the similarity between compartmental and cellcentre equations. In fact, equations (2.7) –(2.9) can be viewed as a special case of equations (2.4) –(2.6), in which nodes are replaced by shells, and the concentrations within each shell represent the averages over the conglomerate of biological cells it contains.

2.4. Radially symmetric continuum model spatial transport coefficient k0, which has the units of permeability and is estimated in table 1 according to data provided in [6].

2.3. Radially symmetric compartment model A simpler way to augment the binding model with a spatial component is to exploit the shell-like nature of tumour cords, the geometrical property that cells are broadly arranged in concentric circles around a central blood vessel (see figures 1 and 4a). It is again assumed that variations along the vessel are negligible. We now assume that the rate of transport of drug between neighbouring shells is proportional to the shared interface area (denoted by Ai for the interface between shells i and i þ 1) and the difference in concentration across the interface. Under these assumptions, the spatial variation in the radial direction can be included by adapting equations (2.1) –(2.3) to give

d1 Vi

d 2 Vi

and

d2 Vi

dt

@ C1 ¼ Dr2 C1 þ ak1 (C2 C1 ), @t

@ C2 ¼ ak1 (C1 C2 ) d2 k2 C2 (C0 C3 ) þ d2 k2 C3 @t @ C3 d2 ¼ d2 k2 C2 (C0 C3 ) d2 k2 C3 , and @t

d2

(2:12) (2:13) (2:14)

where the factors d1 and d2 are defined as before, and D is a diffusion coefficient related to the permeability k0. The precise nature of this relationship can be derived by noting that ð þ Dr2 CdV ¼ n DrCdS (2:15) Vi @ Vi Xð n DrCdS (2:16) ¼

(2:7)

dC(i) (i) 2 ¼ ai k1 (C(i) 1 C2 ) dt dC(i) 3

d1

j[N i @ Vij

dC(i) (iþ1) 1 C(i) C(i) ¼ Ai1 k0 (C(i1) 1 1 ) þ Ai k0 (C1 1 ) dt (i) þ ai k1 (C(i) 2 C1 ),

The final model, the geometry of which is illustrated in figure 4b, is derived by disregarding the cellular structure of the tumour cord and assuming that transport of molecules occurs by isotropic Fickian diffusion in a continuum. A model for the tumour cord is then readily obtained as the system of partial differential equations (PDEs) [10,12]

(i) (i) d2 Vi k2 C(i) 2 (C0 C3 ) þ d2 Vi k2 C3

(2:8)

(i) (i) ¼ d2 Vi k2 C(i) 2 (C0 C3 ) d2 Vi k2 C3 ,

(2:9)

for i ¼ 1, . . . , n, where n is the number of shells, and ai is the cellular surface area within the ith shell. In this work, n ¼ 9 is chosen so that each shell can be identified with a layer of biological cells. The superscript corresponds to the shell number, and this index increases with distance from the blood supply (see figure 4a). The boundary conditions are imposed at the vessel wall by setting C(0) 1 ¼ Cv (t) to be a predefined PK profile

X Aij (j) D(C1 C(i) 1 ), d ij j[N

(2:17)

i

in which n represents the unit normal pointing outwards from volume Vi. Note that the direction of n varies over the surface of Vi. Hence, the relationship between transport coefficients of discrete and continuum models is D ¼ dk0 ,

(2:18)

where d is a length parameter, taken here to be the average cell diameter, 20 mm. The parameter a in equations (2.12) and (2.13) is easily derived from equations (2.7) or (2.8) and represents the ratio of the cellular surface area within a region to that region’s volume, a ¼ a/V (table 1). Given these definitions, integrating equations (2.12) – (2.14) over a cell or shell returns the equations of the multi-


The factors d1 and d2 are defined as in the cell-centre model, so that d1Vi and d2Vi are, respectively, the extracellular and intracellular volumes in the ith layer. The interface area between shells i and i þ 1 is

5


in the blood vessel and replacing k0 by kv at this interface. A no-flux boundary condition is imposed at the outer boundary of the cord by setting An ¼ 0 (equivalent to replacing k0 by zero at this interface). The volumes Vi of the shells are readily determined from the geometry: assuming a shell thickness d and a vessel radius l, the volume (again omitting the factor of the vessel length owing to the assumed uniformity along the vessel) of the ith shell is


(b)

6

L

L

Vn ··· ··· V2 V1

V

l


(a) An A2 A1 A0

l

d

Figure 4. Schematics of the two radially symmetric models: (a) the shell-like arrangement for the compartmental model and (b) the continuum model, for the tumour cord. (Online version in colour.)

dimensional cell-centre and radially symmetric compartmental models, respectively. The boundary conditions imposed are the same as for the other two models. The flux of free drug across the vessel wall is assumed to be proportional to the difference in concentration between vessel and interstitium, so [15] DrC1 n1 ¼ kv (Cv C1 jr¼l ),

A ¼ 7:46 102 l1

a ¼ 2:69 103 s1

B ¼ 2:49 103 l1

b ¼ 2:83 104 s1

C ¼ 5:52 104 l1

g ¼ 1:18 105 s1 :

(2:19)

that is, the normal flux at the vessel wall is proportional to kv, the vessel permeability and the difference in free drug concentration across the vessel wall (Cv being the drug concentration in the blood, determined by the prescribed PK profile). A no-flux condition is imposed at the outer boundary DrC1 n2 ¼ 0 at r ¼ L:

of the values given in table 2 of reference [29]. This gives (to three significant figures)

(2:20)

The respective unit normal vectors, at the vessel and the outer boundary, pointing out of the intervening tissue, are denoted by n1 and n2. A spectral method [26,27] is used in this work for the spatial discretization of equations (2.12) –(2.14). Note that, because radial symmetry is assumed, the Laplacian term in equation (2.12) is actually of the form 2 @ C1 1 @ C1 Dr2 C1 ¼ D þ : (2:21) r @r @ r2

2.5. Pharmacokinetic profiles The clinical pharmacokinetics of doxorubicin are well characterized in the literature [28]. Doxorubicin concentrations are known to decay in a tri-exponential manner following intravenous (IV) bolus or infusion and typical parameters are available in the literature (e.g. [29], in which doxorubicin is administered as an IV bolus, not infusion). The first PK profile considered here is based on the data provided by Robert et al. [29], in which Cv(t) is assumed to decay as a tri-exponential (see also [8]), i.e. D0 A at B C Cv (t) ¼ (e 1)eat þ (ebt 1)ebt þ (egt 1)egt , t a b g (2:22) for t t, where t is the infusion time, D0 is the dose and parameters A, B, C, a, b and g are estimated by taking averages

The rapid initial infusion of the drug is modelled by taking D0 A B C (1 eat ) þ (1 ebt ) þ (1 egt ) , (2:23) Cv (t) ¼ t a b g for t , t, which lifts the concentration to the appropriate value at t ¼ t. The duration of the perfusion for the total dose injected, t ¼ 180 s, was also taken from [29], and the total dose D0 ¼ 1.19827102 mmol was calculated to give an ‘area under the curve’ of ð1 AUC ; Cv (t)dt ¼ 104 mM s 2:78 mM h , (2:24) 0

which is typical of what one might find in a patient [30,31]. The AUC is an important parameter, because the area under the plasma drug concentration–time curve (AUC) is considered to reflect the actual tumour (cellular) exposure to drug after administration of a drug dose, and to correlate with toxicity—though it is more difficult to correlate with clinical efficacy [32]. Two further PK profiles, both constructed to give the same AUC, are also simulated, to investigate their influence on the drug distribution. 0

— A mono-exponential profile, Cv (t) ¼ A0 ea t with A 0 ¼ 50 mM and a 0 ¼ 0.005 s21. — A uniform (steady-state) profile, Cv (t) ¼ A00 ¼ 3:85802 102 mM up to t ¼ 72 h (and zero afterwards), representing prolonged infusion. This takes the form of a rectangular pulse.

3. Results 3.1. Model comparison The first set of numerical results are generated to address the key question ‘Does each of these models give similar results


d


Section 2.5 described three distinct PK profiles, one representing administration by an IV bolus with a tri-exponential decay (characteristic of doxorubicin), a simplified, mono-exponential, approximation to this decay profile and a uniform profile

3.3. Comparing binding affinities The final set of numerical experiments investigates the effect that changing the binding affinity of the intracellular drug, parameters k2 and k22 in our model, has on the exposure of the cells to bound drug. This attempts to address one aspect of the question ‘How does the administration schedule and cell response affect drug delivery?’ As in §3.2, the compartmental model (§2.3) was used to produce the numerical results presented here. No significant differences were seen when the same tests were carried out with the other models (data not shown). The qualitative changes caused by adjusting k2 and k22 depended most significantly on the ratio k22/k2, so only results for different values of k2 (the association rate) are shown in figure 9. All other parameters take the values shown in table 1.

7


3.2. Comparing pharmacokinetic profiles

representing infusion of the drug. All profiles had the same AUC, because this is a standard measure of cellular toxicity. Because the numerical results presented in §3.1 showed similar results for all three models of the tumour cord, the simplest—the compartmental model described in §2.3—was chosen to illustrate how the distribution of the drug is influenced by the PK profile of the supplied drug, Cv(t). The model parameters used are those given in table 1. The other models have been run with the same parameters, but the data are not shown because they contain no significant differences. Figure 6 shows the evolution of the extracellular and bound drug distributions, C1 and C3, as a function of distance from the drug source for each of the three PK profiles described in §2.5. The time variation of the concentrations of both free and bound drug in a given cell layer generally follows that of the PK profile in the vessel, though there is a time-lag which increases the further away from the vessel a cell layer is. For both exponential profiles, the concentration increases initially to a peak value ( particularly rapidly for the mono-exponential profile), then decreases monotonically, but for the uniform profile, it increases monotonically for the duration of the experiment. This is confirmed by examining the temporal variation of the concentrations at specific points in the domain. Figure 7 shows this at the centres of the first (innermost), fifth (middle) and ninth (outermost) shells of the compartmental model. Note that the extracellular drug concentrations at very early times (less than 1 h) extend far beyond the maximum value on the vertical axis of the graph for the IV bolus profiles: the true maximum values are given in table 2. These peaks become less extreme further away from the vessel. Figure 8 shows the variation of the exposure of the cells to ÐT the bound drug (AUC, given by 0 C3 dt) as a function of distance from the blood vessel, after T ¼ 24 h and T ¼ 72 h. Early in the simulations, the exponential profiles give a far higher exposure than the uniform profile, but as time progresses, the differences between the profiles reduce. However, results at later times should be interpreted carefully, because elimination of drug is not included in the model, and the only drug clearance is due to the drug returning to the vessel. This effect will be addressed in future models, through the inclusion of elimination mechanisms such as cellular metabolism, sequestration/binding to the extracellular matrix and drug efflux. Drug clearance owing to lymphatics may be considered for larger tumour volumes, though functional lymphatic vessels are not thought to be prevalent in tumours [3].


for the variation in drug concentration in the tumour cord for similar parameters?’ In order to assess this, we compare the predictions for a set of parameter values, shown in table 1, common to all three models. These simulations are carried out for the tri-exponential (IV bolus) PK profile described in §2.5. All numerical experiments were carried out on the same computational domain, comprising a single circular blood vessel of radius 16 mm at the centre of a circular cord of tumour cells of radius 200 mm. The vessel radius approximates the average radius of arterioles and venules rather than larger vessels. For the radially symmetric compartmental model, the region between the vessel and outer boundary was divided into nine shells of equal width (approximately the diameter of a biological cell). The two-dimensional cell-centre results were generated using 60 separate configurations, to assess the effect that small random variations in the distribution of the node positions and radii have on the drug distribution. Each configuration is derived from a different randomly generated set of node positions and radii (see §2.2) and contains approximately 400 nodes, the number required to fill the domain with cells of the given radius. The nature of the binding model means that the node spacing does not have to match the biological cell size, because each compartment contains both intracellular and extracellular concentrations. However, future developments may treat the cells as separate entities, immersed in interstitial fluid, so this is a useful length scale at which to investigate the model. The impact on the results of changing the spacing is similar to that of changing the resolution in an approximation to a continuum model. Increasing the spacing will effectively increase the rate of diffusion, because the model assumes that, at each time step, any drug that passes into a compartment instantaneously equidistributes its concentration throughout that compartment. We have conducted numerical experiments with different cell sizes and, although there are small quantitative differences, we have found no evidence that the qualitative behaviour might be changed. In order to visualize the two-dimensional simulation results, both mean values and standard deviations are plotted, after clustering the nodes from all 60 configurations into 20 bins, according to their distance from the blood vessel. The mean distance is plotted against the mean concentration for each of these bins, and the standard deviations of both variables are illustrated by horizontal (distance) and vertical (concentration) bars. The compartmental model is illustrated using solid dots plotted at the centres of the cylindrical shells, and the continuum model is represented by a solid line. Figure 5 shows snapshots of the concentration profiles, as a function of distance from the source of the drug, for C3 after 1, 6, 24 and 72 h, generated using the tri-exponential PK profile from §2.5 and the parameter values shown in table 1. Note that the results for C1 and C2 (in the electronic supplementary material) are almost indistinguishable. This similarity is common to all the tests we have run and consistent with very rapid transport of drug across the cell membrane.


bound intracellular drug, C3 (t = 1 h)

(a) 1.5


(b)

0.5

0

20

40 60 80 100 120 140 160 180 200 distance from vessel centre (mm) bound intracellular drug, C3 (t = 24 h)

(c) 1.5

0

20


(d)

cell-centre model continuum model compartment model concentration (mM)

40 60 80 100 120 140 160 180 200 distance from vessel centre (mm)

cell-centre model continuum model compartment model

1.0

0.5

0

20


0

20


Figure 5. Dependence of the concentrations of bound drug, C3, on distance from the drug supply at times t ¼ 1 h (a), t ¼ 6 h (b), t ¼ 24 h (c) and t ¼ 72 h (d). The solid line represents the continuum model, the large solid dots represent the compartmental model and the smaller dots with vertical and horizontal bars represent the cell-centre model. The tri-exponential (IV bolus) PK profile has been used as input. (Online version in colour.)

4. Discussion In §2, three different approaches to modelling the delivery of drug from a blood vessel to a surrounding tumour cord were presented. One of these, a discrete cell-centre model which identifies computational nodes with individual biological cells, is genuinely multi-dimensional and could be applied to more complex geometries than the one investigated here, such as those modelled using a finite-element discretization of a reaction – diffusion system in [33]. The remaining two (a discrete, compartment-based, model and a continuum, PDE-based, model) are tailored to the specific problem of drug delivery from a single vessel to a homogeneous tumour cord, assuming radial symmetry. All were built on a binding model involving extracellular drug and free and bound intracellular drug, which extends those of

[12] and [17] by allowing drug binding to be saturable and reversible. The two radially symmetric models are much simpler and therefore computationally much faster than the multidimensional model, but this leads to the question ‘Does each of these models give similar results for the variation in drug concentration in the tumour cord?’ Figure 5 shows a representative comparison of the bound drug variations of the three models for a PK profile derived from in vivo data for an IV bolus. The concentrations predicted by the three models differ by less than 15% throughout the simulation. In fact, the plotted standard deviations for the multi-dimensional cell-centre model show that the variation between the randomly generated two-dimensional configurations (particularly when there is a steep gradient in the drug concentration close to the central


concentration (mM)

1.0




8


(a) (i) 0.020

(b) (i) 2.0

0.014

1.4

0.012

1.2

0.010

1.0

0.008

0.8

0.006

0.6

0.004

0.4

0.002

0.2

0

0

(ii)

(ii) 2.0

0.020 0.018

concentration (mM)

1.8

t=1h t=6h t = 24 h t = 72 h

0.016

t=1h t=6h t = 24 h t = 72 h

1.6

0.014

1.4

0.012

1.2

0.010

1.0

0.008

0.8

0.006

0.6

0.004

0.4

0.002

0.2

0

0

(iii)

(iii) 2.0

0.020 0.018

t=1h t=6h t = 24 h t = 72 h

0.016 concentration (mM)

t=1h t=6h t = 24 h t = 72 h

1.6


concentration (mM)

0.016

1.8

1.4

0.012

1.2

0.010

1.0

0.008

0.8

0.006

0.6

0.004

0.4

0.002

0.2

0

20


t=1h t=6h t = 24 h t = 72 h

1.6

0.014

0

20


Figure 6. Dependence of the concentrations of extracellular drug, C1 (a), and bound drug, C3 (b), on distance from the drug supply for three PK profiles: triexponential/IV bolus (i), mono-exponential/simplified IV bolus (ii), uniform/infusion (iii). Each graph shows profiles at times t ¼ 1, t ¼ 6, t ¼ 24 and t ¼ 72 h. (Online version in colour.) vessel) is typically larger than the difference between the average multi-dimensional results and the radially symmetric results. This observation is supported by further simulations,


1.8

t=1h t=6h t = 24 h t = 72 h

0.018

9

some of which are shown in the electronic supplementary material. In all cases, the qualitative features of the concentration profiles are similar, whichever model is used.


(a) (i)


1.8

r ª 26 mm r ª 108 mm r ª 190 mm

0.010

1.6 1.4

0.008

1.2 1.0

0.006

0.8 0.004

0.6 r ª 26 mm r ª 108 mm r ª 190 mm

0.4

0.002

0.2 0

0 (ii)

(ii)

0.012

concentration (mM)

2.0 1.8

r ª 26 mm r ª 108 mm r ª 190 mm

0.010

1.6 1.4

0.008

1.2 1.0

0.006

0.8 0.004

0.6 r ª 26 mm r ª 108 mm r ª 190 mm

0.4

0.002

0.2

(iii)

0

0

0.012

(iii) 2.0 1.8

r ª 26 mm r ª 108 mm r ª 190 mm

0.010

concentration (mM)


concentration (mM)

10

(b) (i) 2.0

0.012

1.6 1.4

0.008

1.2 1.0

0.006

0.8 0.004

0.6 r ª 26 mm r ª 108 mm r ª 190 mm

0.4

0.002

0.2 0

10

20

30 40 time (h)

50

60

70

0

10

20

30 40 time (h)

50

60

70

Figure 7. Dependence of the concentrations of extracellular drug, C1 (a), and bound drug, C3 (b), on time for three PK profiles: tri-exponential/IV bolus (i), monoexponential/simplified IV bolus (ii), uniform/infusion (iii). Each graph shows profiles at the centres of the first, fifth and ninth of the model’s shells (distances 26, 108 and 190 mm from the centre of the vessel). (Online version in colour.)


Table 2. Maximum concentrations of free extracellular and bound intracellular drug, to 3 significant figures, and the experiment times at which they occur, to the nearest second (for three distances from the centre of the blood vessel and all three PK profiles). 5 108

9 190

tri-exponential/IV bolus maximum C1 (mM)

0.382

0.140

0.112

time of maximum C1 (h:m:s) maximum C3 (mM)

00:03:56 1.06

00:08:21 0.838

00:10:01 0.822

time of maximum C3 (h:m:s)

21:18:33

59:48:33

70:23:33

mono-exponential/IV bolus maximum C1 (mM)

1.80

0.523

0.406


00:01:15 1.93

00:04:25 0.983

00:06:14 0.872


00:41:02

01:02:28

31:45:39

0.00830

0.00681

0.00650


72:00:00 1.16

72:00:00 0.904

72:00:00 0.853


72:00:00

72:00:00

72:00:00

uniform/infusion maximum C1 (mM)

In figure 5, the two radially symmetric models agree more closely with each other than with the multi-dimensional model. This is not generally the case for all parameter sets. We note that the discrete models would be expected to converge to the continuum model asymptotically if the sizes of the ‘cells’ were allowed to tend to zero, though this is not biologically realistic, because continuum models of this type are designed for use at much larger length scales. By contrast, although the discrete models operate on realistic cell sizes, in doing so, they implicitly assume that diffusion/mixing within a cell happens instantaneously. We have conducted a range of numerical simulations for the tumour cord geometry, with different parameters, without finding any systematic differences between the results which might suggest that one model is consistently more accurate than the others. All three models give qualitatively and quantitatively similar results, and we do not yet have experimental data to enable us to assess whether one model is better or worse than another. The multi-dimensional model is, computationally, more expensive and therefore inefficient for the radially symmetric tests considered here, but we include it because it would be used for more complex geometries and heterogeneous tissue. Continuum models are very commonly used, but are based on the assumption that the differential equation is valid at every point in space. The radially symmetric, compartmentbased, model only assumes that the differential equation is valid in an integral-averaged sense (and is, in effect, a finite volume discretization of the continuum model [34], an integration of the continuum model over volumes chosen to be on the scale of the biological cells), leading to a very natural framework for simulating mass balance processes. The compartmental model is limited to radially symmetric problems (cords, multicell spheroids), but this is a common constraint imposed in computational modelling for (i) efficiency of computation and (ii) direct comparison with mathematical analysis, which is often limited to such quasi-one-dimensional geometries.

We therefore choose to use the radially symmetric, compartment-based model to investigate further the effect of varying the supplied PK profile and the behaviour of the underlying binding model, and note that the similarity of the results provides some validation of the more complex approaches. This gives us confidence that they could be used reliably in the more complex scenarios for which they are designed. A more comprehensive validation would involve comparison with measurements of drug distribution from an in vitro tumour cord model system. Simple modifications to the model geometry would also allow validation against experimental data for multicell spheroids. In both cases, knowledge of heterogeneity in the system could be readily incorporated in the multi-dimensional approach described in §2.2. Drug delivery to tumours is dependent upon a number of factors, the principal ones being the dose and schedule of administration, delivery of the drug via the blood vessels, the flux or distribution of drug through avascular tissue and consumption of drug by the cells or the extracellular matrix [3]. For example, increasing the diffusion rate, k0, tends to make the distribution of the drug more homogeneous ([35] and the electronic supplementary material), because the drug can be transported further before it is bound. However, in this paper, we focus on other factors. In §3.2, three different delivery profiles were compared, all of which provide the same overall dose. Figures 6 and 7 show significant differences between the distribution of drug in the tissue at any given point in time. When the PK profile in the blood vessel represents an IV bolus, a sharp peak in free drug concentration occurs in the first hour as the drug diffuses rapidly into the surrounding tissue and is transported into the cells. The profiles at t ¼ 1 h in figure 6 are similar in shape to the experimentally measured gradients in [23]. The concentration of free drug then drops steeply as it is bound until an approximate equilibrium is reached. This is followed by a slower decrease once the concentration of drug in the vessel drops


1 26


layer index distance from vessel centre (mm)

11


(a)


45

35

90 80 70 60

25

50

20

40

15

30

10

20

5

10 20


0

tri-exponential mono-exponential uniform 20


Ð Figure 8. Dependence of the exposure to bound drug ( C3 dt) on distance from the supply of drug at t ¼ 24 h (a) and t ¼ 72 h (b) for three PK profiles: tri-exponential/IV bolus, mono-exponential/simplified IV bolus, uniform/infusion. (Online version in colour.)


250

k2 = 0.4 M−1 s−1

200

exposure (mM h)

k2 = 4 M−1 s−1 k2 = 40 M−1 s−1

150

100

50

0

20


Ð Figure 9. Dependence of the exposure to bound drug ( C3 dt) on distance from the supply of drug at t ¼ 72 h for three different binding affinities where the association rate constant only varies over the range: k2 ¼ 4 1027, k2 ¼ 4 1026 and k2 ¼ 4 1025 mM21 s21. The tri-exponential/IV bolus PK profile was used as input. (Online version in colour.)

below that of the free drug in the tissue, and the net flux is of drug returning to the vessel. Because the binding process acts more slowly than the initial influx of drug, the concentration of bound drug changes less rapidly and only cells close to the vessel experience an initial peak in concentration (most extreme for the mono-exponential PK profile). As a consequence, early in the simulation, when the concentration in the vessel is high and changing rapidly, the drug concentration close to the vessel is much higher than in cells further away, because drug has not had sufficient time to reach and bind to

12

cells at a relatively large distance from the vessel. Later, the drug concentration in the vessel is still decreasing, but slowly enough, relative to the rates at which it is transported through the tissue and binds to the cells, for the drug distribution in the tissue to remain uniform throughout. When the PK profile represents infusion over a longer period, the concentrations in the tissue steadily increase during the period of infusion, after an initial rapid increase in free drug concentrations. The spatial distribution of the concentration is quite even, though slightly higher close to the vessel, because the concentration of supplied drug is not varying in time. Given the large differences in concentration profiles, it might be expected that the exposure to bound drug also depends critically on the supplied PK profile. The results shown in figure 8 suggest that, for the model of binding proposed in §2.1, all of the PK profiles give a similar spatial distribution of exposure to bound drug. After 72 h, infusion gives a significantly lower exposure than IV bolus (40–50% lower than the mono-exponential profile), but this difference is less significant than after 24 h and continues to reduce over longer time scales. However, our binding model contains no explicit elimination or decay term: drug can only leave the system through free drug returning to the vessel when the concentration in the vessel is below that in the adjacent tissue. A more sophisticated binding model could account for this and include a representation of the cell cycle, which will have an influence over longer time scales. These would need to be included before investigating the dependence of exposure on PK profiles over longer time scales. In order to assess how the binding affects the delivery of the drug, a final set of numerical simulations investigated the influence of binding affinity. Figure 9 shows that the exposure of the cells to bound drug depends strongly on the ratio of the intracellular drug association and dissociation rates (b ¼ k22/k2). For the original parameter set (table 1) b 16 mM, for which the exposure of the cells to bound drug was fairly uniformly distributed between the blood vessel and the


30

0



tri-exponential mono-exponential uniform

40

exposure (mM h)

(b)


Funding statement. This work was funded by Yorkshire Cancer Research Project grant no. B208, with additional support from the Leeds Cancer Research UK Centre.

References 1.

2.

3.

4.

5.

6.

7.

Vaupel P. 2008 Hypoxia and aggressive tumor phenotype: implications for therapy and prognosis. Oncologist 13(Suppl. 3), 21 –26. (doi:10.1634/ theoncologist.13-S3-21) Denekamp J, Dasu A, Waites A. 1998 Vasculature and microenvironmental gradients: the missing links in novel approaches to cancer therapy? Adv. Enzyme Regul. 38, 281–299. (doi:10.1016/S00652571(97)00015-0) Minchinton AI, Tannock IF. 2006 Drug penetration in solid tumours. Nat. Rev. Cancer 6, 583 –592. (doi:10.1038/nrc1893) Jain RK. 1999 Transport of molecules, particles, and cells in solid tumors. Annu. Rev. Biomed. Eng. 1, 241–263. (doi:10.1146/annurev.bioeng.1.1.241) Pitt-Francis J et al. 2009 Chaste: a test-driven approach to software development for biological modelling. Comput. Phys. Commun. 180, 2452–2471. (doi:10.1016/j.cpc.2009.07.019) Evans CJ, Phillips RM, Jones PF, Loadman PM, Sleeman BD, Twelves CJ, Smye SW. 2009 A mathematical model of doxorubicin penetration through multicellular layers. J. Theor. Biol. 257, 598–608. (doi:10.1016/j.jtbi.2008.11.031) Evans ND, Errington RJ, Shelley M, Feeney GP, Chapman MJ, Godfrey KR, Smith PJ, Chappell MJ. 2004 A mathematical model for the in vitro kinetics of the anti-cancer agent topotecan. Math.

8.

9.

10.

11.

12.

13.

Biosci. 189, 185–217. (doi:10.1016/j.mbs. 2004.01.007) El-Kareh AW, Secomb TW. 2000 A mathematical model for comparison of bolus injection, continuous infusion, and liposomal delivery of doxorubicin to tumor cells. Neoplasia 2, 325–338. (doi:10.1038/sj. neo.7900096) Lankelma J, Fernańdez Luque R, Dekker H, Schinkel W, Pinedo HM. 2000 A mathematical model of drug transport in human breast cancer. Microvasc. Res. 59, 149–161. (doi:10.1006/mvre.1999.2218) Eikenberry S. 2009 A tumor cord model for doxorubicin delivery and dose optimization in solid tumors. Theor. Biol. Med. Model. 6, 16. (doi:10. 1186/1742-4682-6-16) Ribba B, Marron K, Agur Z, Alarcon T, Maini P. 2005 A mathematical model of doxorubicin treatment efficacy for non-Hodgkin’s lymphoma: investigation of the current protocol through theoretical modelling results. Bull. Math. Biol. 67, 79 –99. (doi:10.1016/j.bulm.2004.06.007) Jackson TL. 2003 Intracellular accumulation and mechanism of action of doxorubicin in a spatiotemporal tumor model. J. Theor. Biol. 220, 201 –213. (doi:10.1006/jtbi.2003.3156) Chapman S, Shipley R. 2010 Multiscale modeling of fluid and drug transport in tumors. Bull. Math Biol. 72, 1464–1491. (doi:10.1007/s11538-010-9504-9)

14. Swanson KR, Alvord Jr EC, Murray JD. 2002 Quantifying efficacy of chemotherapy of brain tumors with homogeneous and heterogeneous drug delivery. Acta Biotheor. 50, 223 –237. (doi:10.1023/ A:1022644031905) 15. Alarcon T, Byrne HM, Maini PK. 2005 A multiple scale model for tumor growth. Multiscale Model. Simul. 3, 440–475. (doi:10.1137/040603760) 16. Thurber GM, Weissleder R. 2011 A systems approach for tumor pharmacokinetics. PLoS ONE 6, e24696. (doi:10.1371/journal.pone.0024696) 17. Dordal MS, Ho AC, Jackson-Stone M, Fu YF, Goolsby CL, Winter JN. 1995 Flow cytometric assessment of the cellular pharmacokinetics of fluorescent drugs. Cytometry 20, 307–314. (doi:10.1002/cyto.990200406) 18. El-Kareh AW, Secomb TW. 2003 A mathematical model for cisplatin cellular pharmacodynamics. Neoplasia 5, 161 –169. 19. El-Kareh AW, Secomb TW. 2005 Two-mechanism peak concentration model for cellular pharmacodynamics of doxorubicin. Neoplasia 7, 705–713. (doi:10.1593/neo.05118) 20. Nicholls GM, Clark BJ, Brown JE. 1992 Solid phase extraction and optimised separation of doxorubicin, epirubicin and their metabolites using reversedphase high performance liquid-chromatography. J. Pharm. Biomed. 10, 949 –957. (doi:10.1016/ 0731-7085(91)80104-H)

13


sufficiently general to be applied to other drugs and other cell lines if the appropriate data are available. This may require the design of new binding models. There is no sufficient evidence to suggest that any of the three models of drug transport proposed in this paper is better than the others, so the simplest was chosen to investigate the influence of the delivery profile and the cell biology. However, the comparison has validated the more flexible multidimensional model, which therefore provides a framework that can be used to gain insights into progressively more complex situations in which the influence of the characteristics of the tumour microenvironment on the PK delivery of the drug and the effects of spatial heterogeneity can be investigated. It is important to emphasize that all models are approximations to reality and the models described here are clearly significant simplifications of complex biology and geometry. The value of models in biology and medicine lies in their role in the iterative development of a quantitative, logical, predictive framework, placing them at the heart of ‘modelbuilding’; the need to write down equations describing the biological mechanisms demands assumptions and yields predictions which can be tested and measurements which have to be made. This, in turn, leads to improved models and initiates a further cycle of experimentation and model building. We also believe that these models have the ability to demonstrate, at least semi-quantitatively, the relative efficacy of some aspects of therapeutic protocols.


outer boundary of the domain (200 mm from the centre of the vessel), decreasing only slightly with distance from the vessel. When b is increased—the net affinity for the drug is reduced—the spatial distribution remains uniform, but the exposure of each layer of cells is reduced. When b is decreased, the exposure of the cells close to the vessel increases, but the exposure further away from the vessel actually starts to decrease. If k2 is increased further than shown in figure 9, then this behaviour becomes more pronounced, to the point where almost no drug gets beyond 100 mm from the supply, because it is all consumed by the cells close to the vessel. This ‘binding site barrier’ [35,36] is reduced in tissue which allows more rapid interstitial drug diffusion (see numerical results in the electronic supplementary material). This suggests that, for this model of binding, there is an optimal binding affinity which allows optimal exposure to doxorubicin: if the binding is too weak, then none of the cells gains enough exposure; if it is too strong, then the cells distant from the vessel are pharmacokinetically resistant, because the closer cells consume the drug. If true, then this might have implications for the use of additional agents which can alter the binding behaviour of doxorubicin, e.g. by occupying binding sites, to adjust drug penetration. It also suggests that variations in the tumour microenvironment could influence the effectiveness of the drug. The model presented in this paper has been tailored to simulate a tumour cord geometry, and the parameters have been estimated based on the binding of doxorubicin to colorectal adenocarcinoma cells (DLD-1). However, it is


26. 27. 28.

30.

31.

32.

33.

34.

35.

36.

doxorubicin in patients with small cell lung cancer. Clin. Pharmacol. Ther. 53, 555–561. (doi:10.1038/ clpt.1993.69) Phillips RM, Bibby MC, Double JA. 1990 A critical appraisal of the predictive value of in vitro chemosensitivity assays. J. Natl. Cancer Inst. 82, 1457– 1468. (doi:10.1093/jnci/82.18.1457) Mo¨nnich D, Troost EGC, Kaanders JHAM, Oyen WJG, Alber M, Thorwarth D. 2011 Modelling and simulation of [18F]fluoromisonidazole dynamics based on histology-derived microvessel maps. Phys. Med. Biol. 56, 2045–2057. (doi:10.1088/00319155/56/7/009) Tannehill JC, Anderson DA, Pletcher RH. 1997 Computational fluid mechanics and heat transfer. Washington, DC: Taylor & Francis. Rippley RK, Stokes CL. 1995 Effects of cellular pharmacology on drug distribution in tissues. Biophys. J. 69, 825–839. (doi:10.1016/S00063495(95)79956-8) Fujimori K, Covell DG, Fletcher JE, Weinstein JN. 1989 Modeling analysis of the global and microscopic distribution of immunoglobin G, F(ab’)2, and Fab in tumors. Cancer Res. 49, 5656– 5663.

14


29.

xenografts by dynamic contrast-enhanced MRI. J. Magn. Reson. Imaging 26, 1033–1042. (doi:10. 1002/jmri.21110) Boyd JB. 2011 Chebyshev and Fourier spectral methods. New York, NY: Dover. Trefethen LN. 2000 Spectral methods in Matlab. Philadelphia, PA: SIAM. Speth PA, van Hoesel QG, Haanen C. 1988 Clinical pharmacokinetics of doxorubicin. Clin. Pharmacokinet. 15, 15 –31. (doi:10.2165/00003088198815010-00002) Robert J, Illiadis A, Hoerni B, Cano J-P, Durand M, Lagarde C. 1982 Pharmacokinetics of Adriamycin in patients with breast cancer: correlation between pharmacokinetic parameters and clinical short-term response. Eur. J. Cancer Clin. Oncol. 18, 739–745. (doi:10.1016/02775379(82)90072-4) Greene RF, Collins JM, Jenkins JF, Speyer JL, Myers CE. 1983 Plasma pharmacokinetics of Adriamycin and Adriamycinol: implications for the design of in vitro experiments and tretment protocols. Cancer Res. 43, 3417–3421. Piscitelli SC, Rodvold KA, Rushing DA, Tewksbury DA. 1993 Pharmacokinetics and pharmacodynamics of


21. Harada H, Xie XJ, Itasaka S, Zeng LH, Zhu YX, Morinibu A, Shinomiya K, Hiraoka M. 2008 Diameter of tumor blood vessels is a good parameter to estimate HIF-1-active regions in solid tumours. Biochem. Biophys. Res. Commun. 373, 533– 538. (doi:10.1016/j.bbrc.2008.06.062) 22. Moore JV, Hasleton PS, Buckley CH. 1985 Tumor cords in 52 human bronchial and cervical squamous-cell carcinomas: inferences for their cellular kinetics and radiobiology. Br. J. Cancer 51, 407–413. (doi:10.1038/bjc.1985.55) 23. Primeau AJ, Rendon A, Hedley D, Lilge L, Tannock IF. 2005 The distribution of the anticancer drug doxorubicin in relation to blood vessels in solid tumors. Clin. Cancer Res. 11, 8782–8788. (doi:10. 1158/1078-0432.CCR-05-1664) 24. Lee C, Kim JS, Waldman T. 2004 PTEN gene targeting reveals a radiation-induced size checkpoint in human cancer cells. Cancer Res. 64, 6906 –6914. (doi:10.1158/0008-5472.CAN04-1767) 25. Vestvik IK, Egeland TAM, Gaustad J-V, Mathiesen B, Rofstad EK. 2007 Assessment of microvascular density, extracellular volume fraction, and radiobiological hypoxia in human melanoma

Cancer evolution: mathematical models and computational inference.

Computational and Mathematical Methods in Cardiovascular Diseases.

Mathematical methods and models in system biomedicine.

Antiretroviral treatment, preexposure prophylaxis, and drug resistance in sub-Saharan Africa: a consensus among mathematical models.

Mathematical Models of College Myopia.

Presurgical diagnosis of adnexal tumours using mathematical models and scoring systems: a systematic review and meta-analysis.

Mathematical Models of Tuberculosis Reactivation and Relapse.

Mathematical models of gene amplification with applications to cellular drug resistance and tumorigenicity.

Mathematical models of breast and ovarian cancers.

Computational models of planning.

Mathematical and Computational Modeling in Complex Biological Systems.

Mathematical and computational modeling in biology at multiple scales.

Mathematical and computational modelling of skin biophysics: a review.

Mathematical models for skin toxicology.

HIV treatment, preexposure prophylaxis, and drug resistance: reconciling conflicting predictions from mathematical models.

Mathematical models for translational and clinical oncology.

Computational analysis of drug transport in tumor microenvironment as a critical compartment for nanotherapeutic pharmacokinetics.

Reproducibility in Computational Neuroscience Models and Simulations.

Mathematical models in drug delivery: how modeling has shaped the way we design new drug delivery systems.

Computational models of heart disease.

Antibody transport and internalization into tumours.

Mathematical models of repopulation and reoxygenation in radiotherapy.

Mathematical models of human african trypanosomiasis epidemiology.

Special issue on Mathematical Models of Inflammation.