Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2022.017257

ARTICLE

Performance Analysis of Magnetic Nanoparticles during Targeted Drug Delivery: Application of OHAM

Muhammad Zafar1,#,*, Muhammad Saif Ullah1,#, Tareq Manzoor2, Muddassir Ali3, Kashif Nazar4, Shaukat Iqbal5, Habib Ullah Manzoor6, Rizwan Haider1 and Woo Young Kim7,*

1Institute of Energy and Environmental Engineering, University of the Punjab, Quaid-e-Azam Campus, Lahore, 54590, Pakistan
2Energy Research Centre, COMSATS University Islamabad, Lahore, 54000, Pakistan
3Department of Energy Engineering, Faculty of Mechanical and Aeronautical Engineering, University of Engineering and Technology, Taxila, 47080, Pakistan
4Department of Mathematics, COMSATS University Islamabad, Lahore, 54000, Pakistan
5School of Systems and Technology, University of Management Sciences and Technology, Lahore, 54770, Pakistan
6Department of Electrical, Electronics and Telecommunication Engineering, University of Engineering and Technology, Lahore (FSD Campus), Faisalabad, 38000, Pakistan
7Department of Electronic Engineering, Faculty of Applied Energy System, Jeju National University, Jeju Special Self-Governing Province, Jeju-si, 63243, Korea
*Corresponding Authors: Muhammad Zafar. Email: zaffarsher@gmail.com; Woo Young Kim. Email: semigumi@jejunu.ac.kr
#These authors contributed equally to this work
Received: 26 April 2021; Accepted: 03 August 2021

1  Introduction

Nanotechnology is effectively playing its significant role in numerous fields such as structure, environment, molecular physics, chemistry, biology, material sciences, computer sciences, engineering, measurements, imaging, and many other disciplines of science and technology [1–6]. Likewise, the application of nanoparticles for drug delivery is also one of the critical forefronts of medical sciences. The use of nanoparticles or nanotubes, as detectors or biosensors, can be exploited as they suffer changes in their respective electrical properties upon application. Moreover, under the influence of the magnetic field, such particles can be called as magnetic nanoparticles.

Magnetic nanoparticles typically consist of nickel, cobalt, or iron and are clusters of magnetic particles where several individual particles are present. By the application of synthesized Fe3C, magnetic hyperthermia can be studied. Similarly, the intrinsic loss power value and a specific absorption rate can also be determined. It is a well-known fact that these aspects are much superior at lower magnetic nanoparticle concentrations [7].

Since blood can be considered a bio-magnetic fluid, prone to be affected upon the magnetic field’s application, the exploitation of magnetic fluid properties can help improve drug targeting through the transportation of magnetic particles in the blood [8]. The application of magnetic nanoparticles in biomedicine offer advantages of stability and size controllability of magnetic nanoparticles. The magnetic nanoparticle applications have been investigated in several diseases such as cardiovascular, endovascular, drug targeting, cancer tissues, and hyperthermia treatment for malignant tumor cells [9–11]. It has been investigated that the blood’s viscosity with magnetic particles caries from the blood under the magnetic field’s influence [12]. The magnetization of fluid under the magnetic field’s control, without the electric field’s induction, affects bio-magnetic fluid flow [12–15]. However, most of these studies did not account for the instability of blood flow [16]. Various mathematical models considered that the blood flow followed non-Newtonian fluid properties. The study of non-Newtonian blood flow characteristics in micro-tube under the magnetic field’s influence is carried out by Shaw et al. [17]. However, Bandyopadhyay et al. [18] studied the Newtonian model for investigating the blood flow through a constricted channel. More interestingly, Shit et al. [19] studied the pulsatile blood flow in an artery environment under the influence of magnetic dipole. Considering Newtonian fluid properties, various studies are conducted to investigate the transportation of magnetic nanoparticles, capture efficiency in a tube for targeted drug delivery, and the effect of nanoparticle augmentation [20–25]. The majority of these studies are focused on pressure gradient and electrokinetic force, which necessitated the further exploration of magnetic nanoparticles and blow, together, in a channel, under the influence of the magnetic field, electrokinetic force, and pressure gradient.

Recently, a semi-analytic approximate method emerged for handling time-dependent partial differential equations, known as the optimal homotopy asymptotic method (OHAM). Compared to numerical methods, the optimal homotopy asymptotic method (OHAM) is powerful and straightforward. It achieves outcomes more quickly while maintaining acceptable results and good agreement with numerical methods. This method combines the advantages of the homotopy principle with an efficacious computational algorithm that gives a thorough and straightforward approach for controlling solution convergence [26–28]. The main advantage of OHAM is that it is free of parameters, does not require identifying the ℏ-curve and any guess initial like in other methods. Moreover, the OHAM consists of inherent convergence criteria like the homotopy analysis method (HAM), with greater flexibility [29]. The OHAM is not only applicable for small parameters, but it also demonstrates its reliability and high accuracy to solve nonlinear systems in engineering and science. There are also few disadvantages of the OHAM. One of the disadvantages is that a set of nonlinear algebraic equations must be solved in each order of approximations. Another limitation is that OHAM consists of several unidentified convergence-control variables, making computations time-consuming [30]. Several nonlinear coupled partial and ordinary differential systems can be accurately solved using this approach. Many recent studies have considered using this technique to solve a variety of nonlinear problems [31–33]. However, to the best of the authors’ knowledge, no attempt is made to investigate magnetic particles and blood movement through a micro-tube under the influence of external forces using the OHAM.

The present study provides insight into the movement of magnetic nanoparticles in a tube with blood flow under the collective effect of the electrokinetic force, magnetic field, and pulsatile pressure gradient. It is well-established that the ability to predict the capture of nanoparticles under the influence of a magnetic field in vivo is essential. In this study, we solved the governing equations to predict blood flow and magnetic particle velocity using the optimal homotopy asymptotic method approach. The parametric study is used here to compute the solutions for the non-dimensional velocity of magnetic particles and blood. These studies’ findings may further explore magnetic therapy viability with magnetic nanoparticles to help design medical devices and drug delivery systems for blood clot removal, infection treatment, and tumor cell treatment.

2  Methodology

2.1 Mathematical Formulation

In the current study, a two-dimensional mathematical formulation based on an unsteady and incompressible blood flow together with magnetic nanoparticles in a cylindrical vessel under the presence of a uniform external magnetic field and an axial electric field is considered. The schematic illustration of the physical problem is presented in Fig. 1, where the r∗ is the radial coordinate perpendicular to flow direction. The z∗ axis is taken along the vessel’s axis. The length and radius of the vessel are taken as L and h, respectively. It is assumed that the flow is axisymmetric and developed fully. A coherent magnetic field is employed perpendicular to the direction of the flow. It is also assumed that the wall zeta potential is constant in an electrolyte solution of blood near the vessel wall. The pressure gradient (∂p∗/∂z∗) is applied along the axis of the vessel and a no-slip boundary on the wall of the vessel. It is assumed that the micro-vessel length is significantly larger than its radius, and the magnetic nanoparticles are uniformly dispersed in the bloodstream. For this case, the magnetic Reynolds number (Rem=vreflref/vm, where vref is the reference velocity, lref is the reference length scale and vm is the magnetic diffusivity) is very low (Rem<<1), and hence the effect of the induced magnetic field can be neglected during the present study.

Figure 1: Schematic representation of the physical problem (Reproduced with permission from A. Mondal and G. C. Shit, Journal of Magnetism and Magnetic Materials; published by Elsevier B. V., 2017)

2.1.1 Distribution of Electrical Potential

For at any point (r∗,z∗) in the vessel, the total electric potential (ϕ∗) is written as [16]

ϕ∗(r∗,z∗)=ψ∗(r∗)−z∗Ez∗ (1)

where ψ∗(r∗) and Ez∗ are the potential distribution due to the presence of an electrical double layer (EDL) and applied electric field, respectively. Electrical potential distribution within the cylindrical coordinate is defined using the famous Poisson equation form

1r∗∂∂r∗(r∗∂ϕ∗∂r∗)+∂2ϕ∗∂z∗2=−ρeϵ (2)

where ρe and ϵ is the net charge density and permittivity of the medium, respectively.

The overall charge density distribution can be given as the sum of all cations and anions in the solution, assuming the equilibrium Boltzmann distribution,

ρe=∑izien0exp⁡(−zieψ∗kBTav) (3)

where zi, e, no, kB, and Tav are the valence of type i ions, a charge of an electron, bulk ionic concentration, Boltzmann constant, and absolute temperature, respectively. Assume that the thermal energy of the ions is much larger than the electric potential, i.e., zveψ∗kBTav≪1, for which the Debye-Hückel linearization principle can be estimated as sinh(zveψ∗kBTav)≈zveψ∗kBTav. According to this assumption, the overall charge density is the electric charge density can be further expressed as

ρe=−2n0zv2e2ψ∗kBTav (4)

Accordingly, Eq. (2) can be rewritten as

1r∗ddr∗(r∗dψ∗dr∗)=ψ∗λD2 (5)

where 1/λD is known as Debye length (EDL thickness) and is described as 1/λD=(2n0zv2e2ϵkBTav)0.5. The corresponding boundary conditions are applied to ψ∗(r∗) are

ψ∗=ξatr∗=h (6)

Now, introducing a normalized electro-osmotic potential function ψ consisting of zeta potential ξ of the medium and the non-dimensional coordinates in the form

ψ=ψ∗ξ,r=r∗h,z=z∗h (7)

where h is the characteristic radius of the micro-vessel.

In terms of the non-dimensional variables described in Eqs. (5), (7) gives

d2ψdr2+1rdψdr−k2ψ=0 (8)

where k=h/λD is known as the electro-osmotic parameter. The non-dimensional form of the corresponding boundary conditions (6) is expressed as

ψ=1atr=1 (9)

The electric potential distribution solution from the Eq. (8) by Eq. (9) can be given in the form of the modified Bessel function of the first kind as

ψ(r)=I0(kr)I0(k) (10)

where Io is known as the zeroth-order modified Bessel function of the first kind.

2.1.2 Flow Investigation

The flow of blood together with the magnetic nanoparticle velocity inside a micro-vessel (assuming cylindrical polar coordinates) is expressed by the governing equation of electron magnetohydrodynamic [22,34,35]

∂u∗∂t∗=−1ρ∂p∗∂z∗+v(∂2u∗∂r∗2+1r∂u∗∂r∗)+ksNρ(v∗−u∗)−σB02u∗ρ+ρeEZ∗ρ (11)

We assumed that a spherical magnetic nanoparticle motion in a viscous carrier fluid is governed by the applied magnetic and electric field effect. The viscous drag and magnetic forces during slow motion are responsible forces, and thus ignoring the Brownian motion, the movement of the nanoparticles is given by the second law of motion (Newton’s) [35–38].

m∂y∂x=Fm+Ff+Fe (12)

where m(=4/3πrp3ρp with a particle radius rp and particle density ρp) is a single nanoparticle mass, v∗(r∗,t∗) is a single magnetic nanoparticle velocity, u∗(r∗,t∗) is the axial velocity of blood, Bo is the strength of the applied magnetic field, p∗ is the fluid pressure, ρ is the blood density, v is the kinematic viscosity coefficient, σ is the electrical conductivity, N is the magnetic nanoparticles number per unit volume and Ks=6πμEhyd,pVr (with an effective hydrodynamic radius of the particle Ehyd,p and magnetic drift velocity Vr) is Stokes constant [36]. The terms on the right side of the Eq. (11) are the forces resulting from the relative motion between magnetic nanoparticles and fluid, the Lorentz force due to the magnetic field and the net electrical body force. The terms Fm, Ff, and Fe on the right side of Eq. (12) are the magnetic, fluidic, and electric forces, respectively.

If the mass of a particle is likely to be negligible, the inertia term m∂v∗∂t∗ can be neglected, and therefore the motion of the particle must fulfill the force balance Fm+Ff+Fe=0.

By using Stoke’s law for the drag on a sphere, the fluidic force encountered by a particle (considering low Reynolds number) is estimated by

Ff=−6πμEhyd,pVr(v∗−u∗) (13)

If the force acting on the particle is fluidic in one dimension [22,36], i.e., Ff expressed by Eq. (13), the Eq. (12) becomes

m∂v∗∂t∗=Kr(u∗−v∗) (14)

Along the axial direction the pulsatile pressure gradient can be given as [18]

−∂p∗∂z∗=A0+A1cos(ω′t∗) (15)

where ω′ is the pulsation frequency, A0 is the pressure gradient constant amplitude, and A1 is the pulsatile component amplitude that increases periodic pressure in a micro-vessel.

The boundary conditions for magnetic nanoparticles and fluid velocity can be put mathematically in the form of

u∗=v∗=0att∗=0, (16)

The following non-dimensional variables can be introduced as

u=u∗UHS,v=v∗UHS,t=t∗vh2,p=p∗hUHSμ (17)

with UHS=−ϵEz∗ξ/μ is the classical electro-osmotic flow velocity, also known as the Helmholtz-Smoluchowski velocity [39]. In this case, non-dimensional variables described in Eq. (17), the Eqs. (11) and (14) take the following form:

∂u∂t=(A0+A1cosα2t)+(∂2u∂r2+1r∂u∂r)+R(v−u)−Hα2u+kψ (18)

G∂v∂t=(u−v) (19)

and non-dimensional boundary conditions take the following form:

u=v=0att=0, (20)

where Ha=B0hσ/μ is known as the Hartmann number, α2=ω′h2/v is the Womersley number, R=ksNh2/μ is the particle concentration parameter, and G=mμ/ρh2ks is defined as the particle mass parameter.

The dimensionless mathematical expression for the wall shear stress (WSS) τw and volumetric flow rate Q is given as, respectively

τw=(−∂u∂r)r=1 (21)

Q=2π∫01u(r,t)rdr (22)

2.1.3 Numerical Method

The optimal homotopy analysis method for the equation is as follows [40]:

with boundary conditions

2.1.4 Principle of OHAM

This section describes the OHAM principle, and the following differential equation can be used to demonstrate the basic concept of the OHAM.

L(u(τ))+N(u(τ))+g(τ)=0 (23)

It depends on the following boundary condition: B(u)=0.

where g, B, N, and L are the known analytical function, boundary, nonlinear, and linear operators, respectively.

First of all, build a series of equations using the OHAM [41,42]

(1−p)[L(ϕ(τ,p))+g(τ)]−H(p)×[L(ϕ(τ,p))+g(τ)+N(ϕ(τ,p))]=0 (24)

with the boundary condition

B(ϕ(τ,p))=0 (25)

where p∈[0,1], ϕ(τ,p), and H(p) are the embedded parameter, unknown function, and nonzero auxiliary function for p≠0 and H(0)=0, respectively. Thus, at p=0 and p=1, it follows that

ϕ(τ,0)=u0(τ),ϕ(τ,1)=u(τ) (26)

Accordingly, as p rises from 0 to 1, the initial solution u0(τ) moves closer to solution u(τ). Now,

The auxiliary function can be considered now, which can be expressed in the following way:

H(p)=pC1+p2C2+p3C3… (27)

where C1, C2, C3 are the constants, which will be later decided.

In the following way, expanding ϕ(τ,p) in a series in terms of p

ϕ(τ,p,Ci)=u0(τ)+∑k≥1uk(τ,Ci)pk,i=1,2,3… (28)

A series of differential equations with boundary conditions are derived by putting Eq. (28) into Eq. (24), gathering the same power of p and equating the coefficient of like power p to zero. As a result, by solving a differential equation with boundary conditions, we get u0(τ), u1(τ,C1), u2(τ,C3), … etc. In general, the solution of Eq. (23) can be approximated in the following way:

u~(m)=u0(τ)+∑k≥1uk(τ,Ci) (29)

The following expression for residual is obtained by substituting Eq. (29) in Eq. (23)

R(τ,Ci)=L(u~(m)(τ,Ci))+g(τ)+N(u~(m)(τ,Ci))=0 (30)

If R(τ,Ci)=0 then u~(τ,Ci) represents the precise solution to the problem. In most cases, especially in nonlinear problems, this does not occur, but we can reduce the functional

J(C1,C2,C3,…Cn)=∫abR2(τ,C1,C2,C3,…Cm)dτ (31)

The constants Ci are obtained using the following conditions:

∂J∂C1=∂J∂C2=…=0 (32)

These constants are used to determine the approximate solution (of order m) from Eq. (29).

2.1.5 Application of OHAM

This section describes the OHAM application to the differential Eq. (3). Moreover, we can build homotopy of the Eq. (3) according to the OHAM [39] as follows:

(1−p)d2ψdr2−(C1p+C2p2+C3p3){d2ψdr2+1rdψdr−k2ψ}=0 (33)

Considering ψ as follows:

ψ=ψ0(r)+pψ1(r)+p2ψ2(r) (34)

then putting the value of ψ from the Eq. (5) into Eq. (4), and assuming some simplifications and reorganizing according to powers of p-terms, it forms as

p0:ψ0′′(r)=0 (35)

p1:k2C1ψ0(r)−C1ψ01(r)r−ψ0′′−C1ψ0′′(r)+ψ1′′(r)=0 (36)

p2:k2C2ψ0(r)+k2C1ψ1(r)−C2ψ0′(r)r−C1ψ1′(r)r−C2ψ0′′−C1ψ1′′(r)−ψ1′′(r)+ψ2′′(r)=0 (37)

with the conditions

ψ0′(0)=0, ψ0(1)=1

we get from Eq. (35)

and similarly with the conditions

we get ψ1 from Eq. (36)

Now for ψ2(r), with the conditions

we get ψ2 from Eq. (37)

Hence the value of ψ(r) is determined from the Eq. ψ(r)=ψ0(r)+ψ1(r)+ψ2(r) is

ψ(r)=−k2(−1+r2)C1+124k2(−1+r2)(−24+k2(−5+r2))C12+12(2−k2(−1+r2)C2) (38)

By the method of least squares, the expression for residual becomes

Herewith k=2

J=∫01R2dr (39)

The values of constants C1 and C2 are obtained by setting

E1=∂J∂C1=0;E2=∂J∂C2=0 (40)

which give the following values of constants

Assuming the values of C1=−0.4503218 and C2=−0.1239597, based on these values, the calculated ψ is

ψ(k=2)=0.0337983(4(r2−5)−24)(r2−1)+1.80129(r2−1)+12(0.495839(r2−1)+2) (41)

Now for the Eqs. (18) and (19) with boundary conditions

with boundary conditions

with

u=u0(r,t)+pu1(r,t)+p2u2(r,t)+p3u3(r,t) (42)

and

v=v0(r,t)+qv1(r,t)+q2v2(r,t)+q3v3(r,t) (43)

According to OHAM for Eq. (42)

p0:∂u0∂r2(r,t)+A0+A1cos(α2t)=0 (44)

p1:C1Hα2u0(r,t)+C1Ru0(r,t)−C1Rv0(r,t)+C1∂u0∂t(r,t)−C1∂2u0∂r2(r,t)−C1∂u0∂r(r,t)r−∂2u0∂r2(r,t)+∂2u1∂r2(r,t)−A0−C1A0−A1C1cos(α2t)−A1cos(α2t)−C1kψ=0 (45)

p2:C2Hα2u0(r,t)+C1Hα2u1(r,t)+C2Ru0(r,t)−C1Ru1(r,t)−C2Rv0(r,t)−C1Rv1(r,t)+C2∂u0∂t(r,t)−C2∂2u0∂r2(r,t)−C2∂u0∂r(r,t)r+C1∂u1∂t(r,t)−C1∂2u1∂r2(r,t)−C1∂u1∂r(r,t)r−∂2u1∂r2(r,t)+∂2u2∂r2(r,t)−A0C2−A1C2cos(α2t)−C2kψ=0 (46)

Now the formulation of OHAM for Eq. (43)

q0:v0(r,t)−u0(r,t)=0 (47)

q1:−C1G∂v0∂t(r,t)+C1u0(r,t)−C1v0(r,t)+u0(r,t)−u1(r,t)−v0(r,t)+v1(r,t)=0 (48)

q2:−C2G∂v0∂t(r,t)−C1G∂v1∂t(r,t)+C2v0(r,t)+C1u1(r,t)−v1(r,t)−C2v0(r,t)−C1v1(r,t)+u1(r,t)−u2(r,t)+v2(r,t)=0 (49)

with the conditions

we get from Eq. (44)

now with the conditions

we get from Eq. (45)

for u2(r,t), with the following conditions

we get from Eq. (46)

now the equation

Becomes

Now using the relationship

we get

Again using the relationship

gives

Similarly with

we get

Now with the equation

gives