University of Huddersfield Repository Drug delivery in a tumour cord model: a computational simulation

The tumour vasculature and microenvironment is complex and heterogeneous, contributing to reduced delivery of cancer drugs to the tumour. We have developed an in silico model of drug transport in a tumour cord to explore the effect of different drug regimes over a 72h period and how changes in pharmacokinetic parameters affect tumour exposure to the cytotoxic drug doxorubicin. We used the model to describe the radial and axial distribution of drug in the tumour cord as a function of changes in the transport rate across the cell membrane, blood vessel and intercellular permeability, ﬂow rate, and the binding and unbinding ratio of drug within the cancer cells. We explored how changes in these parameters may affect cellular exposure to drug. The model demonstrates the extent to which distance from the supplying vessel inﬂuences drug levels and the effect of dosing schedule in relation to saturation of drug-binding sites. It also shows the likely impact on drug distribution of the aberrant vasculature seen within tumours. The model can be adapted for other drugs and extended to include other parameters. The analysis conﬁrms that computational models can play a role in understanding novel cancer therapies to optimize drug administration and delivery.

The tumour vasculature and microenvironment is complex and heterogeneous, contributing to reduced delivery of cancer drugs to the tumour. We have developed an in silico model of drug transport in a tumour cord to explore the effect of different drug regimes over a 72 h period and how changes in pharmacokinetic parameters affect tumour exposure to the cytotoxic drug doxorubicin. We used the model to describe the radial and axial distribution of drug in the tumour cord as a function of changes in the transport rate across the cell membrane, blood vessel and intercellular permeability, flow rate, and the binding and unbinding ratio of drug within the cancer cells. We explored how changes in these parameters may affect cellular exposure to drug. The model demonstrates the extent to which distance from the supplying vessel influences drug levels and the effect of dosing schedule in relation to saturation of drug-binding sites. It also shows the likely impact on drug distribution of the aberrant vasculature seen within tumours. The model can be adapted for other drugs and extended to include other parameters. The analysis confirms that computational models can play a role in understanding novel cancer therapies to optimize drug administration and delivery. In §2, we first outline the underlying binding model, which is parameterized by experimental data for doxorubicin, and then describe the full multidimensional tumour cord model and its discretization. Section 3 contains a summary of the three PK profiles we compare. In §4, a thorough comparison of the predicted effects owing to each of the PK profiles is presented, and the consequences of varying model parameters are investigated. These results are then summarized and discussed in § §5 and 6.

Models
The work in this paper builds on the one-dimensional, radially symmetric, compartment model we previously proposed [9]. Three distinct approximations of the radial variation in a longitudinally uniform tumour cord geometry were compared and we concluded that the quantitative and qualitative differences between the simulated results were small enough, relative to the uncertainty in the modelling process and in the experimental data that were used to parameterize the underlying binding processes, to select the simplest and fastest (referred to in [9] as the radially symmetric compartment model) as the most appropriate for predictive simulations.
In this paper, we retain the same binding model, which is recapitulated in §2.1 for completeness, and recall that no explicit elimination or decay of drug is included beyond clearance via the central blood vessel which supplies the drug initially. This was assumed owing to the lack of functional lymphatics within tumour tissue, and the lack of data with which we might parameterize any additional generic clearance within living tissue. Because this omission is likely to render the model invalid for longer time periods, we limit our simulations to a maximum of 72 h.
In §2.2, we generalize the radially symmetric compartment model of spatial variation of the drug concentration in a cylindrical tumour cord, centred on a single vessel supplying the drug [9], to allow for variation along the direction of the supplying vessel. The radial symmetry is retained, but the drug transport through the tissue is modelled in a two-dimensional plane, in the radial and axial coordinates. The advective transport of drug within this vessel is also modelled.

Binding model
The interaction of the chemotherapeutic agent with the microenvironment of cells is restricted to drug binding only, and described by a three-compartment model, composed of extracellular space (volume V 1 ) with a concentration C 1 of free drug and intracellular space (volume V 2 ) with concentrations C 2 and C 3 corresponding to free and bound drug (where we understand the term bound to include both DNAintercalated drug and drug bound to the cell in other ways). Binding is described by a simple kinetic model with drug binding reversibly to sites within the cell, which represents the principal form of binding for doxorubicin in living cells [11].
Applying the principle of mass action leads to three coupled ordinary differential equations which describe the system [9]: and in which k 1 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 cells), k 2 and k −2 are the drug binding and unbinding rates, respectively, and C 0 is the concentration of binding sites within the cell. This model is illustrated by the schematic in figure 1. Values for the kinetic rate constants for the binding process, given in table 1, have been derived from a bespoke experimental binding assay, outlined in [9].
Remark. The kinetic rate parameters were derived under the assumption that the binding was reversible. In the limit of irreversible binding, i.e. the limit as β = k −2 /k 2 → 0 in equations (2.1)-(2.3), the mass-conserving steady-state solution is given by  Figure 1. A three-compartment model of drug distribution in tissue. C 1 represents extracellular drug concentration, C 2 is free intracellular drug concentration and C 3 is bound intracellular drug concentration.  [8,9], and 'histology' indicates estimation from histological tissue images, such as those illustrated at www.virtualpathology.leeds.ac.uk. The parameter k 0 has been estimated based on the value of r 1 (transport rate between cell layers) in the multilayer model of [8]      in which δ = V 1 /V 2 and C 1 (0) is the initial concentration of extracellular drug. In other words, the drug continues to bind until either all of the free drug is used up or the binding sites become saturated. These equations predict a piecewise linear dependence of the steady-state concentrations on C 1 (0), which is clearly not mirrored by the experimental data used to fit the model parameters in [9]. Hence, the data support the assumption that the binding is reversible when using this model to represent the dynamics.

Spatial model
Previously [9], the binding model was augmented with a spatial component by exploiting the shell-like nature of tumour cords, the geometric property that cells are broadly arranged in concentric cylindrical shells around a central blood vessel. In this work, a generalization of that one-dimensional model will be considered, which retains the cylindrical symmetry but allows variation in the axial direction as well as the radial direction. This enables the transport of drug along the vessel to be modelled as well as transport through the tissue. In our two-dimensional model, we distinguish between the coordinate directions in the definition of the geometry by using subscripts r and z to represent, respectively, the radial and axial coordinates.  In tissue, the spatial variation of the discrete drug concentrations are now given by and for indices i = 1, . . . , n and j = 1, . . . , m, where n is the number of shells in the radial direction, m is the number of discs in the axial direction and a i is the cellular surface area within any compartment in the ith shell: the volumes of these compartments are independent of the axial index j. The superscripts correspond, therefore, to the shell and disc indices of the discrete compartments, as illustrated by compartment i,j in figure 2. In this work, n = 9 and m = 25 are chosen with d = 20 µm so that each shell (and disc) can be approximately identified with a layer of biological cells. The geometry of the compartments is defined by in which d r = d z = d gives the dimensions of each compartment. Choosing d r = d z means that the diffusive transport is isotropic when the permeability k 0 is the same in both coordinate directions. Axial uniformity is assumed in the compartment geometry, so none of the local volumes or areas depends on the axial index, j. Note that the half-indices represent interfaces between shells.
At the vessel wall boundary the term A 1/2 is the concentration of drug in the adjacent part of the supplying vessel. At the outer boundary in the radial direction, the term A n+1/2 r k 0 (C n+1,j 1 − C n,j 1 ) is replaced with zero (also for j = 1, . . . , m). A no-flux boundary condition is also applied at the boundaries in the axial direction, where the terms A i are both replaced with zero along the whole boundary, i.e. for i = 1, . . . , n.
The drug concentration in the one-dimensional representation of the vessel, denoted here by C v , is modelled using advective transport. This is coupled with the tissue by a term that allows for passage of drug through the vessel wall into the adjacent layer of cells, and we do not model the plasma binding as a separate compartment. The variation along the vessel is, therefore, governed by where the index j corresponds to the axial index in the adjacent tissue compartment and velocity λ drives flow along the vessel and no boundary condition is required at the outflow end of the vessel segment. The geometry is defined by where d z = d is the width of a disc in the axial direction. The full geometry of the two-dimensional approximation is illustrated in figure 2, in which the vessel is represented by the central cylinder, and the values for the geometric and transport parameters are given in table 1.

Pharmacokinetic profiles
The clinical PKs of doxorubicin are well characterized in the literature [20]. Doxorubicin concentrations decay in a tri-exponential manner following intravenous (IV) bolus or infusion and typical parameters are available in the literature [21], in which doxorubicin is administered as an IV bolus. Clearance of the drug is modelled via the time dependence of this concentration profile. The first PK profile considered here, subsequently denoted as PK1, is based on the data [21], in which C v (t) is assumed to decay as a tri-exponential (see also [22]) after a short infusion, i.e.: where τ is the infusion time, D 0 is the dose and parameters A, B, C, A , B and C are estimated by taking averages of the values given in table 2 of [21]. This gives (to three significant figures) The duration of the perfusion for the total dose injected, τ = 180 s, was also taken from Robert et al. [21] and the total dose D 0 = 1.19827 × 10 2 µmol was calculated to give an 'area under curve' (AUC) of which is typical of what one might find in a patient [23,24]. The AUC is considered to reflect the actual tumour (i.e. cellular) exposure to drug, and to correlate with toxicity, i.e. to a lesser extent with clinical efficacy [25]. Two further PK profiles, both constructed to achieve the same AUC, are also considered in this work: -a profile, denoted by PK2, in which there are three short infusions (to mimic the repeat drug administration suggested in [10]), each of one-third the dose of C PK1 v (t), delivered at 24 h intervals, i.e.: -a uniform profile, denoted by PK3, representing prolonged exposure at constant concentration over a period of 3 h, given by These three profiles and the accumulation of their AUCs over a period of 72 h are illustrated in the electronic supplementary material, figures S1 and S2.

Results
Numerical experiments were carried out to investigate the effects of changing the PK profile of the supplied drug (to imitate different modes of delivery) and the model parameters.   Figure 3 illustrates the spatial variation of the free extracellular (C 1 ) and bound intracellular (C 3 ) exposure profiles at 72 h for two sample simulations. The top plots show the results obtained for the parameters given in table 1; exposure decreases away from the vessel, but changes in the axial direction are hardly visible over the 500 µm length of the simulated domain. The values of the tissue binding and transport parameters, k 0 , k 1 , k 2 , k −2 and C 0 , have been estimated for specific cancer cell lines [8,9], but the vessel-related parameter values, k v , l and λ, are likely to be highly variable in the leaky, chaotic, vasculature that is characteristic of cancerous tissue. Figure 3c,d illustrates how this might affect the spatial variation of the exposure by simulating a narrower vessel, with a less obstructive wall and slower flow; variation in the axial direction is now visible for this particular choice of modified vasculature, suggesting a higher risk of pharmacokinetic resistance owing to the non-uniformity of the distribution of drug.

Varying vessel properties
A more thorough examination of the effects of changing parameters is described in the next section. Each model parameter was varied independently, with all remaining parameters retaining the values given in table 1. All simulations were run to t max = 72 h and the output predicted by the model is illustrated using total exposure to bound drug at specific points in the tissue, We treat this as a determinant of the effect of the drug, in that increasing the exposure to bound drug of a cell should signify an increased likelihood of efficacy (and/or toxicity).  in §3. The model output leads to the following observations which are characteristic of all simulation results obtained.

Varying model parameters
-The exposure of tissue to bound drug far from the supply is invariably lower than that close to the supply. -The exposure due to a single short infusion (PK1) is lower than that due to a uniform profile (PK3) with the same AUC, typically by between 5 and 20%. The profile with three short infusions (PK2) gives even lower exposures. We emphasize here that the simulations are run to 72 h and that neither of the short infusion profiles delivers its full AUC within this timescale (as illustrated in the electronic supplementary material, figure S2). When the drug concentration throughout the vessel is held at the source concentration, C v (t), instead of modelling its escape into the surrounding tissue using equation (2.9), there is no variation in the axial direction and the concentration in the tissue depends only on the radial coordinate, indexed by i. This situation is closely approximated when simulations are run with the parameters given in table 1, as illustrated in figure 3a,b. In this case, as long as the prescribed, time-dependent profile satisfies C v (t) → 0 as t → 0, it can be shown (see appendix A) that in which i represents the radial index of the computational shell. In other words, the AUCs for the free drug compartments are not only independent of the shape of the supplied PK profile but also equal to the AUC of the supplied profile. For the bound drug so its AUC does vary in space and depends on the shape of the supplied profile. However, when the concentration of bound drug is far from saturating the available binding sites, i.e. when C i 3 C 0 ∀i, it follows that the AUC for the bound drug compartment only depends weakly on the shape of the supplied PK profile, because -Close to the vessel, drug exposure increases as β = k −2 /k 2 decreases, i.e. with stronger binding and/or weaker unbinding. However, at greater distances from the vessel there is an optimal value of β, below which the exposure of the tissue to bound drug decreases. This is because, when binding is relatively strong, the drug is captured by the cells close to the vessel and not transported to the more distant cells.
The value of β given in table 1 (indicated by the vertical dashed line in the figure) gives a very low exposure, relative to what might be achieved close to the vessel with stronger binding. However, it also gives a fairly uniform distribution through the tissue (figure 3b), so all regions of the tissue receive similar amounts of drug. In fact, the exposure furthest from the supply is close to the maximum achievable with the other parameters fixed. This is more clearly visible in the magnified plot shown in figure 5, and in figure 6b, where a ridge indicates the values of the binding rate k 2 , and the unbinding rate k −2 , for which exposure at a distance from the supply is maximized. Figure 6 also shows that (except for very weak unbinding) the exposure depends on the ratio, β, of binding to unbinding rates and not the separate values of k 2 and k −2 . This observation is supported by the mathematical analysis in appendix A, which shows that, in the quasi-onedimensional case where there is no variation in the axial direction This does not necessarily hold for small values of β in the fully two-dimensional case because strong binding causes the free drug concentration in the vessel to vary significantly in the axial direction. -Exposure increases with k 1 , the rate of transport across the cell membrane. However, the value for k 1 given in table 1 (indicated by the vertical dashed line) is in a region of parameter space where the exposure is insensitive to its value. This confirms the observation made when fitting the rate parameters to experimental data in [9], in which the value of k 1 was chosen to be close to the minimum that could be taken without significantly changing the fit to the experimental assay data. -When the diffusive transport rate in the tissue (represented here by the permeability k 0 ) is varied, the effect on exposure depends on the distance from the supplying vessel. Close to the vessel the exposure actually decreases as the diffusion rate increases because, when transport is obstructed, the drug accumulates close to the vessel. Since it inhibits drug from moving away from the vessel, a low diffusion rate corresponds to low exposure further from the vessel. When the  diffusion rate is high enough, the drug is uniformly distributed through the modelled region of tissue, so the exposure is the same throughout. -Increasing k v , the permeability of the vessel, increases the exposure throughout the tissue. Very similar behaviour is seen when l, the radius of the supplying vessel, or λ, the velocity of flow in the vessel, is increased.

Varying dose
In this paper, we choose predominantly to present the model predictions in terms of exposure to bound drug, treating it as an indicator of treatment efficacy. Having investigated its dependence on model parameters, we now analyse the relationship of drug exposure with the dose of drug, D 0 in equation (  of §3, for each of the PK profiles presented in that section. Figure 7 shows the dependence of exposure on dose at positions close to and far from the drug supply. -For low doses, all three PK profiles give very similar exposures, even though the time variation of C 3 is very different for each profile. The electronic supplementary material, figure S4, provides an example of how C 3 can vary in time at different distances from the supply. -For high doses, the uniform profile (PK3), representing prolonged subjection to a constant concentration, gives significantly lower exposure than the other two. The cause is illustrated in figure 8. For high doses the binding sites rapidly become saturated, i.e. C 3 ≈ C 0 : in all cases the uniform profile (PK3) drops to zero after 3 h, so the drug starts to return to the vessel, but the two short infusion profiles (PK1 and PK2) continue to supply drug at levels which maintain saturation for the full 72 h simulation time. It is clear from equation (4.2) that, even when variation is only allowed in the radial direction, the shape of the supply profile might have a significant effect on exposure when saturation is approached, though we acknowledge that profile PK3 has been chosen to represent a specific, exaggerated circumstance which might be created in vitro but is not physiologically realistic. -The exposure to bound drug is fairly uniform throughout the computational domain for all doses when the parameter values from table 1 are used (e.g. figure 3a,b), so very little difference is seen between the exposures close to and far from the supply in figure 7a,b. However, when the values of the binding and unbinding rates, k 2 and k −2 , respectively, are replaced with those corresponding to the far right corner of the domain in figure 6 (an extreme situation, chosen to  illustrate the consequences of a microenvironment which supports strong spatial heterogeneity in drug distribution), significant differences appear between different parts of the tissue, which are visible in figure 7c,d. The cells furthest from the drug supply now require a dose two orders of magnitude higher than cells adjacent to the supply if they are to receive the same exposure. It can also be seen in figure 7 that the exposures of all three PK profiles exhibit similar dependence on dose when the stronger binding is used. In this case, saturation at high doses is maintained for longer than with weaker binding when the uniform profile (PK3) is administered, so the exposures become close to those of both short infusion profiles. For the specific parameter values chosen here, the exposure due to the profile with three short infusions (PK2) is lowest, but this does not always hold, e.g. the profiles shown in figure 7a,b indicate that PK2 gives the greatest exposure of the three for high doses when the original set of parameters (table 1) is used.
Remark. Instead of using exposure to bound drug as an indicator of the response to drug, it could be assumed that the damage done to a cell by the drug depends monotonically on its exposure to bound drug; an individual cell might survive if the damage (exposure) is below a particular threshold value or die if it exceeds that threshold. In models without spatial heterogeneity, the exposure could be converted into a survival fraction using, for example, a Hill equation [26]; we can imitate this by applying a threshold to the exposure to bound drug in each computational compartment (the sizes of which have been chosen to match the size of a biological cell). Applying the same threshold to every compartment is analogous to using a Heaviside function in place of the Hill equation-it represents the limit as the Hill exponent tends to infinity (and the higher the Hill exponent, the lower is the phenotypic heterogeneity of the population). In practice, we do not have the data to allow us to estimate a value for this threshold (and certainly not its variability between cells), so we provide figure 9 for illustrative purposes only. It shows the dependence of survival fraction, in the sense defined above, on dose for the two sets of parameters used to produce figure 7. A low survival threshold was chosen, for which the profile with three short infusions (PK2) clearly gives lower efficacy, which corresponds to the higher doses required to achieve this exposure threshold value. However, it can be seen from figure 7 that the order of effectiveness of the three PK profiles depends on both the threshold value and the choice of model parameters. An extreme example of this would be applying a threshold of 1.5 × 10 5 µM h when the parameters from table 1 are used: in this case three short infusions is the most effective mode of delivery and prolonged exposure at constant concentration kills no cells at all. Allowing the threshold to vary between cells would tend to smooth out the transition regions visible in figure 9.

Varying distance from source
To assess PK resistance for regions where the vasculature is more sparse, we consider a modified geometry in which an avascular tumour spheroid is surrounded by vascularized tissue in which the drug concentration is that of the specified PK profile. This geometry and the numerical approximation are described in detail in the electronic supplementary material, S1.3. Figure 10 shows the dependence of the exposure to bound drug on the radius of the spheroid, sampling at both the edge and the centre of the spheroid. It can be seen that exposure close to the surface of the spheroid has very little dependence on the size of the spheroid. However, the exposure to bound drug of the centre of the spheroid is strongly influenced by its distance from the supply: for the parameter values given in table 1, the exposure at the centre of a 2 mm radius spheroid has dropped almost to zero. This region would typically also be nutrient-deficient and possibly necrotic.

Discussion
A model of drug transport in a tumour cord has been developed to explore the effect of different drug delivery schedules and drug-binding parameters on tumour exposure to doxorubicin over a period of up to 72 h. In particular, we have used the model to assess the benefits of sequential administration of drug doses and at steady state. In common with all computational approaches, the model makes a number of simplifications, including neglecting changes in drug clearance beyond that reflected by the concentration-time curve of drug levels in the central blood vessel. Model parameters have been determined a priori from a combination of in vitro experimental findings and the literature. Importantly, all model parameters (apart from the transport rate across the cell membrane, for which only a minimum value could be estimated from the experimental assay data) are within or close to the regions of parameter space in which the exposure is sensitive to their values (indicated by the proximity of the vertical dashed lines to the regions of the graphs in which exposure is changing rapidly in figures [4][5][6]. Likewise, small modifications of some of the parameters, representing tumour vascularization, for example, can potentially have a significant impact on drug delivery and, therefore, certain drug combinations should be administered carefully. Subject to the caveats above, the model generates several relevant findings. Exposure to bound drug is lower further away from the supply, as illustrated in figure 3. This is intuitively obvious, but mathematical analysis of the model (see appendix A) confirms further that exposure to bound drug is similar for all three PK profiles, except when saturation of binding sites is approached; at that point the uniform profile (prolonged subjection to a constant concentration) gives significantly lower exposures. This is likely to be true not only for tumours but also for healthy tissues in relation to toxicity. Interestingly, some clinical studies have suggested that continuous infusion of doxorubicin is associated with less cardiac toxicity than bolus administration [27,28]. This has been attributed to lower peak drug concentrations, but the model suggests that overall exposure may also be lower with prolonged infusion of doxorubicin. The effect of drug exposure on distance from the blood supply, shown in the context of a tumour spheroid in figure 10, demonstrates the challenge of delivering drug to the avascular, hypoxic core of a tumour. Such regions are recognized as being chemo-resistant and the model demonstrates the role of PK drug resistance in this context and offers the prospect of modifying the model parameters to improve delivery to such regions and guide the development of novel agents with the necessary characteristics [1,[29][30][31].
From a therapeutic perspective, the most effective mode of administration depends on the threshold of exposure that must be exceeded to kill tumour cells (see figures 7 and 9), as well as the model parameters ( figure 4). As the binding becomes stronger relative to the unbinding, exposure to bound drug increases close to the supply; further from the supply there is, however, an optimal binding ratio, above which the exposure to bound drug starts to decrease (shown in figure 4a and in figure 5). While we have not included an explicit model of cell-kill, this finding implies that drugs which bind avidly may have reduced clinical effectiveness, as they would fail to reach cells distant from the vessel. An example may be the monoclonal antibody trastuzumab, which binds to the cell surface HER2 receptor that is overexpressed in about 15% of breast cancers. Experiments in spheroids suggest that penetration of trastuzumab deeper into the spheroid may indeed be limited by binding to the superficial cell layers [32].
Exposure to bound drug also increases as: (i) the transport rate across the cell membrane increases; (ii) the permeability of the wall of the supplying vessel increases; (iii) the radius of the supplying vessel increases; and (iv) the velocity of the blood flow increases (as illustrated in figure 4). The tumour vasculature is recognized as being disorganized, with vessels that are more permeable, immature and tortuous, with inconsistent diameter and impaired blood flow [2]. Our model suggests that these characteristics of the tumour vasculature may have differing effects on drug exposure. 'Normalizing' the tumour vasculature, as has been shown with the microtubule-targeted cytotoxic eribulin, may enhance the delivery of anticancer drugs [33]. Likewise, some studies demonstrate that bevacizumab increases the number of small vessels and tumour perfusion, allowing a better distribution of paclitaxel when given in combination [34].
When diffusion through the tissue is very rapid, the distribution of drug is uniform (and so is the exposure to bound drug). When diffusion is slow, exposure to bound drug is high close to the supplying vessel but rapidly drops away from the vessel (figure 4c). The clinical significance of this observation is unclear but it is noteworthy that one of the explanations proffered for the very limited efficacy of chemotherapy in pancreatic cancer is the presence of extensive fibrosis, common within such cancers, that is postulated to affect drug delivery [35]. The dependence of uniformity of exposure on transport rate suggests that the raised interstitial pressure observed within tumours will indeed influence drug penetration and distribution in tumours [4]. Similarly, uniform exposure throughout the tissue (and hence good drug penetration) is facilitated by fast transport within the vessel, across the vessel wall and across the cell membrane ( figure 4).
The ultimate aim of predictive modelling of the type presented here is the optimization of drug delivery, but this requires additional knowledge of the dependence of cell survival on the time course of drug concentration. To illustrate the behaviour of our model, we have chosen to indicate the local effectiveness of the treatment by the exposure of the tissue to bound drug. However, the toxicity of doxorubicin is not solely dependent on AUC, particularly at low concentrations, and there is some evidence that peak concentration might be a better indicator of cell-kill ( [22,26] and references therein). It would be straightforward to conduct our study using either AUC or peak values of any of the concentrations in our model, but a better understanding of these dependencies would be required before the model could be deemed to be truly predictive, and therefore effective in aiding treatment optimization. We are also aware that, unlike the peak value, the AUC would continue to increase after the end of the 72 h time window investigated here. Although we believe that we would need a more realistic model of drug clearance for our model to be valid for longer times, we note that running our simulations beyond 72 h did not affect the qualitative behaviour of the results.
The model can be adapted to other cytotoxics or molecularly targeted cancer drugs and extended to include, for example, tumour heterogeneity and different geometries; it can also be applied to other clinical indications. Particularly attractive is the potential to incorporate the model in drug development to select agents most likely to achieve uniform distribution of drug at a relevant concentration throughout the target tissue or tumour. Changing the physico-chemical properties of the drug alters the associated parameters incorporated in the model, such as binding rates, but the values can be identified from drug-specific data. This may allow the identification, from a range of candidate molecules, of the one(s) most likely to address the issue of PK resistance. This could have important implications in terms of the time and expense of preclinical evaluation, as well as addressing the principles of the 3Rs (replacement, reduction and refinement) that have been developed as a framework for humane animal research and subsequently become embedded in guidelines and legislation regulating the use of animals in scientific procedures.

Conclusion
Our model of drug transport in a tumour cord explores the effect of different drug delivery schedules and tumour-binding parameters on tumour exposure to doxorubicin, using parameters determined from experimental data. While acknowledging that the model remains a simplification of a complex system, we believe that such computational models can play an important role in understanding existing and novel cancer therapeutics by generating quantitative and testable insights that inform further experiments in the clinically important area of optimizing drug delivery.
Data accessibility. All data can be accessed at the Dryad Digital Repository http://dx.doi.org/10.5061/dryad.j5s12 [36]. We now know that all of the compartment concentrations tend to zero as t → ∞, so we can consider ∞ 0 δ 1 dM 1 dt + δ 2 dM 2 dt + δ 2 dM 3 dt dt The left-hand side of this equation is equal to zero, because the zero concentrations at t = 0, ∞ give It immediately follows that By a similar argument, Finally, for 1 ≤ i ≤ n, It follows from this that