A multi-scale analysis of drug transport and response for a multi-phase tumour model†

In this article, we consider the spatial homogenisation of a multi-phase model for avascular tumour growth and response to chemotherapeutic treatment. The key contribution of this work is the derivation of a system of homogenised partial differential equations describing macroscopic tumour growth, coupled to transport of drug and nutrient, that explicitly incorporates details of the structure and dynamics of the tumour at the microscale. In order to derive these equations, we employ an asymptotic homogenisation of a microscopic description under the assumption of strong interphase drag, periodic microstructure, and strong separation of scales. The resulting macroscale model comprises a Darcy flow coupled to a system of reaction–advection partial differential equations. The coupled growth, response, and transport dynamics on the tissue scale are investigated via numerical experiments for simple academic test cases of microstructural information and tissue geometry, in which we observe drug- and nutrient-regulated growth and response consistent with the anticipated dynamics of the macroscale system.


Introduction
As a current major cause of death in western countries, cancer represents one of the key challenges for healthcare in the 21st century. For many decades, there has been an extensive effort by experimentalists and clinicians to determine the precise nature of tumour growth and evolution, with the purpose of developing suitable means of treatment. In recent years, there has also been much interest in this field from the mathematical and computational modelling communities (Araujo and McElwain, 2004). Traditionally, this has focussed primarily on the development of mathematical models that describe growth dynamics (and other associated phenomena) so as to enhance qualitative understanding of a given mechanism. More recently, however, a significant effort has been made towards making quantitative predictions for how the disease will progress for an individual patient and, moreover, how successful, or otherwise, a given treatment protocol may be.
Cancer is a disease of mutation (Michor et al., 2004) characterised by a rapid proliferation of damaged cells (Laird, 1964;Folkman and Cotran, 1976). The dynamics and effects of this proliferation are dependent on highly complex physical and biological phenomena across a range of spatial and temporal scales (see, e.g., Alarcón et al. (2005) for a detailed discussion) such as genotypic and phenotypic heterogeneity (Gallaher and Anderson, 2013), the tumour micro-environment (Casciari et al., 1992), and the structure of the surrounding tissue (Rejniak et al., 2013). As a result, among the myriad technical challenges in mathematical and computational modelling of solid tumour growth, the question of how mechanisms occurring on multiple scales may be coupled in a meaningful and efficient manner has garnered much recent interest, see Alarcón et al. (2003Alarcón et al. ( , 2004Alarcón et al. ( , 2005Alarcón et al. ( , 2006, Frieboes et al. (2007), Macklin et al. (2009), Owen et al. (2009Owen et al. ( , 2011, Perfahl et al. (2011), Powathil et al. (2015, , and references therein. In particular, understanding this multi-scale dependence is crucial in predicting the effect of a chemotherapeutic agent on a given tumour. The efficacy of a particular drug is highly dependent on its ability to reach proliferating cells in sufficient quantity to bring about cell damage and a reduction in tumour size (Jain, 1989). This transport is most appropriately modelled on the lengthscale associated with the full extent of the tumour tissue; however, the transport properties of the drug and growth dynamics of the tumour are highly dependent on the evolving microstructural properties of the tumour and its micro-vasculature , Perfahl et al., 2011. In this work, we extend the multi-scale analyses of O' Dea et al. (2015), , and  by employing a multi-phase model of tumour growth, of the type developed in Breward et al. (2002) and , and exploited in the context of vascular tumour growth in Breward et al. (2004) and Hubbard and Byrne (2013), in which we explicitly incorporate microscale dynamics into a system of homogenised partial differential equations (PDEs) for tumour growth, response, and transport.
It is beyond the scope of this article to present a thorough review of the extensive literature associated with the mathematical and computational modelling of solid tumour growth and drug transport. As such, we reference the review papers Araujo and McElwain (2004), Jain (2001), Lowengrub et al. (2010), Preziosi and Tosin (2009b), Roose et al. (2007), and Tracqui (2009), and limit our discussion here, predominantly, to multiphase fluid dynamics models for solid tumour growth of the type considered in this work. The first two-phase models of tumour growth (Breward et al., 2002 decompose the tissue into tumoural and extracellular phases considered as viscous, inertialess fluids. Interphase interactions comprise exchange of mass and momentum via reaction processes regulated by a passive solute and interphase drag, respectively. In subsequent works, additional phases have been added to incorporate tumour vasculature, and enable differentiation between multiple cell species (Breward et al., 2004, Hubbard and. It has been demonstrated in these works that simple multi-phase models are able to reproduce many of the features that characterise solid tumour growth, such as the development of a necrotic core behind a proliferating rim of viable tumour cells. Models of similar structure have also been employed in the study of tumour capsules (Lubkin and Jackson, 2002), the effect of external loading on tumour growth , development of tumour cords and fibrosis Tosin, 2009a, Tosin andPreziosi, 2010), and cell migration (Byrne and Owen, 2004). Multi-phase models have also been applied in the field of tissue engineering, where, typically, viscous fluid flows through a fixed solid phase corresponding to a porous scaffold (see, e.g., Lemon et al. (2006), and Lemon and King (2007)). These models were extended in O' Dea et al. (2008Dea et al. ( , 2010 to study tissue growth in a perfusion bioreactor, in which viscous fluids are driven by externally imposed pressure gradients; and in Ambrosi and Preziosi (2009), where a more general framework of viscoelastic behaviour is introduced to account for cell adhesion effects.
The multi-phase models discussed above are applied at a tissue scale and do not incorporate any microstructural information regarding the tumour or its vasculature (though via volume averaging one can show that multi-phase models of the type described above are physically meaningful, see Drew (1971Drew ( , 1983). There is a growing literature associated with hybrid models of tumour growth, in which discrete populations of cells are simulated with continuous fields of nutrient, and possibly drug and angiogenic growth factor, see Alarcón et al. (2003Alarcón et al. ( , 2004Alarcón et al. ( , 2005Alarcón et al. ( , 2006, Chaplain et al. (2006), Frieboes et al. (2007), Macklin et al. (2009), Owen et al. (2009), Perfahl et al. (2011), and Powathil et al. (2015. The hybrid nature of these models allows the incorporation of effects that occur on multiple spatial and temporal scales. In the wider field of transport in biological tissues, however, there is much interest in the application of homogenisation techniques as a means to incorporate effects over multiple scales. These techniques allow the derivation of effective ordinary or PDEs at the tissue scale, that directly incorporate explicit information regarding the microscale in a mathematically precise manner. Examples employing these methods across a range of biological applications include Band and King (2012), Fozard et al. (2010), O'Dea and King (2011, Ptashnyk and Chavarría-Krauser (2010), Roose (2010), andTurner et al. (2004). However, of particular interest to the current work is , in which flow and transport equations in a solid tumour and its vasculature are homogenised to obtain a macroscale description of drug transport; and the more recent works of O'Dea et al. (2015) and Penta et al. (2014) which consider the homogenisation of models for growing tissues.
The analysis in the current work extends the classical homogenisation of flow and transport phenomena along similar lines to the analyses in O' Dea et al. (2015) and Penta et al. (2014), wherein the growth, flow, and transport dynamics of porous tissues are investigated. Here, however, we additionally consider the dynamics of a multi-phase mixture and, as a means of incorporating interstitial growth, we further consider the dynamics of a free interface on which phase transition may occur. This article is organised as follows. In Section 2, we introduce a generic, multi-phase PDE model governing drugand nutrient-regulated tumour growth at the microscale, together with transport of the aforementioned passive solutes. In Section 3, we present the derivation of a corresponding macroscale continuum model via a multi-scale homogenisation method. In Section 4, we present an illustrative numerical experiment of the microscale system, for an academic example of the geometry. In Section 5, we describe a specific model of tumour growth and response, and present a series of numerical experiments for representative model data on academic examples of geometry, parameterised using the data obtained in the microscale numerical experiment in Section 4. Finally, in Section 6, we make some concluding remarks and highlight ongoing and future work.

Model description
In this section, we present an idealised model of microscale tumour growth; for the sake of generality and clarity of presentation, the model given in this section, as well as the analysis presented in Section 3 and numerical experiments on the microscale presented in Section 4, is intentionally generic. Detailed biological motivation is therefore kept to a minimum, except where necessary for rationalising specific model development choices. Given this generality, our analysis may be applicable to fields of study other than solid tumour growth (such as phase transition in geophysical flows, see, e.g., the introduction to Morland and Sellers (2001) and references cited therein). Our biological motivation for the underlying microscale structures that we consider in Section 4, arises from the use of similar structures elsewhere in the chemotherapeutic drug transport literature, e.g., Rejniak et al. (2013) which considers drug transport through ovarian tumour tissue. For the numerical experiments in Section 5, we return to a biological setting by considering a specific model for avascular tumour growth, discussed therein.
We consider a region of tissue containing a tumour as an idealised porous medium in R d , d = 2, 3. The medium is represented as a highly connected material with a spatially periodic microstructure consisting of a multi-component mixture, whose dynamics is regulated by the concentrations of passive solutes (representing a chemotherapeutic drug and nutrient, respectively), and saturated with a viscous Newtonian fluid which we refer to as interstitial fluid. We assume that the material may be characterised by two distinct lengthscales: , corresponding to the scale of the periodicity of the microstructure, and L corresponding to the scale of the tissue; we refer to these as the microscale and macroscale for the remainder of this article. We denote by Ω the microscale domain corresponding to periodic microstructure, and its boundary by ∂Ω.
We further partition Ω into two subdomains Ω 1 and Ω 2 corresponding to the multicomponent mixture (henceforth referred to as the mixture for brevity) and interstitial fluid, respectively, with the interface between Ω 1 and Ω 2 denoted by Γ ; the unit inward normal to Ω 1 on Γ is denoted by n and the unit tangent(s) to Γ are denoted by τ . Growth is represented by movement of the interface Γ .
A schematic diagram of the porous material described here is shown in Figure 1. We denote by Ω L the macroscale domain associated with the tissue. Finally, we assume that there is a strong separation of the two lengthscales, i.e., we introduce a small dimensionless parameter 0 < 1 defined by = L . (2.1) We remark that macroscale variation within Ω is permitted, though we do not pursue this in the current work.
There is an extensive literature related to the classical homogenisation of flow and transport in porous media in both physical and biological applications, see, e.g., Keller (1980), Mei and Auriault (1991), Rubinstein (1987), and Tartar (1980). The analysis we present here represents an extension of the traditional homogenisation of flow and transport phenomena and the recent attempts to apply these ideas to growing material in O' Dea et al. (2015) and Penta et al. (2014), by explicitly incorporating the dynamics of the multi-phase mixture in Ω 1 , the free interface Γ , and allowing phase transition on the singular surface Γ as a means of incorporating interstitial growth.

Governing equations on the microscale
In this section, we define the systems of equations in Ω governing the motion of the mixture, interstitial fluid, and their interface, as well as the passive transport of drug and nutrient.

Equations governing flow and transport in Ω 1
The structure of the model for the mixture presented in this section follows closely the structure of the multi-phase model of tumour growth employed in Hubbard and Byrne (2013). We consider Ω 1 as a mixture of N θ interacting fluid phases. In our biological context, each phase represents a biological material, e.g., cell species, water, or extra cellular matrix. We denote the volume fractions, velocities, pressures, densities, kinematic viscosities, and bulk viscosities of the different phases by θ i , u 1,i , p 1,i , ρ 1,i , μ 1,i , and λ 1,i , respectively, for 1 6 i 6 N θ , where we consider θ i , u 1,i , and p 1,i as spatially and temporally varying dependent variables, and ρ 1,i , μ 1,i , and λ 1,i as constant model parameters. Conservation of mass for each of the cell species yields the following PDE for the volume fraction θ i : where S i represents the net mass transfer from phase j to phase i, defined as and s ij denotes the volume flux from phase j into phase i. We specify that mass is conserved across all phases, thus yielding the condition (2.4) and that there are no voids in Ω 1 , i.e., In the absence of inertia, conservation of momentum may be written as where σ i,1 and F i denote the stress and momentum sources, respectively, for each phase. We remark that the assumption that we may neglect inertial terms is not required a priori for the proceeding analysis. Under the dimensionless scalings set out in Section 2.2, we see that inertial terms are in fact relegated to O 2 , see O'Dea et al. (2015). However, we present the derivation in this manner for the sake of concision, and for consistency with similar multi-phase formulations, Hubbard and Byrne (2013) for instance. The equation for the composite velocity for the mixture is obtained by summing (2.2) over all N θ phases and applying (2.5) to yield We constitutively specify the stress tensors σ 1,i by σ 1,i = −p 1,i I + μ 1,i D(u 1,i ) + λ 1,i ∇ · u 1,i I , for 1 6 i 6 N θ , (2.8) where D denotes the rate of strain tensor We further specify the momentum sources constitutively as where d ij denotes the drag coefficient associated with the relative motion of phases i and j. To simplify our microscale formulation, we now consider the limit in which the interphase viscous drag is large (d ij 1). This assumption is motivated by works such as O' Dea et al. (2008); further, we refer to Franks and King (2003), for example, for historic use of such a model in the same broad context as this work. Under this assumption, each phase has the same velocity u 1 , pressure p 1 , and stress σ 1 . As such, equations (2.2), (2.6), (2.7), and (2.8) simplify to where we have implicitly assumed that each phase has kinematic viscosity μ 1 and bulk viscosity λ 1 . We remark that this adoption of the large drag limit enables the multi-scale analysis that follows, through its simplification of the momentum and mass conservation equations (2.12) and (2.13). A corresponding analysis incorporating distinct phase velocities forms important future work. In addition to considering the motion of the mixture in Ω 1 , we consider the transport of two passive solutes, whose concentrations are denoted c and n. In our biological context, we consider c to correspond to a drug, e.g., doxorubicin or cisplatin, and n to correspond to a nutrient, e.g., glucose or oxygen. The equations governing the evolution of the concentration of these transported species in Ω 1 are given by and ∂n ∂t + ∇ · (u 1 n) = ∇ · D 1,n ∇n − γ 1,n n, (2.16) where D 1,c denotes diffusivity of c through the mixture, γ 1,c denotes the consumption rate of c, and ξ c denotes its decay rate. Analogous notation is employed for parameters related to the transport of n, but are denoted with an appropriate subscript. In general, the diffusivity D, and rates γ and ξ, may have functional dependence on θ i , c, and n; in the following, however, we assume for simplicity that the diffusivities and decay rates are constant model parameters, but the consumption rate may be dependent on θ i , c, and n.

Equations governing flow and transport in Ω 2
We model the interstitial fluid in Ω 2 as a viscous Newtonian fluid and denote its velocity, pressure, density, and kinematic viscosity by u 2 , p 2 , ρ 2 , and μ 2 , respectively, where we consider u 2 and p 2 as spatially and temporally varying dependent variables, and ρ 2 and μ 2 as constant model parameters. As such, we assume the motion of the interstitial fluid is governed by the incompressible Navier-Stokes equations, given by Here, for reasons of convention, we retain inertial terms in (2.6). This serves to highlight their subsequent relegation to O 2 under the dimensionless scalings given in Section 2.2 (see also Section 2.1.1). Once more, we consider the transport of passive solutes whose concentrations are denoted c and n. The equations governing the evolution of the concentration of these transported species are given by where D 2,c and D 2,n denote the diffusivities through the interstitial fluid of the solutes.

Interface conditions on Γ
In order to fully couple the systems of equations given in Sections 2.1.1 and 2.1.2, we must specify appropriate conditions for stress and velocity of the fluids, mass balance and concentration of the passive solutes, and the motion of the free boundary. To this end, we first define the ± sides of the free interface Γ . Recalling the definition of the unit normal n oriented into Ω 1 ; then for a scalar quantity a defined in both Ω 1 and Ω 2 , we define the quantities a ± by If we let u Γ denote the velocity of the free interface Γ , then the mass balance conditions for c and n across Γ are given by where a ∈ {c, n}. Furthermore, we impose continuity of concentration of the solutes across the interface, i.e., for a ∈ {c, n}. Letting σ 2 denote the stress tensor for the interstitial fluid in Ω 2 , i.e., then, following Morland and Sellers (2001), we impose the following stress condition for a multi-phase mixture with phase transition on the singular surface Γ , (2.24) We note that, consistent with our momentum descriptions (2.6) and (2.17), inertial terms are omitted from the left-hand side of (2.24) but retained on the right-hand side; however, as already noted, these are in any case relegated to O 2 under the scalings given subsequently. Additionally, we impose continuity of tangential velocity across Γ , i.e., for all tangents τ of Γ . The mass balance condition for the fluid moving across the interface Γ is given by where G is the phase transition rate on Γ , and ρ 1 is the volume fraction averaged density in Ω 1 , defined by Finally, assuming that the interface Γ is described by the zero level set of some function F(x, t), its kinematics are then described by the condition For the sake of generality, we do not provide specific constitutive assumptions on G until we specify the tumour growth model employed in the numerical experiments of Section 5. We remark, however, that constitutively specifying G in the interface condition (2.26) significantly simplifies the description of the interfacial motion. This form arises naturally from a boundary layer analysis of the more general multi-phase description, in the large-drag/fast-consumption limit (see King and Franks (2007)), the details of which are beyond the scope of the current work.

Boundary conditions on ∂Ω
In order to complete the system on the microscale, it remains only to specify boundary conditions on ∂Ω. As we have made the assumption of periodic microstructure and strong separation of the micro-and macroscale, we naturally impose periodic boundary conditions for all components of the fluid stresses and velocities.

Non-dimensionalisation
We non-dimensionalise the independent variables x and t with respect to the characteristic lengthscale of the periodic microstructure and advection across the tissue, respectively, i.e., where U denotes the characteristic velocity scale and * denotes the dimensionless variable. This choice of characteristic time is necessary for the proceeding multi-scale analysis, as it ensures that inertial terms are relegated to O 2 , as discussed in Sections 2.1.1 and 2.1.2, thus allowing the Darcy flow assumption of Section 3.1. Importantly, this choice of characteristic timescale ensures that the leading order equations for the flow in Ω 1 and Ω 2 are steady and the free interface Γ is quasi-stationary, thus relegating the growth dynamics to higher order in and simplifying the analysis. We remark that the analysis of O(1) growth dynamics or O( ) differences in material densities (through choice of a timescale t = t * /S, say) significantly complicates the resulting system of equations. In these cases, the spatial homogenisation techniques considered herein will likely need supplementing with a multiple timescale analysis. Such considerations are beyond the scope of this current work.
We non-dimensionalise the dependent variables in Ω 1 and Ω 2 with We non-dimensionalise the phase source terms in Ω 1 with and we further define the Reynolds number for macroscale flow of interstitial fluid by The dependent variables defined in all of Ω are non-dimensionalised with c = Cc * and n = Nn * , where C and N are characteristic concentrations of the passive solutes, respectively. Finally, we define the following dimensionless parameters For the remainder of this article, we work exclusively with the dimensionless variables and parameters 1 and forgo the notation * for clarity of presentation. We proceed by setting out the full system of non-dimensional PDEs.
2.2.1 Flow and transport in Ω 1 We first consider the dimensionless system of equations that describes the mixture in Ω 1 subject to the assumption of strong interphase drag. This is given by The dimensionless form of the transport equations in Ω 1 are given by (2.36)

Flow and transport in Ω 2
The dimensionless equations that describe the motion of the interstitial fluid are given by Note also that the dimensionless form of the constitutive assumption on the fluid stress in Ω 2 becomes The dimensionless form of the equations governing transport of the passive solutes in Ω 2 are given by and Finally, the non-dimensional version of the kinematic condition (2.28) is given by (2.44)

Multi-scale analysis
We now analyse the dimensionless system of equations given in Sections 2.2.1-2.2.3 with the aim of deriving a description of the behaviour at the tissue scale which explicitly incorporates the structure and dynamics of the microscale in a manner analogous to that presented in Burridge and Keller (1981), O'Dea et al. (2015), Penta et al. (2014), and . We define X to be the dimensionless macroscale spatial variable. By equation (2.1), this is related to the microscale variable x via X = x. Under the assumption of strong separation of scales, we may expand all dependent variables Ψ in multiple-scales form via an expansion of the form Under this coordinate transformation ∇ becomes where ∇ x and ∇ X represent differentiation with respect to the microscale and macroscale spatial variables, respectively. We proceed by substituting multiple-scale expansions of the form (3.1) for each of the dependent variables into the non-dimensional system of equations set out in Sections 2.2.1-2.2.3 to obtain a system of PDEs at increasing orders of . We then analyse the systems obtained at each order in , with the aim of obtaining a description of the growth and transport dynamics at the macroscale.

Microscale behaviour
In this section, we derive systems of equations for flow and transport in Ω to O(1) and O( ), and introduce a Darcy flow assumption on the leading-order velocities and first-order pressures.

Equations at O(1)
At leading order in , the equations for the flow in Ω 1 are given by From equations (3.4) and (3.5), we see that p (0) 1 and θ (0) i are locally constant in Ω 1 , i.e., they exhibit no dependence on the microscale variable x and may be written as and where we use the notation(·) to denote microscale independence. If we define the linear transport operator for i ∈ {1, 2}, a ∈ {c, n}, and k ∈ N 0 , the transport equations in Ω 1 reduce to The equations for the flow in Ω 2 reduce to ∇ x p 2 = 0, (3.10) From equation (3.10), we see that p (0) 2 is locally constant in Ω 2 , i.e., p (0) 2 =p 2 (X , t).
Using the definition (3.7), the transport equations in Ω 2 are given by L 0 2,n n (0) = 0. (3.13) The stress condition on Γ is given by and if we substitute the constitutive assumptions for stresses given by equations (2.39) and (2.34) into (3.14), we see thatp and, therefore, the leading-order pressurep(X , t) is constant across all of Ω. Continuity of tangential velocity is given by The kinematic condition reduces to which, given that F defines a proper boundary Γ via its zero level set, implies that u (0) Γ · n = 0, i.e., at leading order in the free interface is stationary. This observation further implies that the conservation of fluid mass across Γ given by (2.43) becomes i.e., at leading order there is no mass flux across Γ . The interface conditions for the transportable species are given by Remark 3.1 We will subsequently show in Section 3.2 that despite phase transition on Γ and mass source terms in Ω 1 not occurring until O( ) under the scaling given in Section 2.2, their effects are evident in the leading-order macroscale equations for both flow and transport.

Equations at O( )
At first order in , the equations for the flow and transport in Ω 1 are given by (3.24) in which we have made appropriate substitutions employing equation (2.31). The equations for the flow and transport in Ω 2 are given by (3.28) The stress condition on Γ reduces to If we substitute the constitutive assumption for σ 2 given by equation (2.39) and σ 1 given by equation (2.34) into equation (3.29), we obtain the condition where D x denotes the rate of strain tensor for differentiation with respect to the microscale variable x. Finally, the interface conditions for the fluid are given by and the interface conditions for the passive solutes are given by for a ∈ {c, n}.
3.1.3 Ansatz for u (0) and p (1) The system of equations given in Section 3.1.1 is not sufficient to specify the flow in either Ω 1 or Ω 2 . As such, we proceed by assuming suitable forms for u (0) 1 , u (0) 2 , p 1 , and p 2 and substituting them into the appropriate equations at O(1) and O( ) to obtain a pair of coupled, steady Stokes problems that may be solved in order to obtain a description of the flow in Ω, see Davit et al. (2013), O'Dea et al. (2015), Rubinstein and Torquato (1989), and .
In light of the linearity of (3.3), (3.11), (3.21), and (3.25), an appropriate form for the leading-order microscale flow may be obtained by making the following assumptions on u (0) i and p (1) t) is the mean value of the first-order pressure in Ω i . The permeability tensor K i and vector a i exhibit both microscale and macroscale dependence, but are independent of time to leading order. This ansatz lets us quantify the microscale dependence on u (0) i and p (1) i via a steady problem on the periodic microstructure for K i and a i . Substituting equation (3.36) into either equations (3.3) and (3.21) or equations (3.11) and (3.25), for flow in Ω 1 or Ω 2 , respectively, and exploiting linearity, we obtain the following steady Stokes problem for K i and a for i = 1, 2. In order to obtain suitable conditions to couple the Stokes problems on Γ , we substitute equation ( Finally, we substitute into the O(1) condition on the tangential velocities given by equation (3.15), yielding Remark 3.2 We note that had we not specified stress and velocity conditions on Γ in our underlying microscale problem, (3.38)-(3.40) would be the natural choice of interface conditions for a pair of coupled tensor Stokes problems with no exchange of fluid mass between Ω 1 and Ω 2 .
The system of PDEs defined by (3.37)-(3.40) together with periodic boundary conditions on ∂Ω is not sufficient to uniquely specify the microscale variables a 1 and a 2 . As such, we specify the additional constraint, Ω1 (a 1 ) i dx + Ω2 (a 2 ) i dx = 0, for 1 6 i 6 d, (3.41) in order to guarantee the uniqueness of a 1 and a 2 . We emphasise that the microscale flow problem described here is steady because, under the choice of dimensionless scaling set out in Section 2.2, the dynamics of the microscale system are relegated to O( ).

Macroscale behaviour
In order to obtain the macroscale flow and transport equations, we take spatial averages of the leading-order microscale flow and transport equations. To this end, we define the following integral averages 42) and the porosity of the material, φ, by (3.43)

Macroscale flow
The leading-order macroscale flow description is obtained by averaging the microscale equations for the velocity in equation (3.36) to obtain Averaging equations (3.20) and (3.26) over Ω 1 and Ω 2 , respectively, and applying the divergence theorem yields If we define the combined (average) velocity in Ω at leading order by (3.49)

Macroscale transport of passive solutes
The macroscale equations governing the transport of drug and nutrient are obtained by averaging the Fredholm alternative solvability conditions for the O( ) transport equations (see O'Dea et al. (2015), , and ) to obtain the following: where R c and R n denote the terms (3.53) In this section, we have derived a macroscale model for flow and transport across Ω L that incorporates information regarding the structure of the microscale. The microscale model consists of a pair of coupled steady tensor Stokes problems for permeability tensors on the subdomains Ω 1 and Ω 2 . These permeability tensors are then employed to parameterise a macroscale model comprising a Darcy system for flow and advectionreaction PDEs for the mixture component volume fractions, and concentrations of the passive solutes. Through the choice of dimensionless scaling, inertial terms and dynamics on the microscale are relegated to higher order, thus yielding quasi-stationary linear flow problems to leading order. Phase transition on the singular surface Γ and the mass sources for the mixture components nevertheless induce a flow at the macroscale to leading order and, moreover, affect the macroscale evolution of the leading-order mixture volume fractions and the concentrations of the passive solutes. In biological applications, we observe that density differences between phases may be small and are typically neglected in the tissue growth literature. However, in the proceeding sections, we retain non-infinitesimal density differences between phases for the sake of generality (that may be observed in physically motivated, as opposed biologically motivated, models), so as to highlight the full resultant dynamics of the proposed macroscale model.

Numerical experiments on the microscale
In this section, we present a numerical experiment demonstrating the numerical approximation of the solution to the microscale system given in Section 3.1.3 for a representative, academic example of the geometry. This example is chosen to highlight suitable means of computing microscale information that will then be employed to parameterise representative numerical experiments of the macroscale system in Section 5, and provide details of appropriate numerical methods for discretising the microscale system. It is, however, beyond the scope of this article to present full details of the discretisation of the system of PDEs that define the microscale model and, as such, we refer to our companion article Collis et al. (2016) which fully describes the discretisations employed here. All numerical experiments were carried out with the AptoFEM finite element package (Antonietti et al., 2015).
In this current work, we consider only academic examples of microscale geometry, as the main contribution of this article is the multiple scales analysis of the generic multiphase model presented in Section 3, and its subsequent application to a model of tumour growth and response in Section 5. Moreover, it is not appropriate to consider a proper treatment of the techniques associated with image recognition and segmentation that are necessary to reliably incorporate physiological microscale data in this current work. For relevant examples of microscale geometry from the field of tissue engineering, we refer the reader to Morris et al. (2014) and Visser et al. (2015).
For computational simplicity, we would, ideally, consider only two-dimensional calculations. However, in order to obtain meaningful average permeabilities for both subdomains, we are required to consider a three-dimensional problem. As such, we consider a threedimensional microscale problem, but restrict the resulting data employed to parameterise the macroscale problem to two spatial dimensions, in order to consider a two-dimensional macroscale problem in Section 5.
The system of tensor Stokes problems is solved by decomposing the tensor problem into three standard Stokes problems for the components K ij and a i for 1 6 i, j 6 d. We then approximate the solution to these problems using a discontinuous Galerkin (dG) finite element (FE) method following Toselli (2002), with an additional penalisation term to impose the no penetration condition (3.39). The domain Ω is discretised into a uniformhexahedral mesh containing 30 × 30 × 30 elements, and the dG FE approximation was computed employing polynomials of degree (2, 2, 2, 1) for the components of K and a, respectively. Figures 3 and 4 show representative numerical approximations for K 33 on the geometry shown in Figure 2. As K is a permeability and appears only in our macroscale model as an averaged quantity, it is difficult to interpret Figures 3 and 4 in a physically meaningful way with respect to both the micro-and macroscale; we therefore include these figures here only to provide the reader with an indication of the variation in K 33 across the microscale domain, and refer to our companion article Collis et al. (2016) for more detail. Here, it suffices to note that due to the symmetry of Ω, Ω 1 , and Ω 2 , the computation of K 13 , K 23 , K 33 , and a 3 is sufficient to fully specify K and a throughout Ω, from which we may calculate the spatially averaged permeabilities as We note that the off-diagonal entries of K 1 1 and K 2 2 are, effectively, zero: under the choice of geometry shown in Figure 2, we would expect K 1 1 and K 2 2 to be isotropic, to within the discretisation error in the dG FE method, and the results (4.2) and (4.3) are consistent with this.

Numerical experiments on the macroscale
In this section, we depart from the rather generic multi-phase setting of Sections 2-4 by utilising a specific model of tumour growth and response that is consistent with the analysis presented in the preceding sections. We subsequently demonstrate the solution of the macroscale system defined in Sections 3.2.1-3.2.3 for this model, employing representative parameter values, boundary and initial data, and geometry as well as microscale data obtained from the numerical experiments detailed in Section 4.

Model of tumour growth and response
We now proceed by specifying the model of tumour growth employed in the numerical experiments. For simplicity, we closely follow that presented in Hubbard and Byrne (2013), with additional terms corresponding to cell response to chemotherapy and removal of a vascular phase that was of relevance to that study. As such, we consider three phases in the mixture, namely normal cells, tumour cells, and extra-cellular material (ECM), denoted by i = 1, 2, 3, respectively. We specify the sources, S i , in the microscale model to be of the form S 1 = ρ 1,3 k m,1 θ 1 θ 3 n n P + n cell birth − ρ 1,1 k d,1 θ 1 n 1 + n n 2 + n nutrient-regulated cell death − ρ 1,1 k c,1 θ 1 c c P + c drug-regulated cell death , (5.1) S 2 = ρ 1,3 k m,2 θ 2 θ 3 n n P + n − ρ 1,2 k d,2 θ 2 n 1 + n n 2 + n − ρ 1,2 k c,2 θ 2 c c P + c , (5.2) S 3 = ρ 1,1 k d,1 θ 1 + ρ 1,2 k d,2 θ 2 n 1 + n n 2 + n + ρ 1,1 k c,1 θ 1 + ρ 1,2 k c,2 θ 2 c c P + c − ρ 1,3 k m,1 θ 1 + k m,2 θ 2 θ 3 n n P + n . (5.3) For the sake of brevity, we refer the reader to Hubbard and Byrne (2013) for a full description of the terms associated with cell birth and nutrient-regulated cell death, and associated parameters. For the tumour response to chemotherapy, we assume the death rate is a monotonically increasing and saturating function of c increasing from zero to a maximal rate k c,i as c → ∞, for i = 1, 2. Furthermore, we specify that the consumption rate of drug, γ 1,c , and nutrient, γ 1,n , are given by γ 1,c =γ 1,c (θ 1 + θ 2 ), (5.4) γ 1,n = k γ,1 θ 1 + k γ,2 θ 2 baseline consumption + k m γ,1 θ 1 + k m γ,2 θ 2 θ 3 n n P + n cell birth consumption . (5.5) Once more, we refer to Hubbard and Byrne (2013) for a complete description of the terms and associated model parameters. In the following numerical experiments, we consider two sub-cases: first (in Section 5.2.1), we specify the constitutive assumption on the phase transition on the singular surface, Γ , as Under this assumption, macroscale flow is driven entirely though differences in density between each phase in the mixture. Then (in Section 5.2.2), we consider an alternative constitutive assumption on G, that will allow us to recover some of the qualitative features observed in multi-phase models of tumour growth elsewhere in the literature. Once more we remark that density differences between phases in the mixture are retained here for generality and in order to emphasise the full emergent dynamics of the proposed macroscopic model. We further recall the strong drag assumption between the phases given in Section 2.1.1. Finally, we specify appropriate boundary and initial conditions for the macroscale problem on Ω L . We assume that the computational domain is surrounded by a region of healthy tissue in which drug and nutrient are uniformly distributed. We further assume that the interstitial pressure outside of the computational domain is known. Under these assumptions, appropriate boundary conditions are given bȳ where Λ -(b) denotes the inflow boundary with respect to a vector field b, Θ denotes the volume fraction of healthy cells in the tissue surrounding the tumour, p ∂ denotes the known pressure in the host tissue, and c ∂ (t) is a time-varying drug concentration corresponding to multiple rounds of treatment. Values of Θ, p ∂ , and c ∂ are specified in Section 5.2.

Numerical experiments
We approximate the solution of the flow problem (3.48) at each point in time using the lowest order, mixed Raviart-Thomas (RT) dG FE method for (u,p) (see Raviart and Thomas (1977) and Brezzi and Fortin (1991)). The solutions to the transport equations (3.49), (3.50), (3.51) are approximated on a shape regular triangulation T h using a dG FE method as described in Houston and Süli (2001), with polynomial degree 0. The timestepping is performed using an explicit Euler method, linearised about the solution at the previous time point. The initial conditions for the volume fractions for the components of the mixture are given bȳ for X ∈ (X 1 , X 2 ) : X 2 1 + X 2 2 6 a 2 , and θ 1 (X , 0) = Θ, where the off-diagonal components in equations (4.2) and (4.3) are sufficiently small that they may be treated as zero.

No interfacial phase transition
In this section, we consider the case of no interfacial phase transition, i.e., we specify that We define the macroscale geometry Ω L = (0, 1) 2 . Ω L is then discretised into a shape-regular conforming triangulation, T h , shown in Figure 5. Furthermore, we specify the boundary condition for the pressure as  Table 1. In Figures 6(a) and (b), we observe an initial growth of the tumour, resulting in a reduction in the volume fractions of healthy tissue and ECM. In Figures 6(c)-(e), we then observe a reduction in the tumour as a result of drug-induced cell death. The presence of drug also results in the death of healthy tissue, though at a lower rate than that of tumour. This cell death results in an increased volume fraction of ECM. We further remark that the net flow from left to right is induced by equation (5.20). In Figures 7(a), (c), (e), and (g) we observe the transport of drug across Ω L . Comparing these to Figure 6, we see that the regions of high drug concentration correspond to regions with decreasing volume fractions of tumour. In Figures 7(b), (d), (f), and (h), we observe the replenishment of the nutrient across the domain as a result of the induced flow. Two features we highlight here are the increased consumption of nutrient by the tumour compared to the healthy tissue as a result of greater rate of mitosis, as evident from Figures 7(b) and (d), and the regions of high nutrient concentration apparent in 7(f) and (h) as a result of reduced consumption resulting from the lower volume fraction of tumour and local compression in the underlying velocity field caused by drug-induced cell death. Figure 8 shows the tumour mass and the boundary data for the drug concentration over the full time of the simulation. In this figure, we clearly see the initial increase and subsequent decrease in tumour mass resulting from the introduction of the drug. We remark that there is a clear lag between the peak in the boundary data and the change in the growth dynamics of the tumour corresponding to the time taken for the drug to be transported from the boundary to the region of high tumour volume fraction, cf. Figures  6 and 7.
We note in passing that, under the assumption of strong interphase drag and no interfacial phase transition, a velocity field consistent with a growing tumour (in the   context of this model) and that transports drug and nutrient from the boundary across the domain, is not readily available. Indeed, we have been unable to obtain such a flow in our numerical experiments (though a more complete parameter survey may reveal a choice for which this is realised); for this reason, we have imposed a flow across the domain via suitable boundary conditions on the pressure. This is due to the advective

Interfacial phase transition
We now consider a model incorporating interfacial phase transition. As such, we are required to specify constitutively the mass transfer function G. To this end, we suppose that there is some equilibrium volume fraction of ECM, denoted θ E , such that there is phase transition across Γ in order to maintain this equilibrium. Under this assumption, a natural choice for the form of G is given by which we adopt for this numerical experiment. As noted above, for G = 0 the requirement to impose a flow on the macroscale system is removed and we therefore specify the boundary condition for the pressure as p ∂ = 0 (5.22) for illustrative purposes (and since p ∂ = 0, we are less restricted in the choice of physically meaningful domains). The boundary data for the drug is still given by equation (5.16). We further consider an alternative macroscale geometry to that employed in Section 5.2.1, namely Ω L = X : X 2 1 + X 2 2 6 0.5 2 , and discretise this into a shape-regular conforming triangulation, T h , shown in Figure 9. Note that there is a discrepancy between the triangulation and the true domain, but we assume that this is small and has no qualitative effect on the results. Figures 10 and 11 show the volume fraction of each of the components of the mixture and the concentration of the passive solutes, respectively, at a range of times for the parameter values given in Table 2. The behaviour exhibited in these figures is qualitatively similar to that presented in Figures 6 and 9; however, we re-emphasise that we are now able to observe transport of drug and nutrient into the domain without inducing a flow via the pressure boundary conditions, for reasonable initial configurations of the mixture.
In Figures 10(a)-(c) (central panels), we observe an initial growth of the tumour; then in Figures 10(c)-(e), we start to observe the effect of the drug on the tumour phase. We remark that, however, in Figures 10(b) and (c), we additionally observe the effect of drug on the healthy tissue.
In Figure 11, we see qualitatively similar behaviour to that described in Section 5.2.1 (see Figure 7). However, we remark that the assumption on G enhances the compressive nature of the induced flow and hence we see higher and lower concentrations of drug and nutrient compared with those observed in Section 5.2.1.

Conclusions
In this article, we have performed a spatial homogenisation of a multi-phase model for avascular tumour growth and response to chemotherapy, based on the model employed in Hubbard and Byrne (2013). The multiple-scales technique that we exploit has been widely employed in homogenisation of flow and transport in porous media in Keller (1980), Mei and Auriault (1991), Rubinstein (1987), and Tartar (1980), and more recently to simple models of growing tissue in O' Dea et al. (2015) and Penta et al. (2014). Here, we seek to extend these ideas to a more complex description of the underlying tissue dynamics, incorporating multiple phases as well as microstructural features. We obtain a tissue-scale description of tumour growth and response, and transport of drug and nutrient comprising a system of advection-reaction PDEs that are partially parameterised by the solution to a tensor Stokes problem on a microscopic periodic unit. As a consequence of the timescale on which we perform our analysis, growth and transport processes are relegated to O( ), and the underlying tissue microstructure is therefore quasi-steady. Nevertheless, they directly influence the macroscale flow, transport, and tumour evolution dynamics.
We demonstrate the dynamics of the macroscale system via numerical examples employing academic examples of micro-and macroscale geometry. We refer the interested reader to our companion article Collis et al. (2016) for full details of the numerical approach.
There are several natural extensions to the work considered in this article. The foremost being that we have employed academic examples of microstructural geometry to illustrate the model behaviour; in order for this model to be of direct relevance to experimentalists or clinicians, it is crucial that realistic geometries are employed for both the micro-and macroscale models. With regards to extensions of the analysis presented here, current ongoing work surrounds the removal of the strong drag assumption (in order to more comprehensively consider interstitial growth) and the supplementation of this current work with a multiple timescales analysis in order to consider O(1) microscale growth in a computationally tractable macroscale formulation. The former may illuminate appropriate forms for G to yield greater realism in the observed macroscale model dynamics, within our simplified formulation.