Magnetic nanodrug delivery in non-Newtonian blood flows
MMagnetic nanodrug delivery in non-Newtonian blood flows
C. Fanelli , K. Kaouri , T.N. Phillips , T.G. Myers , , F. Font , February 10, 2021
Abstract
With the goal of determining strategies to maximise drug delivery to a specific site in the body,a mathematical model for the transport of drug nanocarriers (nanoparticles) in the bloodstreamunder the influence of an external magnetic field is presented. Under the assumption of long(compared to the radius) blood vessels the Navier-Stokes equations are reduced, consistently withlubrication theory. Under these assumptions, analytical results are compared for Newtonian, power-law, Carreau and Ellis fluids, which clearly demonstrate the importance of shear thinning effectswhen modelling blood flow. Incorporating nanoparticles and a magnetic field to the model wedevelop a numerical scheme and study the particle motion for different field strengths. In particular,we demonstrate the importance of the non-Newtonian behaviour: for certain field strengths aNewtonian model predicts that all particles are absorbed into the vessel wall, whilst the non-Newtonian models predict that a small number of particles are absorbed into the wall while therest flow along with the blood and leave the vessel at the outlet. Consequently, models based ona Newtonian fluid can drastically overestimate the effect of a magnetic field. Finally, we evaluatethe particle concentration at the vessel wall and compute the evolution of the particle flux throughthe wall for different permeability values, which is important when assessing the efficacy of drugdelivery applications. The insights from our work bring us a step closer to successfully transferringmagnetic nanoparticle drug delivery to the clinic.
Currently, the main approaches in cancer therapy include surgery, chemotherapy, radiotherapy, andhormone therapy [1]. The first is invasive while the latter three are non-specific. Hence their efficacyis not only reduced but sometimes the cure can be more aggressive than the disease itself. In thelate 1970s, magnetic drug targeting, which consists of injecting and steering magnetic drug carriersthrough the vessel system towards disease locations, was proposed as an alternative and more efficienttreatment 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 for magneticdrug delivery on human patients have been carried out [6, 7, 8]. All these studies demonstrate thatmagnetic forces can attract particles in the region near the magnet but there is a lack of knowledgeon 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 relativelydifficult task due to the interplay of magnetic and hydrodynamic forces acting on the particles and theinherent difficulty of solving the Navier-Stokes equations. One of the main difficulties is to accuratelymodel and simulate the non-Newtonian behaviour of blood. Experiments show that for low shear Departament de Matem`atica Aplicada, Universitat Polit`ecnica de Catalunya, Barcelona, Spain. School of Mathematics, Cardiff University, CF24 4AG, UK Centre de Recerca Matem`atica, Campus de Bellaterra Edifici C, 08193 Bellaterra, Barcelona, Spain. ff[email protected] a r X i v : . [ phy s i c s . f l u - dyn ] F e b ates ( ∼ − ), the viscosity can be as high as 10-11 mPa · s while for shear rates in excess of 1000 s − ,it tends to an asymptotic value of 3-4 mPa · s [10, 11]. In order to simplify calculations, previousstudies 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 continuummodel for the motion of particles subject to a magnetic field, both using a Newtonian flow model forblood. In [12] it was shown, via a simple network model, that it is impossible to specifically targetinterior regions of the body with an external magnetic field; the magnet can be used only for targetsclose 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 amagnetic field using a power-law assumption for the blood flow. Cherry and Eaton [17] developed acomprehensive continuum model for the motion of micro-scale particles using a fitted value for theviscosity from experimental data [20] to model blood shear thinning, and used the model to investigatemagnetic particle steering through a branching vessel.A structural difference within the models developed in the literature concerns the way in whichparticles are described. Instead of describing particle distribution in blood as a continuum and using adiffusion equation to study the evolution of the particle concentration, several studies consider particlesas 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 whichplay a crucial role in magnetic drug targeting, showing how difficult it can be to keep small particlesin the desired region. Recently, in Boghi et al. [15] a numerical simulation of drug delivery in theblood in the coeliac trunk was performed.In the present work we analyse the forces and parameters involved in the process of magnetic drugdelivery and highlight the importance of considering realistic non-Newtonian models for blood flowand the motion of the nanoparticles in it. It is crucial to consider the non-Newtonian behaviour ofthe blood in order to predict if particles are able to reach the desired area. A mathematical model ina two-dimensional channel is introduced in Section 2. In Section 2.1, it is explained how choosing anoversimplified constitutive law for the fluid in the vessel or a wrong value for the viscosity in the centreof the channel can lead to significant errors. In Section 3 we show how magnetic forces act on particlesdepending on their size. Numerical simulations illustrating the influence of some key parameters arepresented in Section 4. As illustrated in Figure 1, we approximate the vessel as a long and thin rectangular channel, consistentwith the sizes of arteries and vessels in the human body (see Table 1). The particles are injected atthe inlet of the vessel and their motion is driven by the field generated by a magnet located at thebottom of the domain as well as the drag force on the particles due to the blood flow.Blood containing nanoparticles can be considered as a nanofluid. Typically, its motion is governedby the Navier-Stokes equations coupled to an advection-diffusion equation for the concentration ofparticles in the fluid [24]. The Navier-Stokes equations describing the flow of an incompressiblenanofluid in a vessel subject to an external magnetic field are ρ (cid:20) ∂ u F ∂t + u F · ∇ u F (cid:21) = −∇ p + ∇ · τ + F mag , (1) ∇ · u F = 0 , (2)2 igure 1: Sketch of the injection of magnetic nanoparticles in a vessel, in the presence of red blood cells andsubject to a magnetic field. Vessel L (m) R (m) N Aorta 0 . . × − . . × − × − . × − . × Capillary 6 × − × − . × Venule 8 × − × − . × Vein 0 . . × − .
22 1 . × − Table 1: Typical values for the various types of vessels in human body: average length ( L ) and width ( R ),estimation of the average number of a specific vessel in the circulatory system ( N ). All values are adapted from[23]. where u F is the fluid velocity, ρ is the fluid density, p is pressure, τ is the stress tensor and F mag themagnetic force acting on the fluid.In this work, we consider that the nanofluid, composed of blood and nanoparticles, is a dilutemixture. 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. Secondly, giventhat 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 ≈ Bodily fluids, such as blood, saliva, eye fluid are invariably non-Newtonian. In particular, blood is aconcentrated suspension of particles in plasma, which is mainly made of water. The three most impor-tant particles that constitute blood are: red blood cells (or erythrocytes), white cells (or leukocytes)and platelets (or thrombocytes). The first ones, which are the most numerous, are mainly responsi-3le for the mechanical properties of the blood [23]. In fact, their tendency to form (and then breakdown) three-dimensional microstructures at low shear rates and to align to the flow at high shear ratescause the blood’s shear thinning behaviour, characterized by the monotonic decrease of the viscositythat tends to some limit for very high shear rates [11]. In the case of blood, the structures lead tosignificant changes in its rheological properties and several models have been developed during thepast 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 beexpressed as linear function of the components of the velocity gradient. The more general relationbetween the stress and rate-of-strain tensors can be written as τ = η ( ˙ γ ) ˙ γ , (3)where η ( ˙ γ ) is the viscosity, ˙ γ = ∇ u + ∇ u T and˙ γ = (cid:32) ∂u i ∂x j + ∂u j ∂x i (cid:33) (cid:32) ∂u i ∂x j + ∂u j ∂x i (cid:33) / , (4)is the generalised shear rate. When η ( ˙ γ ) is constant the Newtonian model is retrieved. A purely shear-thinning fluid will exhibit a monotonic decrease in viscosity with increasing shear rate. Practical fluidsare more likely to exhibit a constant viscosity beyond certain high and low values of shear rate. Inbetween a nonlinear viscosity relation links the two ‘Newtonian plateaus’.The aim of this section is to compare different types of simple, non-Newtonian fluid models. Inparticular, we will compare the power-law model (used in earlier studies of blood flow with magneticparticles), the Carreau model and the Ellis model with the Newtonian model. In Table 2 and Table 3we summarize the expressions for the viscosity and shear stress for each model and the correspond-ing parameter values, respectively. Note that under the assumption of flow in a long thin channel,significant reduction in the governing equations is possible. The shear stress relations in Table 2 areconsistent with this reduction. In the power-law model m is constant. If n p < n p > n p = 1 we retrieve theNewtonian expression. In the Carreau model λ is a constant and η and η ∞ are the limiting viscositiesat low and high shear rates, respectively. In the Ellis model η is the viscosity at zero shear and τ / isthe shear stress at which the viscosity is η /
2. This model does not include the high viscosity plateau,however for most practical situations such high strain rates are never reached and so this plateau valueis irrelevant. 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 with the parameter values of Table 3.For moderate shear rates it is clear that the Ellis and Carreau models show good agreement withdifferences beginning as the shear rate increases above ˙ γ ≈ . s − . In attempting to compensate forthe plateau values the power-law model appears inaccurate over the whole range of shear rates.In the next section we investigate the velocity profiles using the above models and demonstrate theimportance of choosing the right model for blood flowing under the influence of an external magneticfield. 4 odel Viscosity Shear stress Newtonian µ = constant τ yx = µ (cid:12)(cid:12)(cid:12) ∂u∂y (cid:12)(cid:12)(cid:12) Power-law η p ( ˙ γ ) = m | ˙ γ | n p − τ yx = m (cid:12)(cid:12)(cid:12) ∂u∂y (cid:12)(cid:12)(cid:12) n p Carreau η c ( ˙ γ ) = η ∞ + ( η − η ∞ ) (cid:0) λ ˙ γ (cid:1) nc − τ yx = η ∞ + ( η − η ∞ ) (cid:18) λ (cid:12)(cid:12)(cid:12) ∂u∂y (cid:12)(cid:12)(cid:12) (cid:19) nc − (cid:12)(cid:12)(cid:12) ∂u∂y (cid:12)(cid:12)(cid:12) Ellis η e ( τ ) = η (cid:18) (cid:12)(cid:12)(cid:12) ττ / (cid:12)(cid:12)(cid:12) α − (cid:19) − η e ( τ yx ) = η (cid:18) (cid:12)(cid:12)(cid:12) τ yx β (cid:12)(cid:12)(cid:12) α − (cid:19) − Table 2: Newtonian and non-Newtonian viscosity models and corresponding shear stress, assuming ˙ γ ≈ | ∂u/∂y | . Quantity Symbol Value Units References
Newtonian viscosity µ . m n p λ η η ∞ n c η .
056 Pa s [35]Ellis shear stress at η / β α Table 3: Typical parameter values for the blood flow equations for each model, taken from [34, 35, 11].
Expressing equations (1)–(2) into their components, the steady-state governing equations for the fluidbecome ∂u∂x + ∂v∂y = 0 , (5) ∂u∂t + u ∂u∂x + v ∂u∂y = − ρ ∂p∂x + 1 ρ (cid:20) ∂τ xx ∂x + ∂τ yx ∂y (cid:21) , (6) ∂v∂t + u ∂v∂x + v ∂v∂y = − ρ ∂p∂y + 1 ρ (cid:20) ∂τ xy ∂x + ∂τ yy ∂y (cid:21) . (7)which are subject to the no-slip condition and to symmetry conditions u = 0 , v = 0 at y = ± R . (8) ∂u∂y = 0 at y = 0 , (9)5 Figure 2: The viscosity/shear rate plot in logarithmic scale for the power-law (red dotted-dashed line), theCarreau (black solid line) and the Ellis models (green dashed line). The light grey dashed lines represents theconstant limit viscosities η and η ∞ . respectively.In order to understand the order of magnitude of the terms in equations (5)–(7), we proceed tonon-dimensionalise them. We define the non-dimensional variables x = L ˆ x, y = R ˆ y, u = U ˆ u, v = V ˆ v, p = P ˆ p, τ = T ˆ τ . (10)From (5) it is clear that we can balance the equation by setting V = εU , where ε = R/L (cid:28) ∂ ˆ u∂ ˆ x + ∂ ˆ v∂ ˆ y = 0 . (11)Choosing P = L T /R , we can write εU (cid:20) ˆ u ∂ ˆ u∂ ˆ x + ˆ v ∂ ˆ u∂ ˆ y (cid:21) = T ρ (cid:20) − ∂ ˆ p∂ ˆ x + ε ∂ ˆ τ ˆ x ˆ x ∂ ˆ x + ∂ ˆ τ ˆ y ˆ x ∂ ˆ y (cid:21) , (12) ε U (cid:20) ˆ u ∂ ˆ v∂ ˆ x + ˆ v ∂ ˆ v∂ ˆ y (cid:21) = T ρ (cid:20) − ∂ ˆ p∂ ˆ y + ε ∂ ˆ τ ˆ x ˆ y ∂ ˆ x + ε ∂ ˆ τ ˆ y ˆ y ∂ ˆ y (cid:21) . (13)Assuming εU ρ T (cid:28) O ( ε ) or smaller in (5)–(7), the resulting governingequations, in dimensional form, are ∂u∂x + ∂v∂y = 0 , (14) ∂p∂x = ∂τ yx ∂y , (15) ∂p∂y = 0 , (16)which is a much simpler system of equations. 6 .3 Comparison of flow for different viscosity models To determine the most accurate model to describe blood flow, we now compare the fluid velocityprofiles obtained for the four different viscosity laws summarised in Table 2. In all cases we considera Hagen-Poiseuille flow in a channel, which is driven by a constant pressure gradient ∆ p/L and whichcorresponds to a constant volumetric flow rate Q . The results quoted below follow from the modelsdescribed in [35] for flows driven by a pressure gradient and a moving bottom surface. They maybe obtained by setting the velocity of the surface to zero. To switch between the current verticalco-ordinate and that used in [35] we set z = ( y + R ) h/ (2 R ) and the position of the turning point inthe flow, y = 0, then corresponds to z = z m = h/ The Newtonian model:
The solution for the Newtonian fluid is obtained by solving equations(14)–(16) where the shear stress specified in Table 2. Since equation (16) implies that the pressuredoes not vary with y , integrating (15) twice with respect to y and using the boundary conditions(8)-(9) we obtain u ( y ) = 3 Q R (cid:16) R − y (cid:17) , (17)where the relation between ∂p/∂x and Q is established through Q = (cid:82) R − R u dy = constant (see [35]). The power-law model:
Proceeding equivalently to the Newtonian case with the correspondingshear stress from Table 2, the velocity for a non-Newtonian power-law fluid becomes u ( y ) = QR n +1 n (cid:18) n + 12 n + 2 (cid:19) (cid:16) R n +1 n −| y | n +1 n (cid:17) . (18)The main disadvantage of this model is that since ∂u/∂y = 0 at y = 0 the viscosity η ∝ / ( ∂u/∂y ) →∞ . The Carreau model:
Choosing the Carreau model for the viscosity, the equations of the flow arereduced to (14) and (16) coupled with ∂p∂x = ∂∂y η ∞ + ( η − η ∞ ) (cid:32) λ (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂y (cid:12)(cid:12)(cid:12)(cid:12) (cid:33) n − (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂y (cid:12)(cid:12)(cid:12)(cid:12) . (19)This expression cannot be integrated analytically, here 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 expressedanalytically as u ( y ) = 1 η ∂p∂x (cid:34) R −| y | (cid:18) β ∂p∂x (cid:19) α − R α +1 −| y | α +1 α + 1 (cid:35) . (20)The corresponding flux is Q = 1 η ∂p∂x (cid:34) R (cid:18) β ∂p∂x (cid:19) α − R α +2 α + 2 (cid:35) . (21)In order to compare the four models and their velocity field we take typical parameter valuesfrom the literature. The viscosity depends on temperature, so for the whole paper we use parametersconsistent with a typical body temperature, 37 ° C, indicated in Table 3. The geometry parameters are7 a) (b)
Figure 3: Comparison of (a) velocity and (b) viscosity profiles of Newtonian (represented by the letter N andthe blue dotted line), power-law (represented by PL and the red dashed-dotted line), Carreau (represented bythe letter C and the black solid line) and Ellis model (represented by the letter E and the green dashed line). those of an artery with length and width as in Table 1. The flow entering in a small artery can beapproximated by Q ≈ . × − m / s [11].The velocity and corresponding viscosity profiles for blood obtained from the four different modelsare compared in Figures 3(a) and 3(b), respectively. We observe that although the difference invelocity profiles is small for all models, the corresponding viscosities can differ significantly. Giventhat the magnetic particle motion is dependent on the viscosity of the fluid it travels through it isclear that both Newtonian and power-law models are inappropriate. The Ellis and Carreau modelsshow excellent agreement in both velocity and viscosity profiles. In the literature the Carreau model isgenerally preferred due to its ability to predict both Newtonian plateaus. However for in vivo bloodflow it is unlikely that the high shear plateau will be reached suggesting that the Ellis model providesan adequate approximation. This is an important result, since the Ellis model permits an analyticalexpression for the flow. This not only provides a better understanding of the factors affecting the flowbut also considerably simplifies the numerical study. In § The behaviour of the concentration of magnetic nanoparticles in the bloodstream is obtained followingthe continuum model developed by Grief and Richardson [12]. The governing equation describing themotion of magnetic particles in the blood stream is an advection-diffusion equation for the particleconcentration c ( x, y, t ): ∂c∂t + ∇ · (cid:2) ( u F + u p ) c (cid:3) = ∇ · J diff , (22)where u F is the fluid velocity (as discussed in the previous section), u p is the particle velocity and J diff the diffusion flux.The diffusive character of the equation is due to Brownian motion and shear-induced diffusion.Shear-induced diffusion arises due to the fact that the red blood cells suspended in plasma collidewith each other causing random motion with a diffusive character. As recently demonstrated by Liu8nd coworkers [36, 37] via a lattice-Boltzmann based multiscale simulation, the diffusion due to theBrownian motion in the case of nanoparticles transport in a small vessel can be important and, in somecases, even the predominant diffusion process. In particular, for nanoparticles smaller than 100 nmthey found that D Br /D sh (cid:29)
1, where D Br and D sh correspond to the Brownian and shear-induceddiffusion coefficients, respectively, indicating that Brownian diffusion may be dominant for particles inthis size range [37]. Using the Stokes-Einstein equation for the diffusion of spherical particles througha shear thinning fluid [11], we can write the Brownian diffusion coefficient as D Br = k B T6 π η ( ˙ γ ) a , (23)where k B is Boltzmann constant, T the absolute temperature and a the particle radius. On the otherhand, the shear-induced diffusion contribution can be approximated by D sh = K sh ( r RBC ) ˙ γ, (24)where K sh is a dimensionless coefficient that depends on the blood cell concentrations and r RBC isthe red blood cell radius. The coefficient K sh is difficult to measure but the value used in Table 4 isconsidered representative in the literature [12]. Hence, the total diffusion will be D tot = D Br + D sh ,leading to J diff = − D tot ∇ c = − (cid:18) k B T6 π η ( ˙ γ ) a + K sh ( r RBC ) ˙ γ (cid:19) ∇ c. (25)The particle velocity is found by balancing hydrodynamic and magnetic forces. Using the definitionof the Stokes drag, representing the hydrodynamic force on a spherical particle of radius a movingthrough a viscous fluid, we have F St = 6 πa η ( ˙ γ ) u p . (26)The particles reach equilibrium velocity when F St balances the magnetic force F mag and this leads usto the expression u p = F mag πa η ( ˙ γ ) , (27)where F mag = 2 πa µ χ χ/ ∇| H | . (28)We will deal with a constant magnetic force acting perpendicular to the flow, i.e. F mag = F j .Therefore, in combination with (27) we obtain that u p = (0 , v p ( y )) where v p = F / (6 πa η ( ˙ γ )). Finally,since the vertical fluid velocity is negligible, u tot = ( u F ( y ) , v p ( y )) and D tot = D tot ( y ) the concentrationequation (22) becomes ∂c∂t + u F ∂c∂x + ∂ ( v p c ) ∂y = D tot ∂ c∂x + ∂∂y (cid:18) D tot ∂c∂y (cid:19) . (29)We assume that the particles can flow out of the vessel through the walls with a certain vascularpermeability κ . This results in Robin boundary conditions at the top and the bottom of the channel,of the form ( u tot c − D tot ∇ c ) · n = κ c on y = ± R , (30)where n is the unit normal vector at the vessel walls.As predicted by Richardson et al. [13] via matched asymptotic expansions, the parameter κ playsan important role in the determination of the position of the particles deposited onto the vessel wall. Infact, their outer solution shows how small values of the permeability are responsible for the formation9f a boundary layer region in the immediate vicinity of the wall where the advective flux balances thediffusive flux and the thickness of the vessel wall prevents particles from flowing out. In this work, areference value for κ = O (10 − ) is chosen in order to represent a type of wall where particles are ableto pass through but a small boundary layer is observed.Assuming that the flux entering in the channel is constant and equal to the inlet concentration c in for a given interval of time (cid:2) , T inj (cid:3) , we set c ( x, y, t ) | x =0 = (cid:40) c in if 0 ≤ t ≤ T inj , (cid:18) u F c − D tot ∂c∂x (cid:19) (cid:12)(cid:12)(cid:12) x = L = 0 , (32)since the end of the channel is sufficiently far from the magnet that its effect is negligible. The initialparticle concentration in the channel is zero, hence c ( x, y,
0) = 0.
Quantity Symbol Value Units References
Blood cell radius r RBC . × − m [12]Coeff. blood cell conc. K sh × − No. [12]Inlet concentration c in − [16]Particle radius a × − m [12]Boltzmann constant k B . × − m kg s − K − [11]Temperature T 310.15 K [11] Table 4: Characteristic parameter values for the advection-diffusion equation.
To solve (29) numerically we use an explicit Euler scheme in time, with first order upwind approxima-tions for the advection terms and central differences for the second order derivatives (see AppendixB for more details on the numerical scheme). The velocity u F is defined by the results of § v p is calculated via equation (27) witha constant value of the magnetic force F . Advection terms, as usual, will dominate where the flowrate is non-negligible but the diffusive contribution becomes important near the wall of the vessel,where the fluid has nearly zero velocity (see boundary layer discussion in Appendix A.2). For thesimulations shown in the next section, we used a grid of length L = 0 . R = 3 · − mcomposed by n x × n y = 100 ×
100 nodes. The spatial steps are ∆ x = L/ ( n x −
1) = 0 . · − m and∆ y = 2 R/ ( n y −
1) = 0 . · − m, and ∆ t is suitably chosen to satisfy the stability requirements ofthe scheme.At the beginning of the process, there are no particles flowing in the vessel, that is c ( x, y,
0) = 0, butat t = 0 we inject a concentration of particles c ( x , y , t ) = c in . To overcome the numerical difficultyof this abrupt change in concentration, a small time analysis of (29) is performed (see Appendix A.2).Then, instead of the initial condition c ( x, y,
0) = 0 to initialise the numerical scheme we use c (¯ x, y, ¯ t ) = 14 π ¯ t ¯ D tot exp (cid:34) − (¯ x − ¯ u ¯ t ) D tot ¯ t − ( y − δ ¯ v p ¯ t ) D tot ¯ t (cid:35) , (33)10here ¯ D tot , ¯ u F and ¯ v p are the corresponding averaged values of D tot , u F and v p , and ¯ t = t/ε ¯ x = x/ε for a given ε (cid:28) We first analyse the effects that the different fluid models have on the motion of the particles in thebloodstream. As demonstrated in section 2.1 the power-law model gives unrealistic values for theviscosity 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 at x = 0 for five seconds,that is T inj = 5 s.Figure 4 shows the evolution of the concentration of nanoparticles in the vessel when the Newtonianapproximation for blood is considered. The panel on top represents the velocity field, which is given bythe contribution of the drag of the fluid and the magnetic force acting on the system. Since the velocityof the particles depends on the viscosity of the fluid and the Newtonian flow assumes a constant value,the changes in the velocity field are only due to the parabolic profile of u F and, therefore, only in the y -direction. The colour maps show snapshots of the concentration of particles under the influence of amagnetic force F = 1 × − N at six different times. We observe that all the particles entering thebloodstream from the vessel inlet are driven to the lower wall of the vessel and then permeate the wallhalfway down the channel after about 30 s. Moreover, the diffusive contribution plays an importantrole in this region and influence the motion of the particles at this scale. Especially in the snapshotsfor t = 6 s and t = 15 s, we observe a thin boundary layer where particles accumulate in the immediateproximity of the wall (e.g., at t = 15 s large values of c appear near x = 0 .
02 at the lower boundary ofthe vessel), as theoretically predicted in [13].In Figures 5 and 6 the equivalent plots to those of Figure 4 are shown but now using the Carreauand the Ellis models, respectively. As expected both models show very similar behaviour for theconcentration of particles. However, in contrast to the Newtonian case, the nanoparticles are notforced out at the lower boundary. The reduction in concentration observed is a result of diffusion,as the particles spread out. This change in behaviour may be attributed to the high viscosity valuespredicted by the non-Newtonian models in the vicinity of y = 0, due to which the particles find itdifficult to escape the ‘central’ region. This may be observed from the velocity profiles shown inthe top panels. For the Newtonian flow the vertical velocity (proportional to 1 /η ) is significant andfor Ellis and Carreau the higher viscosity values lead to a significantly lower vertical velocity. Thisclearly demonstrates that a Newtonian fluid model overestimates the effect of the magnetic force onthe nanoparticles.Having demonstrated that the non-Newtonian models are more realistic, we now aim to find avalue of the magnetic force F capable of driving all the particles to the lower wall of the vessel, forthe same artery length and radius than the previous simulations and in a time-window <
60 s. InFigures 7 and 8 we show how increasing the magnetic force to F = 8 × − N is enough to attractall particles to the vessel wall in the specified time window: after approximately 40s the concentrationis negligible. Furthermore, looking at the velocity profiles in Figures 7-8 (top panels) we observe howthe particles near the vessel wall experience a much smaller drag force with respect to those at thecenter of the vessel (i.e., the horizontal component of the velocity is smaller near the wall than at thecenter) and, therefore, will react strongly to the magnetic force and deviate from the typical parabolicbehaviour of the fluid.An important aspect for drug delivery applications is the number of particles that will cross theartery wall, which will depend on the evolution of the concentration of the particles at the wall (fora given magnetic force) as well as on the wall permeability. In Figure 9a, we show the evolutionof the concentration at the wall, c | y = − R , using the Ellis model for the case F = 8 × − N and κ = 1 × − m/s. The concentration profiles show that the number of particles at the boundary11 = 0 s t = 4 s t = 6 s t = 10 s t = 15 s t = 30 s Figure 4: Snapshots of the concentration of particles in an artery at six different times, using the Newtonianmodel for blood flow and with a constant magnetic force equal to F = 1 × − N. The top panel representsthe velocity field. = 0 s t = 4 s t = 10 s t = 15 s t = 30 s t = 50 s Figure 5: Concentration of magnetic nanoparticles in an artery at six different times, using the Carreau modelfor blood flow and with a constant magnetic force acting on particles equal to F = 1 × − N. The top panelrepresents the velocity field. = 0 s t = 6 s t = 10 s t = 15 s t = 30 s t = 50 s Figure 6: Snapshots of the concentration of particles in an artery at six different times, using the Ellis modelfor blood flow and with a constant magnetic force equal to F = 1 × − N. The top panel represents thevelocity field. = 0 s t = 2 s t = 6 s t = 15 s t = 30 s t = 40 s Figure 7: Concentration of nanoparticles in an artery at six different times, using the Carreau model for bloodflow and with a constant magnetic force acting on particles equal to F = 8 × − N. The top panel representsthe velocity field. = 0 s t = 2 s t = 6 s t = 15 s t = 30 s t = 40 s Figure 8: Concentration of magnetic nanoparticles in an artery at six different times using the Ellis model forblood flow and with a constant magnetic force F = 8 × − N. The top panel represents the velocity field. igure 9: Left panel: Concentration profiles at the boundary for the case F = 8 × − N and κ = 1 × − m/sat four different times using the Ellis flow model. Right panel: Evolution of the average nanoparticle flux throughthe boundary y = − R as a function of time for F = 8 × − N and four different permeabilities using theEllis model. decreases with time and confirms that after about 40 s the particle concentration at the boundary isnegligible (for t = 44 .
56 s the concentration peak is < .
01 mol/m ). In Figure 9b, we show the averageparticle flux crossing the artery wall as a function of time for several permeability values, which iscomputed using c | y = − R via κ (cid:104) c (cid:105) = κL (cid:90) L c | y = − R dx . (34)Taking as a reference the case κ = 1 × − m/s, we observe that the flux crossing the boundaryinitially increases until t = 5 s due to the injection of particles (note the small kink indicates the endof particle injection). The subsequent peak, at around 35 s corresponds to the time when the attenuatedconcentration peak observed in the snapshots from Figures 7-8 (see transition from t = 30 s to t = 40 s)crosses the boundary. After this, the particle flux starts to decrease rapidly.We note that different tissues may display different permeabilities, and recent experiments indicateeffective permeability values of the order of 10 − m/s or smaller in tumors [38]. As shown in Figure 9b,a decrease in the permeability in our model results in a consistent decrease of the particle flux. We have formulated a model for the motion of magnetic nanoparticles in a vessel subject to an externalmagnetic field in order to optimize the technique of magnetic drug targeting. The model consists ofa system of nonlinear partial differential equations formed by the Navier-Stokes equations for theflow of blood coupled with an advection-diffusion equation for the concentration of nanoparticles. Weconsider Newtonian flow and three different non-Newtonian flows. Assuming a long two-dimensionalvessel, the equations are simplified and the system is solved via analytical and numerical techniques.The aim of this paper has been to account for all the forces acting in the physical process, combinedwith realistic choices for parameters. It has been shown that, in order to correctly simulate the delicatebalance between hydrodynamic and magnetic forces in the vessel, it is crucial to choose an appropriatenon-Newtonian model for blood.The first part of this paper deals with the non-Newtonian behaviour of blood and the importance ofchoosing the right model for the fluid viscosity. The Newtonian approximation proved to be inaccurate17nd the more commonly used power-law model exhibits an unbounded value for the viscosity at thecenter of the vessel. The Carreau and Ellis models are found to be the best approach to simulatingblood behaviour, both providing very similar results. However, the Ellis model has a distinct advantagein that it permits an analytical expression for the fluid velocity.In the second part magnetic particles were introduced into the flow. With such small particlesBrownian motion and shear-induced diffusion effects can play an important role. The key result ofthe paper is that a Newtonian model predicts much greater particle motion than an Ellis or Carreaumodel. Specifically, for one case examined, under the influence of the same external magnetic field, theNewtonian model predicted that all particles would pass through the vessel wall while both Carreauand Ellis models predicted that particles would simply flow out through the vessel outlet. These resultsshow that models based on a Newtonian fluid can drastically overestimate the effect of the magneticfield and therefore we conclude that in order to accurately model nanoparticle drug targeting inrealistic clinical situations the non-Newtonian behaviour of blood needs to be accounted for.In future work this model can be improved by imposing pulsatile flow, including the elasticity ofthe vessels due to the change in pressure and solving in more complex geometries. Furthermore, it canbe combined with detailed experimental studies to optimize the delivery of drugs to specific regionsor other applications.
Acknowledgements
CF acknowledges financial support from the Spanish Ministry of Economy and Competitiveness,through the “Maria de Maeztu” Programme for Units of Excellence in R&D (MDM-2014-0445). TMacknowledges the support of the Spanish Ministry of Science and Innovation Grant No. MTM2017-82317-P. FF acknowledges the financial support of the “Juan de la Cierva” programme (grant No.IJC2018-038463-I) from the Spanish Ministry of Science and Innovation. FF and TM acknowledgethe support of the CERCA Programme of the Generalitat de Catalunya.
A Nondimensionalisation of the diffusion equation and small timelimit
A.1 Nondimensionalisation
Using the non-dimensional variables (10) and introducing new ones t = T ˆ t, c = c ˆ c, u F = U ˆ u F , v p = W ˆ v p , D tot = D ˆ D tot , (35)into (29) we obtain c T ∂ ˆ c∂ ˆ t + U c L ˆ u F ∂ ˆ c∂ ˆ x + W c R ∂ (ˆ v p ˆ c ) ∂ ˆ y = Dc L (cid:32) ˆ D tot ∂ ˆ c∂ ˆ x (cid:33) + Dc R ∂∂ ˆ y (cid:18) ˆ D tot ∂ ˆ c∂ ˆ y (cid:19) . (36)Rearranging the terms and choosing T = L/U to balance the time derivative with the advection termto obtain ∂ ˆ c∂ ˆ t + ˆ u F ∂ ˆ c∂ ˆ x + δε ∂ (ˆ v p ˆ c ) ∂ ˆ y = 1 ε Pe ε (cid:32) ˆ D tot ∂ ˆ c∂ ˆ x (cid:33) + ∂∂ ˆ y (cid:18) ˆ D tot ∂ ˆ c∂ ˆ y (cid:19) , (37)where δ = W/U and Pe = D/ ( RU ) is the P´eclet number which depends on the model chosen.According to the values in Table 4, O ( ε Pe) − ≈ [10 − , − ] depending on the fluid chosen, while ε = O (10 − ), hence both terms on the right-hand side of the equation are small. Therefore, the18dvective terms are dominant and when analysing them, it is important to understand the order ofmagnitude of the fraction δ/ε . In particular we can distinguish three regions in the domain (symmetricwith respect to the center of the vessel): a central region where O ( δ ) < O ( ε ), which is the broadestone 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 narrowboundary layer where O ( δ ) > O ( ε ) and vertical motion dominates. A.2 Small time analysis
Let us consider a constant diffusion coefficient ¯ D tot , where the bar represents the average value of thesum of the Brownian and shear-induced diffusive contributions divided by Pe. Since the interest is tounderstand the solution for small times and close to the inlet, let t = ε ¯ t and x = ε b ¯ x , where ε (cid:28) b is to be determined. Introducing these new variables in equation (37), we have1 ε ∂c∂ ¯ t + 1 ε b ¯ u F ∂c∂ ¯ x + δε ¯ v p ∂c∂y = ¯ D tot ε (cid:32) ε ε b ∂c ∂ ¯ x + ∂c ∂y (cid:33) , (38)which can be rewritten as ∂c∂ ¯ t + ε − b ¯ u F ∂c∂ ¯ x + δ ¯ v p ∂c∂y = ¯ D tot (cid:32) ε − b ∂c ∂ ¯ x + ∂c ∂y (cid:33) . (39)In order to remove the dependence on ε , we impose b = 1, i.e. x = ε ¯ x . With the new scale, the leadingorder terms of equation (39) will be ∂c∂ ¯ t + ¯ u F ∂c∂ ¯ x + δ ¯ v p ∂c∂y = ¯ D tot (cid:32) ∂c ∂ ¯ x + ∂c ∂y (cid:33) . (40)In an unbounded domain, the solution of (40) is c (¯ x, y, ¯ t ) = 14 π ¯ t ¯ D tot exp (cid:34) − (¯ x − ¯ u ¯ t ) D tot ¯ t − ( y − δ ¯ v p ¯ t ) D tot ¯ t (cid:35) , (41)and since for small times the particles are far from the vessel walls, we can use this information forthe initial value of the system and avoid numerical problems due to the initial jump of the continuum. B Finite difference scheme for the diffusion equation
We will use the notation c ni,j := c ( x i , y j , t n ) , u F j := u F ( y j ) , v p j := v p ( y j ) , D T j := D tot ( y j ) . (42)The choice of the direction of the upwind step is made considering that the solution of the velocityof the fluid is always non negative and the direction of the velocity of the particle is always negative(since we have positioned the magnet below the vein). Then, equation (29) can be approximated as c n +1 i,j − c ni,j ∆ t + u F j (cid:32) c ni,j − c ni − ,j ∆ x (cid:33) + (cid:32) v p j +1 c ni,j +1 − v p j c ni,j ∆ y (cid:33) =+ D T j (cid:32) c ni +1 ,j − c ni,j + c ni − ,j ∆ x (cid:33) + D T j + 12 (cid:32) c ni,j +1 − c ni,j ∆ y (cid:33) − D T j − (cid:32) c ni,j − c ni,j − ∆ y (cid:33) , (43)19here D T j ± are evaluated by the arithmetic mean. We also need to specify the inflow boundaryconditions at x = 0 and y = ± R . On the vessel inlet, where the concentration of drugs is injected fora limited interval of time, we have that c n ,j = (cid:40) c in for 0 ≤ t n ≤ T inj , j = 1 , . . . , n y , since we stop injecting nanoparticles at T inj .On the wall of the vessel, that is the upper and the lower side of the rectangle, we approximatedconditions (30), by which particles can pass through the vein with a constant permeability coefficientequal to κ , through the three-point backward difference formula v p ny c n +1 i,n y − D T ny c n +1 i,n y − c n +1 i,n y − + c n +1 i,n y − y = − κc n +1 i,n y , (45)and the three-point forward difference formula v p c n +1 i, − D T (cid:32) − c n +1 i, + 4 c n +1 i, − c n +1 i, y (cid:33) = κc n +1 i, , (46)for i = 1 , . . . , n x . Finally, the vessel outlet is approximated again through the three-point forwarddifference formula u F j c n +1 n x ,j − D T j (cid:32) − c n +1 n x − ,j + 4 c n +1 n x − ,j − c n +1 n x ,j x (cid:33) = 0 , (47)for j = 1 , . . . , n y . References [1] B. Bahrami, M. Hojjat-Farsangi, H. Mohammadi, E. Anvari, G. Ghalamfarsa, M. Yousefi, andF. Jadidi-Niaragh. Nanoparticles and targeted drug delivery in cancer therapy.
Immunologyletters , 190:64–83, 2017.[2] Q.A. Pankhurst, J. Connolly, S.K. Jones, and J. Dobson. Applications of magnetic nanoparticlesin biomedicine.
Journal of Physics D: Applied Physics , 36(13):R167–R181, 2003.[3] C. Alexiou, R. Jurgons, R. Schmid, A. Hilpert, C. Bergemann, F. Parak, and H. Iro. In vitroand in vivo investigations of targeted chemotherapy with magnetic nanoparticles.
Journal ofMagnetism and Magnetic Materials , 293:389–393, 05 2005.[4] M. Muthana, S. Scott, N. Farrow, F. Wright, C. Murdoch, S. Grubb, N. Brown, J. Dobson, andC. Lewis. A novel magnetic approach to enhance the efficacy of cell-based gene therapies.
Genetherapy , 15:902–10, 07 2008.[5] G. Wang, S. Gao, R. Tian, J. Miller-Kleinhenz, Z. Qin, T. Liu, L. Li, F. Zhang, Q. Ma, andL. Zhu. Theranostic hyaluronic acid-iron micellar nanoparticle for magnetic field-enhanced invivo cancer chemotherapy.
ChemMedChem , 13, 10 2017.[6] A. L¨ubbe, C. Bergemann, H. Rieß, F. Schriever, P. Reichardt, K. Possinger, M. Matthias,B. D¨orken, F. Herrmann, R. G¨urtler, P. Hohenberger, N. Haas, R. Sohr, B. Sander, J. Lemke,D. Ohlendorf, W. Huhnt, and D. Huhn. Clinical experiences with magnetic drug targeting: Aphase I study with 4 ´ -epidoxorubicin in 14 patients with advanced solid tumors. Cancer research ,56:4686–93, 11 1996. 207] M. Wilson, R. Kerlan, N. Fidelman, A. Venook, J. Laberge, J. Koda, and R. Gordon. Hepato-cellular carcinoma: Regional therapy with a magnetic targeted carrier bound to doxorubicin in adual MR imaging/ conventional angiography suite-initial experience with four patients.
Radiology ,230:287–93, 02 2004.[8] A.J. Lemke, M.I. Pilsach, A. L¨ubbe, C. Bergemann, H. Rieß, and R. Felix. MRI after magneticdrug targeting in patients with advanced solid malignant tumors.
European radiology , 14:1949–55,12 2004.[9] S. Fiocchi, E. Chiaramello, M. Bonato, G. Tognola, D. Catalucci, M. Parazzini, and P. Ravazzani.Computational simulation of electromagnetic fields on human targets for magnetic targeting ap-plications. In , pages 5674–5677, 2019.[10] Hideki Yamamoto, Takafumi Yabuta, Yusuke Negi, Daiki Horikawa, Kimito Kawamura, EijiTamura, Katsuhiro Tanaka, and Fujimaro Ishida. Measurement of human blood viscosity a usingfalling needle rheometer and the correlation to the modified herschel-bulkley model equation.
Heliyon , 6(9):e04792, 2020.[11] R. Fournier.
Basic Transport Phenomena in Biomedical Engineering, 4th edition . CRC Press, 112017.[12] A.D. Grief and G. Richardson. Mathematical modelling of magnetically targeted drug delivery.
Journal of Magnetism and Magnetic Materials , 293(1):455–463, 2005.[13] G. Richardson, K. Kaouri, and H.M. Byrne. Particle trapping by an external body force in thelimit of large P´eclet number: applications to magnetic targeting in the blood flow.
EuropeanJournal of Applied Mathematics , 21(1):77–107, 2010.[14] P. Yue, S. Lee, S. Afkhami, and Y. Renardy. On the motion of superparamagnetic particles inmagnetic drug targeting.
Acta Mechanica , 223, 03 2012.[15] A. Boghi, F. Russo, and F. Gori. Numerical simulation of magnetic nano drug targeting in apatient-specific coeliac trunk.
Journal of Magnetism and Magnetic Materials , 437:86 – 97, 2017.[16] A. Nacev, C. Beni, O. Bruno, and B. Shapiro. The behaviors of ferromagnetic nano-particlesin and around blood vessels under applied magnetic fields.
Journal of magnetism and magneticmaterials , 323(6):651–668, 2011.[17] E. K. Cherry and J. Eaton. A comprehensive model of magnetic particle motion during magneticdrug targeting.
International Journal of Multiphase Flow , 59:173–185, 02 2014.[18] T. Lunoo and T. Puangmali. Capture efficiency of biocompatible magnetic nanoparticles inarterial flow: A computer simulation for magnetic drug targeting.
Nanoscale Research Letters ,10:426, 10 2015.[19] I. Rukshin, J. Mohrenweiser, P. Yue, and S. Afkhami. Modeling superparamagnetic particles inblood flow for applications in magnetic drug targeting.
Fluids , 2(2):29, 2017.[20] D.E. Brooks, J.W. Goodwin, and G.V. Seaman. Interactions among erythrocytes under shear.
Journal of Applied Physiology , 28(2):172–177, 1970.[21] J. B. Freund and B. Shapiro. Transport of particles by magnetic forces and cellular blood flowin a model microvessel.
Physics of Fluids , 24(5):051904, 2012.2122] H. Ye, Z. Shen, and Y. Li. Computational modeling of magnetic particle margination withinblood flow through lammps.
Computational Mechanics , 62(3):457–476, 2018.[23] L. Formaggia, A. Quarteroni, and A. Veneziani.
Cardiovascular Mathematics: Modeling andSimulation of the Circulatory System , volume 1. Springer-Verlag (Milano), 01 2009.[24] T.G. Myers, H. Ribera, and V. Cregan. Does mathematics contribute to the nanofluid debate?
International Journal in Heat and Mass Transfer , 111:279–288, 2017.[25] M.M. MacDevette, T.G. Myers, and B. Wetton. Boundary layer analysis and heat transfer of ananofluid.
Microfluidics and nanofluidics , 17(2):401–412, 2014.[26] H. Alimohamadi and M. Imani. Transient non-Newtonian blood flow under magnetic targetingdrug delivery in an aneurysm blood vessel with porous walls.
International Journal for Compu-tational Methods in Engineering Science and Mechanics , 15(6):522–533, 2014.[27] P.A. Voltairas, D. Fotiadis, and L. Michalis. Hydrodynamics of magnetic drug targeting.
Journalof biomechanics , 35:813–21, 07 2002.[28] P.D. Ballyk, D.A. Steinman, and C. Ethier. Simulation of non-Newtonian blood flow in anend-to-side anastomosis.
Biorheology , 31:565–86, 09 1994.[29] S. Chien. Shear dependence of effective cell volume as a determinant of blood viscosity.
Science(New York, N.Y.) , 168:977–9, 06 1970.[30] Y. Cho and K. Kensey. Effects of the non-Newtonian viscosity of blood on flows in a diseasedarterial vessel. part 1: Steady flows.
Biorheology , 28:241–62, 02 1991.[31] C.R. Huang and W Fabisiak. Thixotropic parameters of whole human blood.
Thrombosis research ,8:1–8, 06 1976.[32] E. Kaliviotis, J. Sherwood, and S. Balabani. Local viscosity distribution in bifurcating microfluidicblood flows.
Physics of Fluids , 30:030706, 03 2018.[33] J. Sherwood, D. Holmes, E. Kaliviotis, and S. Balabani. Spatial distributions of red blood cellssignificantly alter local haemodynamics.
PloS one , 9:e100473, 06 2014.[34] B.M. Johnston, P.R. Johnston, S. Corney, and D. Kilpatrick. Non-Newtonian blood flow in humanright coronary arteries: steady state simulations.
Journal of Biomechanics , 37(5):709 – 720, 2004.[35] T.G. Myers. Application of non-Newtonian models to thin film flow.
Physical Review E , 72:066302,2005.[36] Z. Liu, J.R. Clausen, R.R. Rao, and C.K. Aidun. Nanoparticle diffusion in sheared cellularblood flow.
Journal of Fluid Mechanics , 871:636–667, 2019.[37] Z. Liu, Y. Zhu, R.R. Rao, J.R. Clausen, and C.K. Aidun. Nanoparticle transport in cellular bloodflow.
Computers and Fluids , 172:609 – 620, 2018.[38] Marc Lim, Vishnu Dharmaraj, Boying Gong, Benson T. Jung, and Ting Xu. Estimating tumorvascular permeability of nanoparticles using an accessible diffusive flux model.