Multi-scale modelling of nanoparticle delivery and heat transport in vascularised tumours

We focus on modelling of cancer hyperthermia driven by the application of the magnetic field to iron oxide nanoparticles. We assume that the particles are interacting with the tumour environment by extravasating from the vessels into the interstitial space. We start from Darcy’s and Stokes’ problems in the interstitial and fluid vessels compartments. Advection–diffusion of nanoparticles takes place in both compartments (as well as uptake in the tumour interstitium), and a heat source proportional to the concentration of nanoparticles drives heat diffusion and convection in the system. The system under consideration is intrinsically multi-scale. The distance between adjacent vessels (the micro-scale ) is much smaller than the average tumour size (the macro-scale ). We then apply the asymptotic homogenisation technique to retain the influence of the micro-structure on the tissue scale distribution of heat and particles. We derive a new system of homogenised partial differential equations (PDEs) describing blood transport, delivery of nanoparticles and heat transport. The new model comprises a double Darcy’s law, coupled with two double advection–diffusion–reaction systems of PDEs describing fluid, particles and heat transport and mass, drug and heat exchange. The role of the micro-structure is encoded in the coefficients of the model, which are to be computed solving appropriate periodic problems. We show that the heat distribution is impaired by increasing vessels’ tortuosity and that regularization of the micro-vessels can produce a significant increase (1–2 degrees) in the maximum temperature. We quantify the impact of modifying the properties of the magnetic field depending on the vessels’ tortuosity.


Introduction
Tumours are cancerous tissues consisting of different cells that interact with each other by heterotypic interaction i.e. involving both cancerous and non-cancerous cells and constituents. They are characterized by different hallmarks that enable them to develop and growth in the body. The cancer cells preserve the chronic proliferative signals and avoid the suppressor genes, which regulate cell duplication. While normal cells replication is limited, cancer cells can grow in an uncontrolled way thus producing the tumour mass. In addition, tumours are capable of obtaining nutrients and oxygen from the angiogenic 334 T. AL SARIRI AND R. PENTA by the Stokes' problem, while we assume that interstitial transport is described by Darcy's law. The governing equations describing drug transport are of advection-diffusion type in the vessels and of advection-diffusion-reaction type in the tumour interstitium. Diffusive and advective heat transport in both compartments is likewise formally represented by a double advection-diffusion-reaction model. The vessel's wall is modeled as a porous semi-permeable membrane, so as to allow the interplay of fluid, mass, and heat between compartments.
In this work we address the sharp length scale that exists between the typical intercapillary distance and the average tumour size by means of the asymptotic homogenisation technique, as summarized e.g. in the works by Bakhvalov & Panasenko (1989); Cioranescu & Donato (1999); Davit et al. (2013); Taffetani et al. (2014); Penta & Gerisch (2017). This strategy provides the upscaled parameters at the macro-scale on the basis of the micro-scale geometry and structural properties. In particular, the resulting macro-scale governing equations encode the crucial role of the micro-scale structure in the coefficient of the model. This approach is motivated by the fact that real-world problems typically involve interactions between a variety of physical three-dimensional systems characterized by different properties. In general, it is practically impossible to resolve all the micro-scale details in three dimensions. Furthermore, experiment measurements are typically performed at the macro-scale. This method has been successfully applied to several scenarios, including heterogeneous porous media, as shown by Penta et al. (2021), vascularised tumours, as in the manuscript by Shipley & Chapman (2010); , and biofilms, e.g. as illustrated by Dalwadi & King (2020).
Our macro-scale results comprise a double Darcy's system of partial differential equations (PDEs) describing fluid transport within and between compartments, and double advection-diffusion-reaction equations for both drug and heat transport. The influence of the micro-structure appears in the hydraulic conductivities, particle diffusion coefficients and thermal conductivities, which can be determined by solving appropriate periodic cell problems. The macro-scale system of PDEs is solved by finite elements in a spherical coordinate setting. The results elucidate the role of tortuosity and absorption rate, as well as their mutual interplay, on heat transport generated by nanoparticles in vascularised tumours. Nabil & Zunino (2016) discussed the hyperthermia cancer treatment using IONP by primarily focussing on the adhesion mechanisms (so that excited particles do not extravasate from the vessels to the tumour). However, this work is different from the one by Nabil & Zunino (2016) as here the tumour is modeled as a three-dimensional domain, which comprises the interstitial spaces and the vessels. The two domains are separated by an interface that represents the vessels' walls. Moreover, there are some differences between the two works related to the macro-scale geometry, computational technique and other concepts such as external boundary conditions. For instance, Nabil & Zunino (2016) represent the tumour at the macro-scale as a cube for the sake of simplifying the numerical computations. In contrast, we modelled the vascularised tumour as a sphere as an analogy with the works by ; ; Shipley & Chapman (2010). This geometry has significant implications on the fluid and drug transport profile and also this particular shape is very convenient when comparing results against experiments, as shown by Jain & Baxter (1988). In addition, the homogenisation technique adopted by Nabil & Zunino (2016) relies on the immersed boundary method, as illustrated by Cattaneo & Zunino (2014a,b). The vessels are dealt with as though they were one-dimensional lines, nevertheless carrying relevant three-dimensional information via appropriate singularities on the boundaries. Our new model retained the three-dimensional character of both the vessels and the interstitial spaces, and asymptotic homogenisation is being used to perform the upscaling and achieve computational feasibility. In this way, the geometrical differences between the vessels and tumour are smoothed out on the macroscale. Moreover, in the present work, we encode information related to the fine-scale structure of the individual compartments, as well as the transport that is occurring across the vessels' walls. The latter is reflected into appropriate sources at the macroscopic scale. In the work by Nabil & Zunino (2016), the interface is not resolved as the vessels are immersed in the three dimensional tumour. However, in their case, information concerning fluid, drug and heat transport across the interface is retained and it appears likewise as a source in the resulting macro-scale model. Furthermore, the vessels temperature is constant in their work, but it varies in our work, as the vessels are represented by a separate compartment in three dimensions, which is described by its own governing equations. In addition, Nabil & Zunino (2016) consider different time steps that depend on the size of nanoparticles. A 40 minutes time interval appropriate for very small nanoparticles and 12 h, 24 h and 48 h for large ones that are called vascular magnetic nanoparticles. In our case, we do not focus on the size of nanoparticles and we focus on a 4 days time interval (and highlight the dynamics which takes place during day 1 by means of 4 different time points at 6, 12, 18, and 24h). Finally, we have assumed micro-scale periodicity, which is a limitation of the present model. It allows us to deal with complicated and potentially tortuous microvessels that are often encountered when dealing with vascularised cancer, as shown by ; Jain et al. (2007). This paper is structured into different sections that are organized as follows. In Section 2 we describe the mathematical model by emphasizing the main assumptions and underlying physical phenomena. These include the differential equations for fluid flow, particle transport and heat convection in both the tumour vessels and the interstitial compartment. We also address the fluid, drug, and heat exchange taking place across the interface via setting up appropriate interface conditions. In Section 3, the differential equations are formulated in non-dimensional form. In Section 4, we apply the asymptotic homogenisation technique and derive the macro-scale results, which are then summarized in Section 5. In Section 6, we briefly discuss how the homogenised coefficients are determined on the basis of a microstructure. In Section 7, the differential equations describing particle and heat transport are written in spherical coordinates and supplemented by corresponding macro-scale initial and boundary conditions. In Section 8, the results obtained via numerical simulations are illustrated and discussed. In Section 9, concluding remarks are presented.

Mathematical modeling
In this work, we address mathematical modelling of cancer hyperthermia therapy carried out via nanoparticles delivery. We represent the vascularised tumor as a three-dimensional domain Ω ⊆ R 3 . The tumour tissue comprises two regions. The interstitium is denoted here by Ω t and the blood vessels' network by Ω v , such thatΩ t ∩Ω v =Ω.
The tumour system under consideration is multi-scale in nature and the typical distance between adjacent blood vessels d ≈ (50 − 100) μm is much smaller than the average size L ≈ 0.5 cm of the cancerous region, as reported by . Therefore, we define the small parameter We are interested in describing heat transport and the subsequent temperature distribution, which is driven by nanoparticles that are considered as being transported as passive scalars. Therefore, we can assume that delivery of nanoparticles is occurring via diffusion through both the vessels and the interstitium, advection due to the fluid flow in both compartments and extravasation across the vessels' walls. In the next section, we illustrate the governing equations for fluid, drug and heat transport in both the vessels and the interstitial space.

336
T. AL SARIRI AND R. PENTA All the variables in this model such as the pressure P, the concentration c, the velocity u and the temperature T are functions of both x and t.

Interstitial fluid transport
The tumour interstitial region is represented as an isotropic porous medium governed by Darcy's law: together with the incompressibility constraint where u t , κ and P t are the interstitial fluid velocity, conductivity and pressure, respectively.

Microvascular flow
The blood flowing in small capillaries is considered an incompressible viscous fluid and it is transported in the body through the vessels. Although non-Newtonian effects may become relevant in small capillaries, we assume that the blood is a Newtonian fluid as a first approximation, as done e.g. by Shipley & Chapman (2010). Therefore, the blood flow in the vessels can be described by the Stokes' problem, which, in absence of body forces, reads where u v , μ and P v are the fluid velocity, viscosity and pressure in the blood vessels, respectively.

Transport of particles
The concentration dynamics of the particles in the vessels c v can be described by the advectiondiffusion equation in Ω v . Absorption of particles in the interstitial compartment is represented by a linear uptake term. Therefore, the governing equations for the interstitial and vessels concentrations c t where D v , D t , and Λ are the particles' diffusivities in the vessels and in the interstitium and the uptake rate, respectively.

Heat convection in the tumour
The heat generated in the tumour is driven by application of a magnetic field that affects the nanoparticles injected into the bloodstream. When the magnetic field is removed, the magnetization of the particles returns to zero. The relaxation can take place either via the particle's rotation around the fluid (Brownian relaxation), or rotation of the magnetic moment around each particle (Neel relaxation), as explained e.g. by Pankhurst et al. (2003). These rotations generate the heat that provides an increase in temperature.
The temperatures in the vessel T v and in the interstitium T t are determined by a system of advectiondiffusion equations. The heat generated by the magnetic field is represented by the heat source αf (c t,v ) and, for the sake of simplicity, here we assume (2.8) The efficacy of the heat produced by the IONP depends on the absorption rate α (Nabil & Zunino, 2016). The latter parameter is in turn related to magnetic field properties (such as intensity and frequency) and particles' shape, although such features are not explicitly taken into account in the present work. The coupled system of PDEs then reads (2.10) Here, K t and K v are the thermal conductivities in the interstitium and in the vessels, respectively. The parameters ρ t and γ t are the tissue density and specific heat capacity, while ρ v , γ v are the blood density and the blood specific heat capacity.

Interface conditions
The interface between the two domains is denoted by Γ = ∂Ω v ∩ ∂Ω t . The fluid flow across the vessel is assumed to be continuous and depending on the pressures' difference between the two domains. We assume that the blood flux across the vessel wall is determined by Starling's law, i.e.
(2.11) 338 T. AL SARIRI AND R. PENTA The vector n represents the unit outward vector normal to the vessels' wall. The parameter L p represents the permeability of the blood vessels, which reflects the leakage of the vessels' wall. In order to close the problem we need to specify a condition for the tangent component of the fluid velocity to account for slip over the porous interface. We assume a Beavers and Joseph condition as done by Shipley & Chapman (2010); , i.e.
where ϕ is a non-dimensional parameter and τ denotes both of the unit vectors tangent to the vessels' walls. The parameter k is the tissue permeability, which is related to the hydraulic conductivity κ by the following relationship (2.13) The particles' flux is assumed to be continuous and proportional to the particles' concentration difference between the vessels and tumour: where p is the diffusive membrane permeability. In this work, transvascular advection across the vessels' membrane, as investigated for instance by Mascheroni & Penta (2017), is neglected for the sake of simplicity, although the theoretical derivation that follows could be readily extended to such contributions. The heat flux is likewise expressed in terms of the temperature difference between the two domains Ω t and Ω v as follows where β is the heat transfer coefficient.

Non-dimensional form of the model
In this section we perform a non-dimensional analysis of the system of PDEs (2.2)-(2.10) supplemented with interface conditions (2.11)-(2.12) and (2.14-2.16) as follows: The index z = v, t denotes either the vessels or the tumour, while C r , T, C, d and L are the reference concentration, temperature, pressure gradient, inter-capillary distance and average tumour size, respectively.
By dropping the primes for the sake of simplicity of notation, the dimensionless PDEs can be written as: with boundary conditions: where the primes have been dropped for the sake of simplicity. The non-dimensional numbers are defined as: Here,κ is the non-dimensional hydraulic conductivity. The coefficients Υ andᾱ are non-dimensional uptake rate and absorption rate. The numbersL,p,β andφ are the non-dimensional vessels' hydraulic and diffusive permeabilities, heat transfer and Beavers and Joseph coefficients, respectively. The nondimensional diffusivities of the particles in the vessels and the tumour are the reciprocal of their corresponding Peclet's numbers, i.e.D where (3.20) The non-dimensional thermal conductivities are given bȳ (3.21) The ε scaling appearing on the right hand side of interface conditions equations (3.11-3.16) is the appropriate one to ensure that blood, drug and heat fluxes inside the tumour stays finite in the limit ε → 0, as observed by . Furthermore, we assume that the parameters appearing in equations (3.17-3.21) are finite in the limit as ε approaches zero. This is done consistently with the approach carried out by  and ensures that both drug and thermal diffusivities, which are well known to play a crucial role in the nanoparticles' dynamics, are captured at leading order. There exist different scaling choices in the literature, see, e.g. the work by Shipley & Chapman (2010), where the authors perform the upscaling of the equations describing fluid and drug transport in vascularised tumours and their choice concerning distinguished limits of the Peclet's numbers results in a suite of reaction-advection models.

The asymptotic homogenisation method
The asymptotic homogenisation method is an upscaling strategy that provides a macro-scale description of a given physical system informed by the microstructure. It has been successfully applied to a large variety of real-world scenarios including previous investigations related to fluid and drug transport in biological tissues and vascularised tumours, as done, e.g. In our model, the application of the multi-scale method is motivated by large difference in sizes between the inter-vessels distance and the tumour radius, as assumed by Shipley & Chapman (2010); . In particular, we assume that these two scales are well separated, so that the small parameter ε defined in equation (2.1) is much smaller than 1. Therefore, we can define two formally independent variables y (the micro-scale) and x (the macro-scale) related by: Following the above-mentioned previous works, we enforce periodicity with respect to the microscale variable y, and we assume that every unknown field can be represented in power series of ε as: where υ collectively represents any variable described in our model namely P z , c z , u z or T z , (with z = v, t).
The differential operators transform according to the chain rule: We equate the same power of ε l , (l = 0, 1, 2, ...) to find suitable differential equations in order to close the problem for the leading order variables P v and T 0 t . As we would like to obtain a system defined on the macro-scale only, for fields that retain a dependence on the micro-scale variable y, we can integrate over the periodic cell. From now on, since we are assuming micro-scale periodicity, Ω v and Ω t represent the vessels and interstitium cell portions, respectively, and the micro-scale cell average is defined by: where |Ω v | and |Ω t | are the vessels and interstitial cell volumes' portions. In particular, we assume that all the fields υ (l) , l = 0, 1, ... are y-periodic.
The cell volume fraction |Ω v |(x), |Ω t |(x), and the surface area of the interface S(x) are defined by, In our case we enforce macroscopic uniformity such that S, |Ω t |,|Ω v | are constant. Analyses of nonmacroscopically uniform structures are beyond the scope of this work; however, they could be relevant in several contexts and alternative approaches to deal with such heterogeneities can be found in the works by Penta et al. (2014; Dalwadi & King (2020), and the references therein.

The macro-scale model obtained via asymptotic homogenisation
By applying the asymptotic homogenisation steps described in the previous section to the system of equations 3.23.9 with boundary condition (3.11)-(3.16) we obtain the macro-scale differential equations for the zero-th order pressures, velocities, concentrations and temperatures P v . These can be summarized as follows, while for the full derivation of the model, which is carried out by generalizing the methodology carried by  to heat transport equations with a concentration-dependent source, the reader can refer to Appendices A and B.
are effective hydraulic, diffusion and thermal conductivity tensors in the vessels' and interstitial compartments, respectively. Here, |Ω v | denotes the vessel volume, |Ω t | is the interstitial volume and S is the vessels' wall surface.
The system of equations (5.2) describes transport in a porous medium with mass transfer between compartments. The leakage of the blood across the vessels is reflected in the mass exchange between the two compartments, which is proportional to the difference between the leading order pressures.
The particles' transport in the vessels' and interstitial compartments depends on the fluid flow and it is represented by the system of coupled advection-diffusion-reaction equations (5.3). Similarly, the system of coupled advection-diffusion-reaction equations in (5.4) describes the heat transport at the macro-scale and the temperatures' profiles depend on both fluid and particles' transport.
Moreover, the macro-scale coefficients, namely hydraulic conductivity tensors, diffusion tensors and thermal conductivity tensors, can be determined by solving the cell problems.
The hydraulic conductivities Y v and the auxiliary tensor Y t are defined as where the tensor W and vector r are satisfied the cell problems: and, The diffusivity tensors F v and F t are given by: and while the auxiliary variables a and b satisfy the cell problems: The thermal conductivities N v , and N t are defined as: For instance, the cell problems related to interstitial fluid flow, drug transport and heat transport are to be closed by a further condition for uniqueness to be achieved (e.g. by assuming the null cell average of the auxiliary variables in the cell), as illustrated by Cioranescu & Donato (1999); .

The effective coefficients and micro-scale cell problems
In order to close the system of PDEs at the macro-scale, we need to compute the effective coefficients by solving appropriate cell problems at the micro-scale. The differential problems that are related to the hydraulic conductivity tensors are discussed by . The authors solved the differential problems numerically and investigated the influence of the vessels' tortuosity on the hydraulic conductivity tensors. Mascheroni & Penta (2017) extended the analysis carried out by  to compute the effective diffusion coefficients by solving the cell problems related to drug transport, i.e. finding the solution for the auxiliary variables that are called a, and b in the present manuscript, cf. (5.16)-(5.17), and (5.18)-(5.19). They also varied the geometrical tortuosity and found its impact on the tensors F v and F t . Changing the vessels' shape or tortuosity implies changes in the interstitial and vessels' volumes. Both the vessels' hydraulic conductivity and particles' diffusivity are affected by the vessels' tortuosity. In particular,  show that the hydraulic conductivity exhibits a nonlinear decreasing profile at increasing tortuosity, while Mascheroni & Penta (2017) show that also diffusion decreases as tortuosity increases, although to a lesser extent. In contrast, the interstitial coefficients are not significantly affected by micro-scale changes in the geometry under consideration.  solved the problem that corresponds to those related to the interstitial fluid, drug and thermal auxiliary variables r, b and e in our work, i.e. (5.12), (5.18)-(5.19) and (5.24)-(5.25), respectively. In particular, these latter problems (5.11)-(5.12), (5.18)-(5.19) and (5.24)-(5.25) are equivalent, so the auxiliary variables r, b and e solve the same problem and ∇ y r = ∇ y b = ∇ y e.
(6.1)  concluded that as long as the vessels' volume fraction is much smaller than the interstitial one, the influence of the micro-scale on the interstitial coefficients is negligible, i.e. they observed that ∇ y r t = ∇ y b t = ∇ y e t ≈ 0. (6.2) As such, from now on we focus on the micro-scale cell problems in the vessels' compartments and account for (6.2), so that, by recalling the definitions (5.6), (5.15) and (5.21), we can assume By following ; Mascheroni & Penta (2017), we enforce invariance with respect to the three orthogonal axes so that the auxiliary tensors in the vessels W, F v , N v are proportional to the identity tensor.
In particular, for the diffusivity F v and thermal conductivity N v , we have: where we can further notice that a and g are actually the solution to the exact same cell problem (5.16-5.17) or equivalently (5. 22-5.23). This leads to the solution of a standard Laplace problem, which reads, e.g. for the component g 1 : ∇g 1 · n = n 1 on Γ (6.5) supplemented by a further condition to ensure uniqueness, e.g.
The analysis that follows is carried out by varying the tortuosity of the microvessels according to . In particular, the center line of every vessel is defined as: where A is the amplitude, ω is the frequency, s is the local parametrization along the branch and l is branch length and we have 0 ≤ s ≤ l. The tumour interstitium is the domain that is complementary to vessels' compartment in the cubic cell. We exploit the solutions of the cell problem 6.4-6.5, which is solved by Mascheroni & Penta (2017), to investigate the role of tortuosity on the homogenised thermal conductivityÑ v by varying the amplitude and spatial frequency ω. The profile of the relative thermal conductivity (and diffusivity) which is based on the results reported by Mascheroni & Penta (2017), is shown in Fig. 2.

346
T. AL SARIRI AND R. PENTA

The mathematical model for a spherical tumour
We assume that the vascularised tumour can be represented in spherical coordinates with radius R. Also, we presume that the symmetric tumour is isolated and it interacts with surrounding environment through the vessels. Assuming radial symmetry, the model reads as follows: where 0 r ≤ R and 0 t T , where T is the time interval under investigation. The macro-scale system describing the fluid transport (5.1-5.2) was solved analytically when accounting for spherical symmetry by . In this case, the system to be solved reads (neglecting from now on the labelling indicating the leading order character of every field (0) for the sake of simplicity of notation): (7.8) The above system of equation was solved by  by accounting for boundary conditions that are consistent with Jain et al. (2007) and the references therein, i.e. those for an isolated tumour with fluid flow driven by the difference between the vascular and the interstitial pressures (the vascular pressure is actually considered constant in Jain et al. (2007) and the references therein). In the above, H v is the vessels' hydraulic conductivity parameter, which, according to , ranges from 2.20 · 10 −4 for a regular microvasculature to 4.89 · 10 −6 for the most tortuous scenario, and satisfies: and H t = 1 in our case by means of (6.2). The solutions of the system (7.3)-(7.8) derived by  are summarized below.

Initial and boundary conditions
We assume that no particle is present in the whole system at t = 0. Also, both the drug and the heat fluxes must vanish in the tumour centre as a consequence of the radial symmetry assumption. We assume a vessels' bolus injection with clearance time ς at the boundary of the macro-scale domain, which means that the concentration of the particles declines exponentially due to body elimination effects in the plasma, as shown by Chou et al. (2013). We also assume the continuity of particles' concentration at the boundary of the interstitial region. The initial temperatures are set to be the standard vessels' temperature 310.15K. Following the approach by Nabil & Zunino (2016), we impose Robin condition on the boundary of the tumour interstitium to account for the heat transfer between the tumour and the vessels' mediated by intermediate layers of tissue (Nabil & Zunino, 2016;Cervadoro et al., 2013;Saeedi et al., 2017;Kaddi et al., 2013). The initial and boundary conditions can be summarized as follows (7.14) and, (7.15) The finite element software Comsol Multiphysics is used to solve the model and the values of the parameters are provided in Table 1.
We have used the finite element commercial software COMSOL Multiphysics, version 4.3a, and both the drug and the heat transport systems (7.1)-(7.2) have been implemented by means of the convection-diffusion module in coefficient form equipped with boundary and initial conditions (7.14)-(7.15) and parameters taken from Table 1. The spatial discretization is carried out by means of P2 elements, while for the discretization in time an implicit backward differentiation formula method is embraced, similarly to Mascheroni & Penta (2017). Although the system is solved in non-dimensional form, the temperatures and the absorption rate are shown in dimensional form in the plots to foster the Reader's clarity in terms of comparison against the previous literature.

Results and discussion
Mascheroni & Penta (2017) studied the macromolecules distribution in both the vessels and the interstitial compartment using the advection-diffusion-reaction equations derived by . The reaction terms are related to the uptake of anti-cancer agents, as well as additional contributions due to the upscaling of transvascular diffusion of particles. The authors presented the result for a spherical tumour, and they discussed the impact of tortuosity on drug transport.
In the present work, we extend the works by ; Mascheroni & Penta (2017) to heat transport and solve the resulting systems of PDEs to obtain the temperature maps that are driven by nanoparticles' transport in the context of cancer hyperthermia. Although also the drug transport analysis carried out here differs from the one by Mascheroni & Penta (2017) in terms of the choice of   here, while zero diffusive interstitial drug flux is assumed by Mascheroni & Penta (2017)), a qualitative comparison concerning the drug transport problem is still possible, and provides a benchmark supporting the reliability of the results presented here. We commence by first presenting our results concerning the solution of the drug transport problem and then show the results concerning temperatures maps against the relative radial position at varying microvessels' tortuosity and absorption rate.  The main results show that geometrical tortuosity can significantly impair heat transport within the tumour and that a higher magnetic field can be required to reach a temperature that is sufficiently high to kill tumour cells by cancer hyperthermia. We provide a detailed and more quantitative description of the results below.

Particle transport
The results displayed in Fig. 3 and Fig. 4 are presented in terms of the leading order concentrations in the tumour and the vessels against the non-dimensional radius within a chosen period of time of 24 hours and 96 hours, respectively.
The solution clearly shows that the nanoparticles manage to reach the tumour center and both concentration profiles are very similar at leading order. This is due to the particles' exchange between the two compartments through the vessels' wall combined with continuity of concentrations on the tumour boundary. Due to the assumption of drug delivered intravascularly via a bolus injection on the boundary, the particles' concentration decreases steadily after 6h (432) at r = R and reaches zero after two days which cause the over all decline on the concentration. However, a fraction of the initial concentration is still able to reach the tumour centre by the end of the time interval under investigation. In the period of time (24h − 72h), the particles' concentration in the center increases from 1% to approximately 7% of the initial concentration. After that, the concentration in the center starts to decrease slightly, i.e. in the period (72h − 96h). In addition, the concentration in the past two days reaches a plateau when moving towards the center. Nabil & Zunino (2016) presented their result in a cubic symmetric setting and they found that the particles' concentration decreases with time. Moreover, the concentration of nanoparticles in the vessels becomes almost uniform at the end of the circulation time they investigate, which is 48 hours.
Nanoparticles and in general drugs are eventually metabolized by tissue. This is done at a specific rate, also referred to as the uptake rate, which depends on the properties of the tissue and drugs at hand, as discussed by Tchoryk et al. (2019). Mascheroni & Penta (2017) compared two specific macromolecules characterized by different uptake rates, with order of magnitudes varying from 10 −11 s −1 to 10 −5 s −1 , as also mentioned by Weinberg et al. (2007). In Fig. 5 we show the influence of high tissue uptake rate on the particle distribution in vascularised tumours, and we then increase the value of the uptake rate from 1.07 · 10 −11 s −1 (see Table 1) to 10 −5 s −1 .
The concentrations in both compartments are decreasing and approximately approaching zero at the center in all periods of time. High uptake rate leads to fast washing out the particles and few of them can reach the center of the tumour. The particles in this case are metabolized very fast by the tumour before they are transported into the tumor center.
The concentrations profiles are qualitatively in agreement with Mascheroni & Penta (2017) and this is shown for the case of the most tortuous vessels' network considered by  and Mascheroni & Penta (2017), i.e. ω = 3 and A = r c , see also Table 2 and Fig. 8.

Heat transport
We now present the major results obtained by solving the full system of macro-scale coupled PDEs (7.1)-(7.2) by finite elements. The tortuosity of the microstructure is varied according to the Table 2 The computational result for the non-dimensional vessels' thermal conductivity with different vessels tortuosity from the analysis by  and Mascheroni & Penta (2017) ∂y 1 v 0 0 8.1· 10 −2 6.149 2.30 2.20 · 10 −4 3.6 · 10 −1 1 r c 7.6 · 10 −2 6.154 2.32 1.69 · 10 −4 3.19 · 10 −1 2 r c 6.9 · 10 −2 6.162 2.57 6.24 · 10 −5 2.13 · 10 −1 3 0.75r c 6.8 · 10 −2 6.162 2.82 2.02 · 10 −5 1.59 · 10 −1 3 r c 6.5 · 10 −2 6.165 3.25 4.89 · 10 −6 0.9 · 10 −1 values reported by  and Mascheroni & Penta (2017) corresponding to five representative geometries, and the two extreme cases are shown in Fig. 8. We have observed that the temperature increases and reaches its maximum after one day, then starts to decline and the maximum temperature varies with vessels tortuosity. As we have also remarked in the Introduction, increasing the tortuosity reduces fluid and particles convection within the tumour, as shown by  and Mascheroni & Penta (2017). As such, this leads in turn to impaired heat convection driving a decline in temperatures. Therefore, the more regular the vessels, the lower magnetic field intensity (which is here encoded in the absorption rate coefficient) is needed to reach the desired target temperature. The plots showing the vessels' and interstitial temperature maps are shown in Figs 6 and 7 at different times, for the first 24h and from day 1 (24h) to day 4 (96h), respectively. These results are related to the most tortuous (i.e. corresponding to the case ω = 3, A = r c reported in Table 2) vessels' microvasculature considered by , see Fig. 8. Figure 6 clearly shows that the temperature increases with time as it reaches its maximum after 24 h. It then starts to decline steadily with time, because the concentration in the blood decreases exponentially according to the bolus injection, cf. initial condition (7.14).
Also, for all period of times under investigation, the temperature in the center is higher than the boundary, as the particles are transported towards the center. This can be explained by the fact that heat trasport is driven by a significant diffusive component as opposed to drug transport, which is mostly Fig. 7. Temperature maps in both the vessels (a) and the interstitium (b) vs radius-(low uptake rate) from day 1 (24h) to day 4 (96h). Fig. 8. The most tortuous micro-vasculature (on the right) vs the regular one (on the left) as setup by , where further details can be found. This Figure (minimally adapted to serve this work) is reprinted from . The role of the microvascular tortuosity in tumor transport phenomena, Journal of Theoretical Biology, 364, page 86, Fig. 4 (a) and (d), with permission from Elsevier. driven by convection instead (cf. thermal conductivities K t or K v vs the particles diffusion coefficients D v or D t in Table 1. In fact, the non-dimensional diffusion coefficients as defined in (3.19) are of order ≈ 10 −4 to 10 −5 as opposed to the non-dimensional thermal conductivities that are of the order of ≈ 10 −1 to unity). This explains the difference between the drug concentration and temperature profiles, despite both phenomena being governed by formally a similar set of advection-diffusion-reaction equations. In fact, the role of advection is more prominent in driving drug transport rather than heat transport, as it can also be observed by the more localized concentration peaks (cf. Fig. 4), as opposed to the smoother and more uniform heat transport process, which is reflected in the temperature profiles as per Figs 6 and 7. At 24 h the temperature in the center is approximately 313K (39.8 • C) where in the boundary it reaches the blood temperature (310.15K), which is prescribed via the boundary condition. This is also in agreement with temperature profile previously reported in other works that address this problem using different modelling frameworks such as those by Bagaria & Johnson (2005); Golneshan & Lahonian (2011);Dutz & Hergt (2014).
The distribution of heat in the tumour is in agreement with Nabil & Zunino (2016) as they reported that the temperature in the center of the cube is higher than in the edges. However, Nabil & Zunino (2016) found that the temperature increases with time (48 h). This discrepancy is related to our different set of boundary conditions. In our case we have an exponential decrease in the particles' concentration, which is directly proportional to the heat source related to the magnetic absorption rate, thus eventually causing a temperature decline over time.

The influence of absorption rate and vessels' tortuosity on the heat distribution
The previous analyses in section 8.2 are related to tumour microvessels, which are most tortuous and leaking vessels with (ω = 3, A = r) as opposed to the healthy ones (ω = 0), see e.g. the works by ; Shipley & Chapman (2010) and Carmeliet & Jain (2000).
The structure of the vessels and their tortuosities are not uniform and they vary from one point to another in the tumour mass, as described by . ; Mascheroni & Penta (2017) discussed the impact of the vessels' geometry on fluid and drug transport, respectively. They deduced that the vessels' tortuosity leads to a relevant decrease in both hydraulic and diffusivity properties of the vessels thus impairing fluid and drug convection within the tumour. Here, we perform a parametric analysis by varying the tortuosity of the vessels' micro-structure and capture its effect on the temperature maps. We make use of the setting that has been exploited by . The data associated with the various parameters involved are reported by ; Mascheroni & Penta (2017). The results show that heat transport is impaired at increasing vessels' tortuosity from the most regular vessels characterized by ω = 0 (representing healthy vessels) to the most tortuous vessels (representing tumour vessels at an advanced stage) with ω = 3, A = r. The temperature varies between approximately 39.8 • C and 40.9 • C as we improve the regularity of the vessels. Also, we have observed (see Fig. 9) that the temperature decreases more remarkably from the vasculature corresponding to ω = 3, A = 0.75r and the most tortuous one (ω = 3, A = r). This is ultimately related to impaired drug and fluid transport, and the latter (especially fluid convection) decreases sharply when the micro-scale fluid profile is no longer parabolic, as shown in Fig. 8 and discussed by .
As the particles are transported smoothly in the healthy vessels and the concentration is high even in the second day, the temperature reaches its maximum 41.5 • C after two days, see Fig. 10. The difference in maximum temperatures between the regular vessels and most tortuous ones is approximately 1.5 degrees.
The temperatures achieved with different tortuosities are very close to the medical and experimental results that show that 42 • C is the appropriate temperature for hyperthermia treatment and 43 • C-44 • C for magnetic hyperthermia treatment, see also the works by Laurent et al. (2011);Silva et al. (2011);Ling-Yun et al. (2013).
Furthermore, the absorption rate (which is proportional to the magnetic field intensity) plays important role on the heat distribution as it mediates the temperature increase that is caused by the nanoparticles' concentration. Therefore, we have varied the value of absorption rate α at increasing tortuosity to detect the impact of these variations on temperature maps. The absorption rate of magnetic nanoparticles is proportional to the square of the magnetic field intensity, as well as its frequency. It  also varies with respect to the nanoparticles size and material, and a range of variation of four orders of magnitude, i.e. from 10 3 W/Kg to 10 7 W/Kg, has been reported by Cervadoro et al. (2013).
We have observed a linear relationship between the absorption rate and the heat distribution for each geometry under investigation, see Fig. 11.
Compared to the values of hydraulic conductivity, diffusivity and thermal conductivities for different vessels tortuosity, we observed that the temperature increases by 6% in the most regular vessels when the absorption rate increases by one order of magnitude. However, the temperature increases by 4.6% for the most tortuous vessels when the same change of absorption rate is applied. Moreover, when the value of the absorption rate is 6 · 10 6 W/Kg, the temperature difference between the tortuous and regular vessels is almost one degree, while the difference is approximately two degrees when the value of absorption rate is 10 7 W/kg.
The absorption rate of nanoparticles can be varied in practice during experiments in order to have the suitable temperature that is required to kill the cancer cells.

356
T. AL SARIRI AND R. PENTA Fig. 11. Temperature maps vs normalized radius at increasing absorption rate for different microvascular geometries after 24 h.

Conclusion
We have derived a new mathematical model which describes the heat transport occurring in vascularised tumours due to IONPs delivered intravascularly, as per current cancer hyperthermia protocols.
We have obtained the results by means of the asymptotic homogenisation technique in terms of a tissue-scale macroscopic description of the coupling between fluid, particles and heat transport, as well as their exchange across the vessels' membranes.
The new coupled system comprises six PDEs describing both interstitial and vascular pressures, concentrations of nanoparticles, and temperatures.
A double Darcy's system describes fluid flow, while the concentration of nanoparticles and heat transport are both governed by double advection-diffusion-reaction system of PDEs.
The impact of the micro-structure is reflected in the effective tensors of coefficients representing the hydraulic and thermal conductivities, as well as particles' diffusion. These effective coefficients can be computed by solving periodic cell problems where the geometry of the micro-vessels is clearly resolved. The role of transvascular mass, heat and particles transport and uptake appears in suitable macro-scale exchange terms that provide the coupling between the governing equations in the vessels and the tumour.
We have solved the full model by means of finite elements, and we have observed that vessels' tortuosity can impair heat transport within the tumour mass, so that regularization of the micro-vessels can produce a significant (1-2 degrees) increase in the maximum temperature that is reached in the tumour center under the same therapeutic conditions (which are here reflected in the tumour absorption rate, which is in turn related to the magnetic field and nanoparticles' properties). Furthermore, we have investigated the impact of a change in the absorption rate for different micro-vessels' geometries, and this analysis can pave the way for informed cancer hyperthermia parameters depending on the geometry of the microvessels, which is ultimately related to the tumour stage. For example, the heat distribution with absorption rate 10 7 fluctuates between approximately 43.5 • C and 41.5 • C, which is aligned with the required temperature to destroy cancer cells, as mentioned in section 8.2.
This analysis is open for improvement and further developments. We have chosen to present the results by means of a spherical coordinate setting as this has enabled us to deduce the semianalytical results that can be readily compared against the current literature. However, our finite element computational platform can be generalized to generic macro-scale geometries depending on the actual tumour shape at hand. This work could also be generalized to include nonlinear heat sources and nonlinear drug uptake, as, given the current scaling assumptions, relevant modifications would only appear at leading order. Different boundary and initial conditions could also be taken into account depending on the interplay between the tumour mass and the surrounding and on the clinical injections conditions at hand. In addition, we have derived our model by considering the same distinguished limit as in ; Mascheroni & Penta (2017) in terms of Peclet numbers and nondimensional thermal hydraulic conductivities. Alternative distinguished limits, which would result in purely convective heat and drug transport contributions at leading order could be considered (see, e.g. Shipley & Chapman (2010) when these are investigated for macromolecules transport). An interesting further development of this work also resides in a comprehensive analysis of admissible distinguished limits that exist for this system, with particular reference to particles' uptake, and diffusion phenomena occurring in different regions of the domain under consideration, see, e.g. Dalwadi et al. (2018); Ptashnyk & Roose (2010), respectively. Furthermore, we did not consider the adhesion between the particles and the vessels' wall, and also inter-particle cohesion, which are e.g. discussed by Decuzzi & Ferrari (2006); Nabil & Zunino (2016). Whenever the adhesion mechanism is primarily driving heat transport (i.e. when the nanoparticles' dimensions prevent them from extravasating into the tumour), the dynamics of the problem changes, but can be as well investigated by means of upscaling techniques such as the asymptotic homogenisation. As such, our results can indeed also provide a basis for future research concerning different transport mechanisms, which will be investigated in future works.
Moreover, in this work, we have defined the transformation of the magnetic field to heat mediated by the absorption rate. The latter is actually related to the properties of the magnetic field, such as its intensity and frequency, as well as the shape and dimensions of the injected nanoparticles, as per the Brownian-Neel relaxation formula (Kaddi et al., 2013;Pankhurst et al., 2003).
Finally, a natural generalization of this work resides in considering the role of tumour growth and deformations (see, e.g. the theoretical works by Penta et al. (2014); Penta & Merodio (2017)), and the corresponding stresses that are due to this interplay, on nanoparticle delivery.

360
T. AL SARIRI AND R. PENTA

A. Appendix: the upscaled governing equations for the vessels
The multi-scale differential equations governing vessels fluid flow, particle and heat transport in the vessels can be obtained from equations 3.23.3, (3.6) and (3.8) with interface conditions (3.11), (3.10), (3.13) and (3.15). We have, by enforcing equation (4.2) and multiplying each equation by a suitable power of ε: The interface conditions are: We now equate the same coefficients for ascending powers of ε. For ε 0 we obtain: Now, we determine the macro-scale relationships for the leading order velocity and pressure u v v and P (0) v . Equation (A.9) implies: This means that the leading order pressure in the vessels is y-constant. Equation (A.17) from the ε 1 conditions, together with (A.10), (A.11) and (A.12) from the ε 0 conditions, we obtain a Stokes' type problem for u  Integrating (A.30) over Ω v leads to the average leading order velocity in the vessels: Here, Equation (A.36) shows that the vessels' fluid flow obeys the Darcy's law with hydraulic conductivity given by relationship (A.37). In order to find the equation for the leading order pressure leading term P v , we take the average of (A.18) and make use of interface condition (A.19), as well as the divergence theorem with respect to y, as follows: Therefore, Thus by means of (A.36), v . The leading order concentration c (0) v can be found by using first equations (A.13) and (A.14), from which we deduce that the zero-th order concentration in the vessels depends on the macro-scale only, i.e.
We can formulate an ansatz for the solution c Integrating (A.25) and using the divergence theorem with respect to y, and subsequently making use of interface condition (A.26) from equating the same power of ε 2 yields: where the additional contribution over the boundary ∂Ω v \ Γ vanishes due to y-periodicity.
Using the ansatz (A.40), we obtain: is the effective diffusivity tensor in the vessels. Equation (A.43) is an advection-diffusion-reaction equation for c (0) v and it describes the macro-scale drug dynamics in the vessels. A macro-scale equation for the heat transport in the vessels can be obtained by following the same steps described above for particle transport. The solution of (A.15) and (A.16) is: Therefore, T v is y-constant. The solution of the problem obtained by collecting (A.23) and (A.24) from the ε 1 conditions can be formulated in terms of the following ansatz for T (1) v : where the auxiliary vector g is the solution of the following cell problem