Magnetic nanodrug delivery in non-Newtonian blood flows

With the goal of determining strategies to maximise drug delivery to a specific site in the body, we developed a mathematical model for the transport of drug nanocarriers (nanoparticles) in the bloodstream under the influence of an external magnetic field. Under the assumption of long (compared to the radius) blood vessels the Navier-Stokes equations are reduced, to a simpler model consistently with lubrication theory. Under these assumptions, analytical results are compared for Newtonian, power-law, Carreau and Ellis fluids, and these clearly demonstrate the importance of shear thinning effects when modelling blood flow. Incorporating nanoparticles and a magnetic field to the model we develop a numerical scheme and study the particle motion for different field strengths. We demonstrate the importance of the non-Newtonian behaviour: for the flow regimes investigated in this work, consistent with those in blood micro vessels, we find that the field strength needed to absorb a certain amount of particles in a non-Newtonian fluid has to be larger than the one needed in a Newtonian fluid. Specifically, for one case examined, a two times larger magnetic force had to be applied in the Ellis fluid than in the Newtonian fluid for the same number of particles to be absorbed through the vessel wall. Consequently, models based on a Newtonian fluid can drastically overestimate the effect of a magnetic field. Finally, we evaluate the particle concentration at the vessel wall and compute the evolution of the particle flux through the wall for different permeability values, as that is important when assessing the efficacy of drug delivery applications. The insights from our work bring us a step closer to successfully transferring magnetic nanoparticle drug delivery to the clinic.


Introduction
Currently, the main approaches in cancer therapy include surgery, chemotherapy, radiotherapy, and hormone therapy [1].The first is invasive while the latter three are non-specific.Hence, their efficacy is not only low but sometimes the cure can be more aggressive than the disease itself.In the late 1970s, magnetic drug targeting, which consists of injecting and steering magnetic drug carriers through the vessel system towards disease locations using an external magnetic field, was proposed as an alternative and more efficient treatment for tumors (see [2] and references therein).At present, there are many promising studies, both in vivo and in vitro [3,4,5] but, to our knowledge, only a few successful trials on human patients have been carried out [6,7,8].All studies demonstrate that magnetic forces can attract particles in the region near the magnet but there is a lack of knowledge on how to quantify and optimise the accumulation of particles [9].
Describing the movement of particles subject to a magnetic field in the bloodstream is a relatively difficult task due to the interplay of magnetic and hydrodynamic forces acting on the particles and the inherent difficulty of solving the Navier-Stokes equations.One further difficulty is to accurately model and simulate the non-Newtonian behaviour of blood.Experiments show that for low shear rates (∼1 s −1 ), the viscosity can be as high as 10-11 mPa•s while for shear rates in excess of 1000 s −1 , it tends to an asymptotic value of 3-4 mPa•s [10,11].In order to simplify calculations, previous studies consider blood as a Newtonian fluid [12,13,14,15] or as a non-Newtonian power-law fluid [16,17,18,19].Grief and Richardson [12] and then Richardson et al. [13] developed a continuum model for the motion of particles subject to a magnetic field, both using a Newtonian flow model for blood.In [12] it was shown, via a simple network model, that it is impossible to specifically target interior regions of the body with an external magnetic field; the magnet can be used only for targets close to the surface.In [13] the boundary layer structure in which particles concentrate was analysed.Using a similar approach, Nacev et al. [16] simulated particle behaviour under the influence of a magnetic field using a power-law assumption for the blood flow.Cherry and Eaton [17] developed a comprehensive continuum model for the motion of micro-scale particles using a value for the viscosity determined from fitting experimental data [20] to model blood shear thinning, and used the model to investigate magnetic particle steering through a branching vessel.
A structural difference within the models developed in the literature lies in how particles are described.Instead of describing particle distribution in blood as a continuum and using a diffusion equation to study the evolution of the particle concentration, several studies consider particles as discrete elements and track the trajectory of each individual particle [14,15,18,19,21,22].Yue et al. [14] implemented a stochastic model for the transport of nanoparticle clusters in a Hagen-Poiseuille flow in order to find an optimal injection point.Using a similar approach, Rukshin et al. [19] developed a stochastic model to simulate the behaviour of magnetic particles in small vessels.Lunnoo and Puangmali [18] used a generalized power-law model to investigate the parameters which play a crucial role in magnetic drug targeting, showing how difficult it can be to keep small particles in the desired region.Recently, in Boghi et al. [15] a numerical simulation of drug delivery in the blood in the coeliac trunk was performed.
In the present work we analyse the forces and parameters involved in the process of magnetic drug delivery and highlight the importance of considering realistic non-Newtonian models for blood flow and the motion of the nanoparticles in it.It is crucial to consider the non-Newtonian behaviour of the blood in order to predict if particles are able to reach the desired area.A mathematical model in a two-dimensional channel is introduced in Section 2. In Section 2.1, we explain how choosing an oversimplified constitutive law for the fluid in the vessel or an incorrect value for the viscosity in the centre of the channel can lead to significant errors.In Section 3 we demonstrate how magnetic forces act on particles depending on their size.Numerical simulations illustrating the influence of key parameters are presented in Section 4. We draw our conclusions in Section 5.

Governing equations for non-Newtonian blood flow
As illustrated in Figure 1, we approximate the vessel as a long and thin rectangular channel, consistent with the size of vessels in the human body (see Table 1).The particles are injected at the inlet of the vessel and their motion is driven by the field generated by an external magnet located at the bottom of the domain and by the drag force on the particles due to the blood flow.
Blood containing nanoparticles can be considered as a nanofluid.Typically, its motion is governed by the Navier-Stokes equations coupled to an advection-diffusion equation for the concentration of particles in the fluid [24].The Navier-Stokes equations describing the flow of an incompressible nanofluid in a vessel subject to an external magnetic field are where u F is the fluid velocity, ρ is the fluid density, p is pressure, τ is the extra-stress tensor and F mag the magnetic force acting on the fluid.
In this work, we consider the nanofluid, composed of blood and nanoparticles, to be a dilute mixture.This allows us to make two assumptions which substantially reduce the complexity of (1)- (2).Firstly, the density of the nanofluid, ρ, which typically depends strongly on the particle concentration [25,24], can be considered approximately constant and equal to the density of blood.Table 1: Typical values for the various types of vessels in human body: average length (L) and radius (R), estimation of the average number of a specific vessel in the circulatory system (N ).All values are adapted from [23].
Secondly, given that blood has negligible magnetization [26] and its overall character is found to be paramagnetic [27], the magnetic field only acts on the nanoparticles.Since we are considering a dilute mixture, the capacity of nanoparticles to modify the flow rheology is negligible and therefore we can assume F mag ≈ 0 in (1)-( 2).Due to these assumptions, the flow equations ( 1)-( 2) are effectively decoupled from the mass conservation (advection-diffusion) equation that will describe the concentration of particles within the vessel (see Section 3).

Blood as a non-Newtonian fluid
Bodily fluids, such as blood, saliva, and eye fluid are invariably non-Newtonian.In particular, blood is a concentrated suspension of particles in plasma, which is mainly made of water.The three most important 'particles' that constitute blood are: red blood cells (RBCs), white cells and platelets.RBCs, which are the most numerous, are mainly responsible for the mechanical properties of blood [23].In fact, their tendency to form (and then break down) three-dimensional microstructures at low shear rates and to align to the flow at high shear rates cause the blood's shear thinning behaviour, characterized by the monotonic decrease of the viscosity that tends to some limit for very high shear rates [11].In the case of blood, the structures lead to significant changes in its rheological properties and several models have been developed during the past 50 years in order to capture the complexity of this behaviour (some examples can be found in [28,29,30,31,32,33]).However, none of those models has been universally accepted.
In mathematical terms, we define a fluid as non-Newtonian if the extra-stress tensor cannot be expressed as a linear function of the components of the velocity gradient.The more general relation between the stress and rate-of-strain tensors can be written as where η( γ) is the viscosity, γ = ∇u + ∇u T and γ =   1 2 is the generalised shear rate.When η( γ) is constant the Newtonian model is recovered.A purely shear-thinning fluid will exhibit a monotonic decrease in viscosity with increasing shear rate.Practical fluids are more likely to exhibit a constant viscosity beyond certain high and low values of shear rate called 'Newtonian plateaus'.In between, a nonlinear viscosity relation links the two plateaus.
The aim of this section is to compare different types of simple, non-Newtonian fluid models.In particular, we will compare the Newtonian model with the power-law model, the Carreau model and the Ellis model.In Tables 2 and Table 3 we summarize the expressions for the viscosity and shear stress for each model and the corresponding parameter values, respectively.Note that under the assumption of flow in a long thin channel, we will significantly simplify the governing equations (Section 2.2).The shear stress relations in Table 2 are consistent with this reduction.In the power-law model m is constant.If n p < 1 the fluid is pseudoplastic or shear thinning and if n p > 1 it is dilatant or shear thickening; for blood n p = 0.357 < 1.If n p = 1 we retrieve the Newtonian expression.In the Carreau model λ is a constant and η 0 and η ∞ are the limiting viscosities at low and high shear rates, respectively.In the Ellis model η 0 is the viscosity at zero shear and τ 1/2 is the shear stress at which the viscosity is η 0 /2.The latter model does not include a high viscosity plateau, however for most practical situations such high strain rates are never reached and so this plateau value is not relevant.In the case of standard blood flow we do not anticipate high shear rates.
In Figure 2 we compare the behaviour of the four models, using the parameter values in Table 3.For moderate shear rates it is clear that the Ellis and Carreau models show good agreement but differences emerge as the shear rate increases above γ ≈ 1.6 s −1 .In attempting to compensate for the plateau values the power-law model seems to be inaccurate over the whole range of shear rates.
In the next section using the above models we determine the velocity profiles and demonstrate the importance of choosing the right model for blood transporting nanoparticles under the influence of an external magnetic field.

Nondimensionalisation and simplification of flow equations
Expressing equations ( 1)-( 2) in components, the governing equations for the fluid become which are subject to the no-slip conditions and to the symmetry condition Figure 2: The viscosity/shear rate plot in logarithmic scale for the power-law (red dotted-dashed line), the Carreau (black solid line) and the Ellis models (green dashed line).The light grey dashed lines represent the limiting viscosities η 0 and η ∞ .
In order to determine the order of magnitude of the terms in equations ( 5)-( 7), we proceed to non-dimensionalise.We define the non-dimensional variables In (5) we can balance the terms by setting V = εU , where ε = 2R/L 1 (since the blood vessels are long and thin).Hence, the non-dimensional continuity equation becomes Choosing P = LT /2R and assuming steady-state, we can write Assuming ρεU 2 T 1 and neglecting terms of O(ε) or smaller in ( 12)-( 13), the resulting governing equations, in dimensional form, are which is a much simpler system of equations, consistent with lubrication theory [36].

Comparison of four viscosity models
To determine the most accurate model describing blood flow, we now compare the fluid velocity profiles obtained for the four viscosity laws summarised in Table 2.In all cases we consider a Hagen-Poiseuille flow in a channel, which is driven by a constant pressure gradient, ∆p/L, and has a constant volumetric flow rate, Q.The results quoted below follow from the models in [35] for flows driven by a pressure gradient and a moving bottom surface, where we set the velocity of the surface to zero.To switch to the current vertical coordinate y from that used in [35] we set z = (y + R)h/(2R).The position of the turning point in the flow, z = z m = h/2, then corresponds to y = 0.
The Newtonian model: The solution for the Newtonian fluid is obtained by solving equations ( 14)-( 16) where the shear stress is specified in Table 2. From equation ( 16), the pressure does not vary with y.Hence, integrating (15) twice with respect to y and using the boundary conditions ( 8)-( 9) we obtain the well-known parabolic profile where the relation between ∂p/∂x and Q is determined using that [35].
The power-law model: Proceeding similarly to the Newtonian case, the velocity for a non-Newtonian power-law fluid with the corresponding shear stress from Table 2 becomes The main disadvantage of this model is that the viscosity η ∝ 1/(∂u/∂y) → ∞ since ∂u/∂y = 0 at y = 0 which is unrealistic.
The Carreau model: Choosing the Carreau model for the viscosity (see Table 2), the equations of the flow ( 14) and ( 16) are coupled with (15), which takes the form This expression cannot be integrated analytically for u(y), so we solve it numerically via the in-built bvp5c function in Matlab.
The Ellis model: Choosing the Ellis model from Table 2, the velocity of the fluid can be expressed analytically as In order to compare the four models and their velocity fields we take typical parameter values from the literature (Table 3).The viscosity depends on temperature, so for the whole paper we use parameters consistent with a typical body temperature, 37ºC.We assume a micro vessel with a radius R = 20 µm and the volume flux is Q = 8 × 10 −13 m 3 /s, consistent with average velocity measurements in capillaries or small arterioles [37].For our simulations, we will use a reference length several times larger than the radius, L = 500 µm.
The velocity and corresponding viscosity profiles for blood, obtained from the four different models, are compared in Figures 3(a) and 3(b), respectively.We observe that although the difference in velocity profiles is small for all models, the corresponding viscosities can differ significantly.The Ellis and Carreau models show good agreement in both velocity and viscosity.In the literature, the Carreau model is generally preferred due to its ability to predict both Newtonian plateaus.However for in vivo blood flow it is unlikely that the high shear plateau will be reached suggesting that the Ellis model provides an adequate alternative approximation.This is an important observation, since the Ellis model yields an analytical expression for the flow.This not only provides a  better understanding of the factors affecting the flow but also considerably simplifies the numerical study.In Section 3.1 we will verify that the Carreau and Ellis models provide similar results for the nanoparticle distribution.The Newtonian model has a constant viscosity and the power-law model cannot predict the Newtonian plateaus.Given that the motion of the magnetic particles is dependent on the viscosity of the fluid it travels through it is clear that both Newtonian and power-law models are not appropriate for modelling blood.

Nanoparticle motion in a non-Newtonian flow subject to an external magnetic field
The behaviour of the concentration of magnetic nanoparticles in the bloodstream is obtained following the continuum model developed by Grief and Richardson [12].The governing equation describing the motion of magnetic particles in the blood stream is an advection-diffusion equation for the particle concentration c(x, y, t): where u F is the fluid velocity (as discussed in the previous section), u p is the particle velocity and The diffusion is due to Brownian motion and shear-induced diffusion and hence D = D Br + D sh .Shear-induced diffusion arises due to the fact that the RBCs suspended in plasma collide with each other causing random motion with a diffusive character.As recently demonstrated by Liu and coworkers [38,39] via a lattice-Boltzmann multiscale simulations, the diffusion due to the Brownian motion in the case of nanoparticle transport in a small vessel can be important and, in some cases, even the predominant diffusion process.In particular, for nanoparticles smaller than 100 nm they found that D Br /D sh 1, that is Brownian diffusion is dominant for particles in this size range [39].Using the Stokes-Einstein equation for the diffusion of spherical particles through a shear thinning fluid [11], we can write the Brownian diffusion coefficient as where k B is the Boltzmann constant, T is the absolute temperature and a is the particle radius.
On the other hand, the shear-induced diffusion contribution can be approximated by where K sh is a dimensionless coefficient that depends on the blood cell concentrations and r RBC is the red blood cell radius.The coefficient K sh is difficult to measure but the value used in Table 4 is considered representative in the literature [12].Hence, The particle velocity is found by balancing hydrodynamic and magnetic forces.Using the definition of the Stokes drag, representing the hydrodynamic force on a spherical particle of radius a moving through a viscous fluid, we have The particles reach equilibrium velocity when F St balances the magnetic force F mag and this leads us to the expression In order to exert a magnetic force on magnetic nanoparticles in the vessels, a magnetic field gradient is required at a distance.We can define the magnetic force F mag on a single particle in a magnetic field B as Nacev et al. [16] have shown that for a magnet held at a long distance compared to the width of the vessel, we can assume the magnetic force is approximately constant in the vertical direction which avoids the need to solve Maxwell's equations.Richardson et al. [13], for example, use a constant value for the magnetic force given by where Υ is the magnetite volume fraction and B g is the gradient magnetic field.This is what we assume in this work.We consider a constant magnetic force acting perpendicular to the flow, i.e.F mag = F 0 j.Therefore, in combination with (26) we obtain that u p = (0, v p (y)) where v p = −F 0 /(6πa η( γ)).Finally, since the vertical fluid velocity is negligible, u tot = (u F (y), v p (y)) and D = D(y), the concentration equation ( 21) becomes We assume that the particles can flow out of the vessel through the walls with a certain vascular permeability κ.This results in Robin boundary conditions at the top and the bottom of the channel, of the form where n is the unit outward normal vector at the vessel walls.
As predicted by Richardson et al. [13] via matched asymptotic expansions, the parameter κ plays an important role in determining the position of the particles deposited onto the vessel wall.In fact, their outer solution shows how small values of the permeability are responsible for the formation of a boundary layer region in the immediate vicinity of the wall where the advective flux balances the diffusive flux and the thickness of the vessel wall prevents particles from flowing out and this may hamper the proposed treatment.In this work, a reference value of κ = O(10 −6 ) is chosen [40].
Assuming that particles enter the channel due to an injection of 3 s duration in the vicinity of x = 0, we have where the concentration c in (y, t) describes an injection in the central region of the vessel (see Appendix 6).At the channel outlet However, this outlet boundary condition is irrelevant since particles never reach the outlet of the vessel in our simulations.The initial particle concentration in the channel is zero, hence c(x, y, 0) = 0.

Simulation of the nanoparticle concentration
To solve ( 29) numerically, we first nondimensionalise the equation (see Appendix 7).Then, we use an explicit Euler scheme in time, with first order upwind approximations for the advection terms and central differences for the second order derivatives (see Appendix 8 for more details on the numerical scheme).The velocity u F is presented in Section 2.3, where the Newtonian and Ellis fluid models yield analytical expressions for u F while the Carreau model velocity is determined numerically.The vertical particle velocity v p is calculated via equation ( 26), with a constant value of the magnetic force F 0 .Advection terms, as usual, will dominate where the flow rate is nonnegligible but the diffusive contribution becomes important near the wall of the vessel, where the fluid has nearly zero velocity (see discussion on boundary layers in Appendix 7).For the simulations shown in the next section, we use a square unit domain with a resolution of 400 nodes in the x and y directions.The time step, ∆t, is suitably chosen to satisfy the stability requirements of the scheme.
At the beginning of the process, there are no particles in the vessel, that is c(x, y, 0) = 0, but at t = 0 we inject a concentration of particles c in (y, t) at the central region of the boundary x = 0. To avoid the numerical difficulty associated with an abrupt change in concentration at this boundary, we assume that the concentration c in (y, t) progressively increases with time until it reaches a maximum value, c 0 , at t = 3 s (see Appendix 6).

Mesh convergence
Here we provide results of a spatial and temporal convergence study before presenting numerical predictions based on the numerical scheme given in Appendix 8. Following the notation from Appendix 8, all variables used in this subsection are dimensionless.The reference numerical approximation, c ref , is generated on a fine mesh with n x = 400 and n y = 400 grid points in the x and y directions, respectively, with n x ∆x = n y ∆y = 1, and time step ∆t 4 where ∆t k = ∆t 0 /2 k−1 , k = 1, . . ., 4, and ∆t 0 = 0.2×10 −5 .Let M n denote the mesh with (25×2 n ) 2 cells with n = 1, . . ., 4.
To investigate convergence with respect to mesh size, the l 2 -norm of the difference between the numerical approximations on meshes M 4 and M n is computed for n = 1, . . ., 3 at the end of the simulation t = 0.8.This is denoted by E n .Since the concentration at the lower boundary y = −1/2 is important in our study, we also compute the l 2 -norm at this boundary denoted by E lb n .These measures of the error are tabulated in Table 5.They confirm second-order convergence of the numerical approximation with respect to mesh spacing as shown in the final column where the rate of convergence p n is computed using The influence of mesh size on the evolution of the average nanoparticle flux through the bottom boundary is shown in Figure 4 in dimensional units for the Ellis model with magnetic force F 0 = 0.5 × 10 −14 N. Also shown in Figure 4 is an application of Richardson extrapolation based on extrapolating the approximations computed on meshes M 2 , M 3 and M 4 .This demonstrates convergence of this quantity with mesh refinement with enhanced accuracy obtained through use of Richardson extrapolation.
Next we investigate convergence with respect to time step.In this temporal convergence study we use the finest mesh M 4 and compute the l 2 -norm of the difference between the numerical approximations obtained using time steps ∆t 4 and ∆t k for k = 1, . . ., 3 at time t = 0.8.These measures of the error are tabulated in Table 6.They confirm first-order convergence of the numerical approximation with respect to time step as shown in the final column where the rate of convergence p k is computed using

Results and discussion
We first analyse the effects that the different fluid models have on the motion of the particles in the bloodstream.As demonstrated in Section 2.1 the power-law model gives unrealistic values for the viscosity and therefore the results for this case will not be discussed.The results presented here, correspond to numerical simulations of equation ( 29) applying an injection for three seconds at x = 0.

Nanoparticle migration in Newtonian blood flows
Figure 5 shows the evolution of the concentration of nanoparticles in the vessel under the influence of a magnetic force F 0 = 0.5 × 10 −14 N when the Newtonian approximation for blood is considered.The panel on top represents the velocity field, u tot .Since the velocity of the particles depends on the viscosity of the fluid and the Newtonian flow assumes a constant value, µ then v p = F 0 /6πaµ, the changes in the velocity field are only due to the parabolic profile u F (y) and only in the y-direction.Furthermore, we observe how the particles near the vessel wall experience a much smaller drag force with respect to those at the center of the vessel (i.e., the horizontal component of the velocity is smaller near the wall than at the center) and, therefore, will react strongly to the magnetic force and deviate from the typical parabolic behaviour of the fluid.The colour maps show snapshots of the concentration of particles at five different times.We observe that the particles entering the bloodstream from the vessel inlet are immediately driven to the lower wall of the vessel, where some of them will permeate the vessel wall.At t = 6.4 s and t = 8 s, we observe the formation of a boundary layer where particles accumulate in the immediate proximity of the wall (e.g., at t = 8 s large values of c appear near x = 0.2 • 10 −3 m, at the lower boundary of the vessel), as theoretically predicted in [13].

Nanoparticle migration in non-Newtonian blood flows
In Figures 6 and 7 the analogous plots to those of Figure 5 are shown for the Carreau and the Ellis models, respectively.As expected, both models yield very similar particle concentrations.However, in contrast to the Newtonian case, not all the nanoparticles are forced out at the lower boundary.This change in behaviour may be attributed to the high viscosity values characterising the non-Newtonian models in the vicinity of y = 0 (see Figure 3b), due to which the particles find it difficult to escape this 'central' region.Since the vertical velocity is proportional to 1/η, in the Newtonian model the vertical velocity is significant and much larger than the corresponding Figure 5: Snapshots of the concentration of magnetic nanoparticles, c/c 0 , in a capillary vessel at five different times (t = 1.6s, t = 3.2s, t = 4.8s, t = 6.4s, t = 8s), using the Newtonian model for blood flow and with a constant magnetic force equal to F 0 = 0.5 × 10 −14 N. The top panel represents the velocity field (equation ( 17)).Other parameter values as in Table 3. one for the Ellis and Carreau models where the higher viscosity values lead to significantly lower vertical velocities in the 'central' region.Near the boundaries the viscosity of the Ellis and Carreau models takes a similar value to the Newtonian fluid and, hence, so does the vertical velocity.This may be observed from the velocity profiles shown in the top panels.The differences found in the concentration maps from Fig. 5-7 are mainly caused by the changes in the viscosity between models and clearly demonstrate that a Newtonian fluid model severely overestimates the effect of the magnetic force on the nanoparticles.
Note the advection term in (21) is the dominant mechanism for nanoparticle transport in all fluid models studied (see discussion in Appendix 7), but diffusion also contributes to the observed reduction in concentration, as particles spread out.This can be observed in all three cases (Figures 5-7) where the colored central region increases size with time.Actually, diffusion allows some reduced number of nanoparticles reach the upper vessel wall (see, for instance, panels for t = 6.4 s in Figures 5-6).

Varying the strength of the magnetic field
Having demonstrated that the non-Newtonian models are more realistic, we now aim to find a value of the magnetic force F 0 capable of driving the same number of particles to the lower wall of the vessel as the Newtonian fluid, for the same capillary vessel length and radius as in previous simulations.For doing so, we analyse the average particle flux crossing the lower vessel wall as a function of time, which is computed using Since the Ellis and Carreau models provide similar results, for simplicity we use the Ellis model as a representative non-Newtonian fluid.In Figure 8 we present the evolution of the average particle flux (33) with time for the Newtonian (solid line) and Ellis (dashed line) fluids taking F 0 = 0.5 × 10 −14 N. We also simulate a third case corresponding to the Ellis fluid (dash-dotted line) with a magnetic force F 0 = 1 × 10 −14 N, i.e., a magnetic force two times larger.As expected, the average flux for the Ellis fluid with F 0 = 0.5 × 10 −14 N is much lower than the Newtonian fluid for the same magnetic force.For instance, at t = 8 s the average flux for the Ellis fluid is a 41% lower than the Newtonian one.However, by doubling the strength of the magnetic field, we observe that the flux for the Ellis model gets closer to the Newtonian case with F 0 = 0.5 × 10 −14 N.
To facilitate the comparison between these two cases, we numerically integrate the corresponding curves over the total simulation time, obtaining 2.4977 × 10 −7 mol/m 2 and 2.4992 × 10 −7 mol/m 2 for the Newtonian with F 0 = 0.5 × 10 −14 N and the Ellis with F 0 = 1 × 10 −14 N fluid, respectively.Hence, increasing the magnetic force to F 0 = 1 × 10 −14 N is enough to attract the same number of particles to the lower vessel wall in the chosen simulation time.
In Figure 9, we present the particle concentration maps for the Ellis model with F 0 = 1×10 −14 N. By comparing them with those in Figure 7, we observe that the number of particles pushed to the lower boundary is clearly larger, consistent with the higher average flux of particles shown in Figure 8.Since increasing F 0 increases the vertical velocity, the formation of the boundary layer in the vicinity of the lower wall, where the y-advection term is dominant, is even clearer than in Figure 7.

Sensitivity under changes of the vessel wall permeability
We note that different tissues may have different permeabilities, and recent experiments indicate effective permeability values of the order of 10 −6 m/s or smaller in tumors [40].As shown in Figure 10, a decrease in the permeability in our model results in a consistent decrease of the particle flux.Also shown in Figure 10 is a sensitivity analysis with respect to the physically meaningful value of the permeability, i.e., κ = 1 • 10 −6 m/s.Increasing or decreasing this reference value of the permeability by 10% with respect the reference value κ = 1 • 10 −6 m/s generates a change in the particle flux of around 5% at the end of the simulation (t = 8 s).This demonstrates that the results are robust with respect to small changes to the value of vascular permeability.Figure 6: Snapshots of the concentration of magnetic nanoparticles, c/c 0 , in a capillary vessel at five different times (t = 1.6s, t = 3.2s, t = 4.8s, t = 6.4s, t = 8s), using the Carreau model for blood flow and with a constant magnetic force equal to F 0 = 0.5 × 10 −14 N. The top panel represents the velocity field (see eq. ( 19)).Other parameter values as in Table 3.
Figure 7: Snapshots of the concentration of magnetic nanoparticles, c/c 0 , in a capillary vessel at five different times (t = 1.6s, t = 3.2s, t = 4.8s, t = 6.4s, t = 8s), using the Ellis model for blood flow and with a constant magnetic force equal to F 0 = 0.5 × 10 −14 N. The top panel represents the velocity field (equation ( 20)).Other parameter values as in Table 3.
Figure 8: The average particle flux at the lower wall for the Newtonian and Ellis models versus time with magnetic force F 0 = 0.5 × 10 −14 N (solid and dashed lines, respectively), and the Ellis model for a doubled force, F 0 = 1 × 10 −14 N (dash-dotted line).

Conclusions
We have formulated a model that describes the motion of magnetic nanoparticles in a blood vessel subject to an external magnetic field in order to optimize magnetic drug targeting.The model consists of a system of nonlinear partial differential equations formed by the Navier-Stokes equations for the flow of blood and with an advection-diffusion equation for the concentration of nanoparticles.We consider Newtonian flow and three different non-Newtonian flows.Assuming a long two-dimensional vessel, the equations are significantly reduced and the system is then solved via analytical and numerical techniques.
We have accounted for all the forces involved in this physical process, combined with realistic choices for the parameter values.It has been shown that, in order to correctly simulate the delicate balance between hydrodynamic (Stokes drag) and magnetic forces in the vessel, it is crucial to choose an appropriate non-Newtonian model for blood.
The first part of this paper examines the non-Newtonian behaviour of blood and the importance of choosing an appropriate model for the fluid viscosity.The Newtonian approximation proved to be inaccurate while the more commonly used power-law model exhibits an unbounded value for the viscosity at the centre of the vessel.The Carreau and Ellis models are both found to be good models for simulating blood behaviour, and they lead to very similar predictions for the velocity of the fluid.However, the Ellis model has a distinct advantage in that it permits an analytical expression for the fluid velocity, so we focus on it for this reason.
In the second part magnetic nanoparticles were introduced in the flow.With such small particles both advection and diffusion effects can play an important role.The key result of the paper is that a Newtonian model predicts greater particle motion than an Ellis or Carreau model under the same external magnetic force.Specifically, for one case examined, a magnetic force had to be applied in the Ellis fluid that was twice as large than in the Newtonian fluid for the same number of particles to be absorbed through the vessel wall.These results show that models based on a Newtonian fluid can drastically overestimate the effect of the magnetic field and therefore we conclude that in order to accurately model nanoparticle drug targeting in realistic clinical situations the non-Newtonian behaviour of blood needs to be accounted for so that a sufficiently strong magnetic field is applied.
In future work the model can be improved by imposing pulsatile flow, including the elasticity of the vessels due to the change in pressure and also by solving in more complex geometries.Furthermore, it can be combined with detailed experimental studies to optimize the delivery of Figure 9: Snapshots of the concentration of magnetic nanoparticles, c/c 0 , in a capillary vessel at five different times (t = 1.6s, t = 3.2s, t = 4.8s, t = 6.4s, t = 8s), using the Ellis model for blood flow and with a constant magnetic force equal to F 0 = 1 × 10 −14 N. The top panel represents the velocity field (equation ( 20)).Other parameter values as in Table 3.In our simulations M = 20.

Appendix B: Nondimensionalisation of the diffusion equation
Using the non-dimensional variables (10) and introducing new ones into (29) we obtain Rearranging the terms and choosing δt = L/U to balance the time derivative with the advection term we obtain where δ = W/U and Pe = 2RU/D is the Péclet number.Pe depends on the model chosen.The time scale δt represents the average time taken by a particle to pass through the vessel in the absence of a magnetic field (e.g., in the Newtonian case δt = 6.67 s).According to our choice of scales and the values in Table 4, O((ε Pe) −1 ) ≈ 10 −1 depending on the fluid chosen, while ε = O(10 −2 ).Hence, both terms on the right-hand side of (38) are small.Therefore, the advective terms are dominant and when analysing them, it is important to understand the order of magnitude of the fraction δ/ε.In particular we can distinguish three regions in the domain which does not affect the results since the simulations are always stopped when the particles are still far from reaching the vessel outlet.This ensures that nanoparticles only leave the vessel through absorption at the boundaries y = ±1/2 during our simulations.
8 Appendix C: Finite difference scheme for the diffusion equation We drop the hats in ( 38)-(41) and define c n i,j := c(x i , y j , t n ), u Fj := u F (y j ), v pj := v p (y j ), D Tj := D(y j ). (42) The choice of the direction of the upwind step is made considering that the solution of the velocity of the fluid is always non negative and the direction of the vertical velocity is always negative (since we have positioned the magnet below the vein).Then, equation ( 38) can be approximated as We also need to specify the boundary conditions at x = 0, 1 and y = ±1/2.On the wall of the vessel, that is the upper and the lower side of the rectangle, we approximated conditions (39), through the three-point backward difference formula

Figure 1 :
Figure 1: Sketch of the injection of magnetic nanoparticles in a vessel, in the presence of red blood cells and subject to a magnetic field.

Figure 4 :
Figure 4: Influence of mesh size on the evolution of the average nanoparticle flux through the boundary y = −R for the Ellis model with magnetic force F 0 = 0.5 × 10 −14 N. Also shown is an application of Richardson extrapolation.

Figure 10 :
Figure 10: Evolution of the average nanoparticle flux through the boundary y = −R as a function of time for F 0 = 1 × 10 −14 N and three different permeabilities using the Ellis model.The shaded area represents the sensitivity of the nanoparticle flux for variations of the permeability within a 10% around κ = 1 • 10 −6 m/s.

Figure 11 :
Figure 11: Profiles of c in at four different times (from bottom to top: t = 0.75 s, t = 1.5 s, t = 2.25 s and t = 3 s).In our simulations M = 20.
(symmetric with respect to the center of the vessel): a central region where O(δ) < O(ε), which is the broadest one where the drag force dominates over the magnetic force; a second region where O(δ) ≈ O(ε) where both advective terms balance; finally, very near to the wall of the vessel, we can find a narrow boundary layer where O(δ) > O(ε) (Stokes drag is very small) and vertical motion due to the magnetic field dominates.The boundary conditions at the vessel wall are

T j± 1 2
are evaluated by the arithmetic mean.