On the effect of mucus rheology on the muco-ciliary transport

Mathematical Biosciences 272 (2016) 44–53

Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

On the effect of mucus rheology on the muco-ciliary transport M.H. Sedaghat a, M.M. Shahmardan a,∗, M. Norouzi a, M. Nazari a, P.G. Jayathilake b a b

Department of Mechanical Engineering, University of Shahrood, Shahrood, Iran School of Mechanical & Systems Engineering, Newcastle University, Newcastle upon Tyne NE1 7RU, United Kingdom

a r t i c l e

i n f o

Article history: Received 26 April 2015 Revised 14 November 2015 Accepted 20 November 2015 Available online 1 December 2015 Keywords: Muco-ciliary transport Mucus Oldroyd-B model Finite difference-lattice Boltzmann method Immersed boundary method

a b s t r a c t A two dimensional numerical model is used to study the muco-ciliary transport process in human respiratory tract. Here, hybrid finite difference-lattice Boltzmann method is used to model the flow physics of the transport of mucus and periciliary liquid (PCL) layer in the airway surface liquid. The immersed boundary method is also used to implement the propulsive effect of the cilia and also the effects of the interface between the mucus and PCL layers. The main contribution of this study is on elucidating the role of the viscoelastic behavior of mucus on the muco-ciliary transport and for this purpose an Oldroyd-B model is used as the constitutive equation of mucus for the first time. Results show that the viscosity and viscosity ratio of mucus have an enormous effect on the muco-ciliary transport process. It is also seen that the mucus velocity is affected by mucus relaxation time when its value is less than 0.002 s. Results also indicate that the variation of these properties on the mucus velocity at lower values of viscosity ratio is more significant. © 2015 Published by Elsevier Inc.

1. Introduction The ventilation of human lung typically ranges between 1000 and 21,000 l a day depending on body size and physical activities [1]. Air which passes into the lungs from the surrounding environment is often polluted with varieties of particles such as fungi and bacteria which could be harmful for the lungs. As such, the airway surface liquid (ASL), which is a fluid layer coating the interior epithelial surfaces of the bronchi and bronchioles, plays an important defensive role against foreign particles and chemicals entering the lungs [2]. Indeed, the ASL exhibits a two-layer structure. The upper layer consists of highly viscous and non-Newtonian mucus which is a nonhomogeneous, viscoelastic fluid containing water, salt and glycosylated mucin proteins [3]. Lai et al. [4] in a valuable research reviewed the biochemistry that governs mucus rheology, the macro- and microrheology of human and laboratory animal mucus. They [4] clearly reveal the physical properties of mucus to advancing the field of drug and gene delivery. The lower layer of ASL is a watery lubricating fluid, or ‘periciliary liquid’ (PCL) layer. The origin of the PCL is not known, but is hypothesized to originate from osmosis across the epithelium [5]. The cilia, which are hair-like structures, are embedded in PCL with only the tips of the cilia contacting mucus layer. The bending of the cilium is produced by an internal arrangement of microtubules [6].

∗

Corresponding author. Tel.: +98 9124737057; fax: +98 2332300258 E-mail addresses: [email protected] (M.H. Sedaghat), [email protected], [email protected] (M.M. Shahmardan), [email protected] (M. Norouzi), [email protected] (M. Nazari), [email protected] (P.G. Jayathilake). http://dx.doi.org/10.1016/j.mbs.2015.11.010 0025-5564/© 2015 Published by Elsevier Inc.

Each cilium have two quite distinct stages to their beat cycle, the effective stroke when the cilium is extended out of the cell surface, and the recovery stroke when it creeps back to the start of its effective stroke. The effective stroke occupies a shorter period of the beat than the recovery stroke [7]. Hence, the work done during the effective stroke is several times the amount of work done performed during the recovery stroke [8]. It is hypothesized that the PCL serves both as a lubricant to allow movement of the viscous mucus layer and as a lower viscosity medium more favorable to the movement of cilia [9]. Considerable researches have been published in the field of modeling muco-ciliary transport. A quick review of the literature shows the Newtonian behavior of mucus has become the center of attention in the literature. The earliest analytical reference on the investigation concerning muco-ciliary flow was performed by Barton and Raynor [10]. They assumed cilium as a rigid rod, which automatically shortens during the recovery stroke. Use of this assumption would limit the accuracy of cilia motion. Lee et al. [11] used hybrid immersed boundary- finite difference method to simulate muco-ciliary transport process. In their study [11] the velocity of the PCL-mucus interface and the force exerted on the fluids by cilia are computed via the Immersed Boundary Method (IBM), which was originally introduced by Peskin [12] for studying blood flow patterns around heart valves [13] and since then has been used by many researches for simulating simple geometries [14–19] and also for biological problems [[11,20– 28]. Lee et al. [11] basically studied the dependency of mean mucus velocity on the mucus viscosity, the cilia beat frequency, the numbers of cilia, the thickness of PCL, and the surface tension at the interface between the mucus and PCL. Their results show that the cilia beat frequency, the number of cilia, and the depth of PCL are the

M.H. Sedaghat et al. / Mathematical Biosciences 272 (2016) 44–53

critical factors affecting the muco-ciliary transport. Jayathilake et al. [27] also used the IBM to develop a three-dimensional numerical model to simulate human pulmonary cilia motion in the PCL layer. Their results indicate the effects of the phase difference between cilia, the cilia beating frequency, the viscosity of PCL, the PCL height, and the ciliary length on the PCL motion. They concluded that the maximum PCL velocity in the stream-wise direction occurs if cilia have phase differences in both stream-wise and spanwise directions. Jayathilake et al. [28] in another numerical investigation studied various abnormalities of the cilia (e.g. cilia beat pattern, ciliary length, immotile cilia, beating amplitude and uncoordinated beating of cilia) on muco-ciliary transport of human respiratory tract. Their results show that the mucus velocity decreases by reducing the beat amplitude. They also represent that the windscreen wiper motion and rigid planar motion of cilia would decrease or almost stop the mucus transport. Chatelin and Poncet [29] modeled muco-ciliary clearance by using a hybrid scheme combining an Eulerian scheme for the Stokes equations and a particle method for transport equations. They [29] considered non-Newtonian aspects of viscosity in their method and study different mechanisms of mucus propelling and found dominant factors in muco-ciliary clearance, particularly in the context of cystic fibrosis and aerosol therapy. Similar investigations were performed to study mucus flow using similar methods [30–33]. Mauroy et al. [34] by considering mucus as a Bingham fluid (gel-like), studied the effects of geometry, mucus physical properties and amplitude of flow rate on the muco-ciliary clearance. Their results showed that both geometry and physical properties of mucus (yield stress and viscosity) play a vital role in mucus flow. It is readily acknowledged that mucus exhibit a variety of rheological complexities which cannot be described even qualitatively using the Newtonian fluid behavior and well-known Navier–Stokes equations. Experiments show that, rheological measurements including viscosity (resistance to flow) and elasticity (stiffness) are often used together to describe the consistency of mucus [4]. A review of the above mentioned researches clearly reveals that simulation of muco-ciliary transport has been investigated extensively considering the Newtonian behavior of mucus. In contrast, much less is known about the viscoelastic behavior of the mucus. Ross [35] in an analytical investigation considered mucus as a nonlinear Maxwell fluid. This study [35] has a few shortcomings that limit its general application. First, the PCL was not modeled in this study and mucus is modeled as a single layer. Secondly, the mucus–cilia interface was modeled as an impermeable wavy wall. Blake [36] stated that the latter assumption i.e. considering the tips of the cilia with a no-slip boundary is not appropriate because the tips of the cilia may penetrate the mucus layer during the effective stroke. As such, he [36] suggested a new approach in which the cilia could be modeled by distributing force singularities along their centerlines. This approach was also modified by some other researches later [37–39]. King et al. [40] used a simple analytical model to study the effect of mucus viscoelasticity on the muco-ciliary transport system. They concluded that mucus transport would intensify when the shear modulus of elasticity is reduced. They assumed that in their model there was no net transport of PCL in the cilia sublayer. However, the experiments of Matsui et al. [41] on human tracheobronchial epithelial culture showed that the transport of PCL approximately equals to the transport of mucus layer. King et al. [40] also assumed that there was a layer of PCL between the top of the cilia sublayer and the mucous layer although micrographs taken by Puchelle et al. [42] shows cilia penetrating into the mucus layer. Smith et al. [43] in an analytical investigation studied transport of mucus and periciliary liquid in the airways by considering mucus as a linearly viscoelastic fluid. In this work [43] the penetration of cilia into the mucus layer was modeled by considering the ASL as three fluid layers separated by flat interfaces. The lower PCL layer was modeled as a Newtonian fluid while the middle layer (traction layer), which represents the region of cilia penetration into the

45

mucus, was modeled by a Maxwell viscoelastic fluid with viscosity μM1 and relaxation time λ1 . The upper mucus layer was modeled by a Maxwell fluid of viscosity μM2 (which is greater than μM1 ) and the same relaxation time λ1 . In addition, the mat of cilia was modeled as an active porous medium. The propulsive effect of the cilia was modeled by time-dependent force acting in a shear-thinned traction layer between the mucus and the PCL. They studied various parameters on the motion of the mucus layer. Their results show that the dependency of mucus transport on the choice of physical parameters would be nonlinear. It was concluded that transport was only significantly disrupted by a reduction in cilia frequency. Their study [43] has two important shortcomings. First, they [43] considered the mucus layer as a linear viscoelastic fluid, while such a time-dependent problem in which the non-Newtonian viscosity plays a dominant role and hence the elastic effects must be also considered in the simulation. However, the elastic effects cannot be studied by linear viscoelastic models properly [44]. The other shortcoming of the model [43] is that they considered the mucus-PCL interface as completely flat. It is true that linear viscoelastic models have the advantage of being simple and convenient to implement in many computer modeling packages. However, although most biological materials are expected to behave linearly for small deformations (< 2–3%), biological materials typically exhibit nonlinear behavior at greater strains [45] and hence it would be important to include non-linear behaviors when modeling those biological systems. As the equations of linear viscoelasticity is not valid for large deformations since they do not satisfy the principle of frame invariance, Oldroyd and others developed a set of frame-invariant differential constitutive equations by defining time derivatives in frames that deform with the material elements. To extend the linear Maxwell model to the nonlinear regime, several time derivatives (e.g. upper convected, lower convected and corotational) are proposed to replace the ordinary time derivative in the original model. The idea of these derivatives is to express the constitutive equation in real space coordinates rather than local coordinates and hence fulfilling the Oldroyd’s admissibility criteria for constitutive equations. The Oldroyd-B model is a simplification of the more elaborate and rarely used Oldroyd 8-constant model which also contains the upper convected, the lower convected, and the corotational Maxwell equations as special cases. Oldroyd-B model is apparently the most popular in viscoelastic flow modeling and simulations [46]. In the current study, a computational model is used to simulate muco-ciliary transport process where mucus is considered as a viscoelastic fluid. The major thrust of this study is to make the role of rheological properties of mucus clearer in the muco-ciliary transport system by using an Oldroyd-B model as the constitutive equation of mucus. The Oldroyd-B model is discussed in term of the convected components of the stress tensor and the metric coefficients of the convected coordinate system. Material constants appeared in this constitutive equation and also temperature history could be included if it is desired to account for nonisothermal effects [44]. To the best our knowledge, this is the first attempt to introduce a proper model of the effects of the rheology properties of mucus in much greater details. The study of the selected properties is one of the most important steps in understanding the structure of mucus distribution in the lungs and to understand how mucus can be moved forward efficiently toward the trachea. In addition, it is believed that prescribing a drug for changing the rheological properties of mucus would be easier than changing the other properties like cilia frequency or length of cilia. We begin with by presenting the governing and constitutive equations describing this class of flows which is shown in Fig. 1 with appropriate boundary conditions in Sections 2 and 3. Section 4 describes the numerical calculation procedure for the muco-ciliary transport model in details. This is followed by a detailed description of the influence of rheological characteristic of mucus on the mucociliary transport problem in Section 5.

46

M.H. Sedaghat et al. / Mathematical Biosciences 272 (2016) 44–53

The first term σ M, N is related to Newtonian part of stress tensor and the second is related to the viscoelastic component of mucus. The Newtonian solvent contribution is given by

σM,N = 2ηM,N D

(3)

And, the viscoelastic contribution is derived based on the upper convected model (UCM) and is given below as ∇

σM,E + λ σM,E = 2ηM,E D

(4)

In these expressions, ηM, N is the viscosity of the Newtonian and ηE is the viscosity of the elastic part of mucus. The total viscosity and the viscosity ratio of mucus is given by Fig. 1. A schematic geometry for the muco-ciliary transport problem.

ηM = ηM,N + ηM,E ,β =

2. Governing equations Governing equations of incompressible fluid flow developed specifically for the geometry considered in this work (as shown in Fig. 1) can be written as follows [11]:

ηM,E . ηM

(5)

In Eq. (4), λ is the relaxation time, and the upper convected derivative of σ M, E is defined as ∇

σM,E =

∂σM,E  M,E − σM,E .∇ u u  .∇σ −∇  T .σM,E . +u ∂ t˜

(6)

 ·u =0 ∇

(1a)

The rate of deformation, D, is related to the velocity gradients as follows:

1 1 ∂ u u  .∇ + ∇ p= ∇ .σ + f +u ∂t ρ ρ

(1b)

D=

f (x, t ) =





 (s, t ))ds F (s, t )δ (x − X

∂ X (s, t )   = U (X (s, t ), t ) = ∂t

 

 (s, t ))dx,  (x, t )δ (x − X u

(1c)

(1d)

 is the velocity vector, ρ is denwhere x is the Eulerian coordinates, u sity, p is static pressure, t is time, σ is the stress tensor and f is the  is the Lagrangian coorboundary force acting on the fluid field. Also X  is the velocity vector of Lagrangian nodes, s is the arc length dinates, U of Lagrangian nodes and F is the boundary force density which con (s, t )) is a Dirac tains cilia and membrane forces. In addition, δ (x − X delta function. Eq. (1(c)) and (d) describes the interaction between the immersed boundaries and the fluid by distributing the boundary force at the Lagrangian points ( ) to the Eulerian points () and interpolating the velocity at the Eulerian points to the Lagrangian points. The geometry considered in the current study is shown in Fig. 1. As indicated in this figure the geometry consists of two layers; the bottom PCL layer where the cilia have a cyclic motion and the upper mucus layer. There is no interface material defined between the PCL and the mucus layers and as such the cilia can penetrate into the mucus layer through the mucus-PCL interface. In addition, the effect of the surface tension at the mucus-PCL interface is added by a force term generated due to an imaginary elastic membrane at the interface. In this study, the flow is assumed to be periodic in the horizontal direction while the bottom boundary is assumed to be a no-slip wall and the top boundary is assumed to be a free-slip boundary [11]. It should be mentioned that for ease of implementation it is assumed that the top boundary remains flat. However, the present numerical model can be easily extended to include a non-flat free surface of the mucus layer.

As remarked earlier, mucus exhibits the so-called viscoelastic fluid behavior. In this section, the constitutive equation of mucus as a viscoelastic fluid is explained by using the Oldroyd-B model. This model is useful for describing the rheological behavior of dilute viscoelastic solutions. Indeed, the stress tensor σ M can be decomposed into two parts as (2)

(7)

The Oldroyd-B constitutive equation can be derived from a molecular model in which the heavy molecule is idealized as an infinitely extensible Hookean spring connecting two Brownian beads [44]. This model has been specially used for predicting simple shear flows that are in qualitative agreement with measurements for certain Boger fluids [47]. 4. Numerical method The immersed boundary-lattice Boltzmann method (IB-LBM) is used to study the schematic geometry shown in Fig. 1. This method uses a fixed Cartesian mesh which is called Eulerian points to represent fluid phase. A set of Lagrangian points are used to represent the boundary immersed in the fluid [48]. In the present study, the Newtonian part of the momentum equation is solved by the lattice Boltzmann (LB) equation in Eulerian points and both the elastic part of the stress tensor and Lagrangian forces containing cilia and the interfacial forces are added to LB equations as the extra force terms. It is worth mentioning that the numerical method validations can be seen in [49] and would not repeated in the present paper. Physical properties of the fluid along with the dimension of the geometry have been defined in the physical domain. The lattice Boltzmann domain is also defined for solving the lattice Boltzmann equations. Scaling considerations are used to extract the pertinent relationship between physical and lattice Boltzmann domains. Each point on the physical grid is an equivalent of a point on the lattice Boltzmann domain. Therefore, transformation coefficients form the physical domain to the lattice Boltzmann domain should be introduced here to simulate real properties of the fluid and the real dimension of the geometry as follows:

Chx =

3. Constitutive equation

σM = σM,N + σM,E

 1 ∇ u + ∇ u T . 2

xphy yphy t phy ρ phy , Chy = , Ct = , Cρ = lb lb lb lb x y t ρ0

(8)

In Eq. (8) xphy and yphy are grid spacing in x and y directions, tphy is the time spacing and ρ phy is the density of the physical domain, xlb and ylb are grid spacing, tlb is the time spacing and ρ0lb is referenced to as density in LBM domain. In addition, Chx and Chy

are related to the space in x and y directions, respectively, Ct and Cρ are time and density transformation coefficients, respectively. In this study “lb” superscript is related to the lattice Boltzmann variables and for simplicity, physical variables are shown without any superscripts.

M.H. Sedaghat et al. / Mathematical Biosciences 272 (2016) 44–53

47

In Eq. (8), by considering x = y and xlb = ylb , we can write Chx = Chy = Ch . The other transformation coefficients can be obtained as follows:

 C 2 P νN Ch2 h = , C = = C ρ P Ct ν lb Ct P lb

 u C = h , Cν =  Ct ulb

Cu =

(9)

where ν N is the kinematic viscosity of the Newtonian part of the fluid and Cu ,Cν and CP are velocity, Newtonian kinematic viscosity and pressure transformation coefficients respectively. The lattice Boltzmann equations proposed by Guo et al. [50] are given in Eulerian points as follows:



fα x + eα t , t + t lb

lb

lb

lb





lb

− fα x , t

lb

  1   lb lb  eq =− fα x , t − fα xlb , t lb + Fα t lb τ  e − u lb e .ulb   1 α α eα · f lb Fα = ωα 1 − + 2τ cs2 cs4 ρ lb =

α

4.1. Cilia forces (10a)

(10b)

eα fα



ρ lb u lb =

α

eα fα +

(10d)

ρ lb lb 2

f t lb

(10d) eq



fα

=ρ

lb



2

2  lb  3u  lb )  lb ) 3(eα · u 9(eα · u ωα 1 + + − 2 4 c 2c 2c2

(11)

Here, the D2 Q9 lattice velocity model is used, in which the discrete velocity vectors are as follows:

⎧ [0, 0] α = 0 ⎪ ⎪      ⎪ ⎪ ⎨c cos (α − 1 )π , sin (α − 1 )π α= 1 − 4 eα =   2  2  ⎪ ⎪ √ (α − 5 )π π (α − 1 )π π ⎪ ⎪ + + α= 5 − 8 , sin ⎩ 2c cos 2 4 2 4 (12)

⎧4 ⎪ ⎪ 9 α= 0 ⎪ ⎨ ωα = 1 α = 1 − 4 9 ⎪ ⎪ ⎪ ⎩1 α= 5 − 8

 ∗lb = u

1

ρ lb

α

eα fα

(16)

And the velocity correction, proportional to the force term is defined as

1 lb lb f t 2

δ u lb =

(17)

Furthermore, the fluid velocity correction can be determined by interpolating the value of the cilia velocity correction

δ u lb =



  lb lb δU lb Di j xlb  sl i j − Xl

(18)

l

− Xllb ) is Delta function defined as [11] where Di j (xlb ij





lb lb lb lb Di j xlb = δ (xlb i j − Xl i j − Xl )δ (yi j − Yl )

(19a)

 δ (r ) =

1 − |r|, |r| ≤ 2 (19b)

0, |r| > 2

 lb , can be calculated as follows: The cilia velocity correction, δU



(20a)

δU1lb , δU2lb , . . . , δUmlb

T

(20b)

(13)

⎡

36

ν lb

(15)

δ u lb

In Eq. (15), is the intermediate fluid velocity and is the  ∗lb is given as velocity correction. The intermediate fluid velocity u

X=

Fig. 2 shows all possible velocities for D2Q9 model on the square lattice. The single relaxation time in LB equations (τ ) is calculated as

cs2 t lb

 ∗lb + δ u  lb = u  lb u

 =B  AX

And the weighting coefficients wi are defined as

τ=

As mentioned before, cilia force is the force which cilia impose on the fluid and moves it forward. Immersed boundary method is used to calculate this force [48]. Here, the interpolated velocity of fluid at the cilia points is enforced to be equal to the velocity of the cilia at the same position at every evolution time step by a set of velocity corrections. Indeed, the velocity of the fluid field is defined as follows:

 ∗lb u

where fα is the distribution function, fα is its corresponding equilibrium distribution function for the discrete velocity eα , Fα is the discrete force term, f lb is the force density action on the fluid field, wα is the weighting coefficient which depends on the selected lattice vesingle relaxation time, ρ lb is the density of fluid locity model, τ is the √ in LB domain, cs = c/ 3 is the sound speed and c = xlb /t lb is the lattice speed where xlb and t lb are the lattice constant and the eq time step size, respectively. The equilibrium distribution function fα is defined as eq

Fig. 2. Lattice arrangements for 2D problems, D2Q9.



+ 0.5

(14)

where ν lb is the kinematic viscosity in LBM domain. lb which is reIn Eq. (10), f lb consists of three different forces: fCilia lated to the force that each cilium impose on the fluid, f lb which is Mem

related to the elastic force on the interface between the mucus and PCL and fElb which is related to the elastic part of the stress tensor in viscoelastic fluids. These forces are explained in details below.

δ11 δ12 · · · ⎢ δ21 δ22 · · · A=⎢ . .. .. ⎣ .. . . δm1 δm2 · · · ⎧  lb ⎫ ⎡ U1 ⎪ δ11 ⎪ ⎪ ⎨U lb ⎪ ⎬ ⎢ δ21 2 B= −⎢ . .. ⎣ .. ⎪ ⎪ . ⎪ ⎩  lb ⎪ ⎭ δm1 Um

⎤⎡ B ⎤ B B · · · δ1m δ11 δ12 δ1n B B B δ21 ⎥⎢δ21 δ22 · · · δ2m ⎥ ⎢ .. ⎥ .. .. ⎥ .. ⎦⎣ ... ⎦ . . . . B B B δmn δn1 δn2 · · · δnm ⎤⎧ ⎫ δ12 · · · δ1n ⎪u ∗lb ⎪ 1 ⎪ ⎪ ⎬ δ22 · · · δ21 ⎥⎨u ∗lb 2 ⎥ .. . . .. .. ⎦⎪ .. ⎪ . . ⎪ ⎩ ∗lb ⎪ ⎭ δm2 · · · δmn n u

(20c)

where, m is the number of the Lagrangian points on the cilia, and  lb (l = 1, 2, . . . , m ) is the n is the number of the Eulerian points. U l  lb (l = 1, 2, . . . , m ) velocity vector of Lagrangian (cilia) points and δU l

48

M.H. Sedaghat et al. / Mathematical Biosciences 272 (2016) 44–53

is the unknown velocity correction vector at the Lagrangian points. − Xllb )xlb ylb and δiBj = Di j (xlb − Xllb )slb In Eq. (20), δi j = Di j (xlb ij ij l

where slb is the arc length of cilia elements. Therefore, according to l Eq. (17) the relationship between the force density, FCilia and the fluid  lb , can be written as velocity correction, δU

 lb 2δU lb = FCilia δt

(21)

By interpolating the value of cilia forces from Lagrangian (cilia) nodes to Eulerian (fluid) points, the cilia force on each node of the fluid can be calculated as

f lb = Cilia



lb lb lb Di j (xlb FCilia i j − Xl sl )

(22)

l

It is noteworthy that all the Lagrangian variables in this part (e.g.  lb ) are related to the cilia.  lb and δU Xllb ,slb ,U l l 4.2. Membrane force The elastic force which represents the surface tension exerted by the membrane interface of PCL and mucus on the fluid is calculated as [11]

∂ FMem = [T (s )τ (s )] ρ ∂s 1

(23)



∂ XMem T (s ) = T0 −1 ∂ s0   ∂ XMem  ∂ XMem  / τ (s ) = ∂s  ∂s 

Ct 2 Ch



lb lb lb Di j (xlb FMem i j − XMem,l sMem,l )

(29)

The interfacial velocity UMem (XMem ) is interpolated from the neighboring grid points. Subsequently, the motion of the membrane interface is updated as (30)

In summary, the numerical implementation of the present scheme can be outlined as follows:

(24b)

i. Knowing the values of physical geometry, x,t and by setting xlb = ylb = t lb = 1 space, time and density transformation coefficients are calculated from Eq. (8). ii. Knowing the Newtonian kinematic viscosity and velocity in all nodes and using velocity and viscosity transformation coefficients (Eq. (9)), kinematic viscosity and velocity are computed in all nodes in LBM domain. iii. Using Eq. (10(a)) to get the pressure distribution function at time level t = tn (initially setting Fα = 0), the macroscopic variables are computed using Eq. (10(c)) and (10(d)). iv. Using velocity transforms coefficient and changing the obtained velocity from Step 3 to physical velocity, the total force term ( f lb ) is obtained from Eq. (29). v. Fα is computed from Eq. (10(b)) and the new velocity is calculated from Eq. (10(d)) at t = tn + 1. vi. The position of the membrane interface is updated as Eq. (30) by finding UMem at t = tn + 1 .

(25)

lb ) from LaBy interpolating the value of interfacial force (FMem grangian (interfacial) nodes to Eulerian (fluid) points, the interfacial force on each node of the fluid can be calculated as

f lb = Mem

f lb = f lb + f lb + f lb Mem E Cilia

(24a)

In Eq. (24(a)) and (b), τ (s) is the unit tangential vector to the PCL– mucus interface. The arc-lengths s and s0 are measured along the current and initial configuration of the interface, respectively. The scalar T0 is the stiffness constant which describes the elastic property of the flexible boundary. Since FMem is defined as the physical variables, it should be multiplied by force transformation coefficient to be transformed into LB force as follows: lb FMem = FMem

The value of the total force which should be added to LB equations as the extra force term is given by

∂ XMem = UMem ∂t

where T(s) and τ (s) are defined as



Fig. 3. The beat cycle derived by Fulford and Blake [37]

(26)

l

5. Results and discussion

4.3. Elastic forces As mentioned in Section 3, the stress tensor of mucus layer,σ M, E , is decomposed into Newtonian and elastic parts. The Newtonian part of σ is solved by LBM as described in the numerical method. Using the physical velocity obtained from LB equations, the elastic part of the stress tensor (Eq. (4)) is solved by finite difference method and the value of σ M, E in each time step is calculated. Therefore, the elastic force of the stress tensor can be calculated as

fE = 1 ∇  .σM,E

ρ

(27)

Like fMem , fE is also defined as the physical variables, so it should be multiplied by force transformation coefficient to be transformed into LB force as follows: 2

f lb = fE Ct E Ch

(28)

The effects of rheological characteristic of mucus on the mucociliary clearance are studied and the numerical results are presented here. Indeed, by considering an Oldroyd-B model as the constitutive equation of mucus, the effects of viscosity ratio (β = ηM,E /ηM ), relaxation time (λ) and viscosity (ηM ) of mucus on the muco-ciliary clearance are studied in details. At the beginning of this section, it is necessary to define the geometry in details and mention the standard parameters (parameters commonly appeared in the literature for a healthy human respiratory system). Cilia have a cyclic motion in PCL as also mentioned in the preceding discussion. In order to calculate the cilia forces, it is straightforward to postulate the cilia beat pattern in each time step. Sanderson and Sleigh [51] observed the ciliary beat cycle in the rabbit trachea using high speed cine-photography and scanning electron microscopy. Folfurd and Blake [37] used Sanderson and Sleigh’s data [51] and

M.H. Sedaghat et al. / Mathematical Biosciences 272 (2016) 44–53

Fig. 4. Arrangement of the cilia for t = 0 and 8Tc /13 in the present study.

49

proposed an analytical representation by the truncated Fourier series for the cilia beat pattern. The 2D cilium beat cycle in this study is the cilia beat pattern reported in Folfurd and Blake [37] at 13 different time step as shown in Fig. 3. In the present study, 13 cilia each having a cyclic motion are considered. The cilia arrangements for t = 0 and 8Tc /13 are plotted in Fig. 4 (Tc is the consecutive beat cycle). The length of cilia is LCilia = 6 μm [43] and the spacing between any two neighboring cilia along the wall is chosen as d = 3 μm [11]. Pedersen et al. [52] declared the depth of the PCL was about 5–10 μm in the upper airways of bronchi, trachea and nasal cavity. In this study, the standard depth of PCL is so considered as LPCL = 6 μm and the depth of mucus layer is as LM = 4 μm [43]. The viscosity of PCL is ηPCL = 0.001Pa.s [43] and the viscosity of mucus layer is ηM = 0.0482 Pa.s[43]. As remarked earlier, ηM is decomposed into viscous and elastic parts. Since a full viscoelastic characterization for mucus is not yet available, it is assumed that the standard viscosity of viscous part of mucus is as the same as the viscosity of PCL i.e. ηM,N = 0.001Pa.s and the elastic part of mucus viscosity is ηM,E = 0.0472 Pa.s. Therefore, the viscosity ratio is β = ηM,E /ηM ≈ 0.98. The standard value of relaxation time of mucus is about λ = .034s [43] and the density of the fluid is ρ = 1000 kg/m3 [43]. The stiffness constant of the membrane is T0 = 32 × 10−3 N/m [11] and the cilia beat frequency is set at σ = 60 rad/s [43]. The computational mesh is chosen as 112 × 29 (x = y ≈ 3.5 × 10−7 ), the number of control points along each cilium is 20 and along mucus-PCL interface is 105. The time step is t = 8.5 × 10−7 s. Our

Fig. 5. The velocity field and the beating cilia at t = 0 and 4Tc /13 for the standard parameter set.

50

M.H. Sedaghat et al. / Mathematical Biosciences 272 (2016) 44–53

Fig. 6. Enlarged location of the mucus-PCL interface at t = 8Tc /13 for the standard parameter set.

Fig. 7. Non-dimensional mean mucus transport against viscosity ratio β for λ = 0.034s and ηM = 0.0482 Pa.s. Table 1 Mean mucus velocity (μm/s) for various mucus relaxation times and viscosity ratios for ηM = 0.0482 Pa.s.

β /λ(s)

0.001

0.002

0.005

0.01

0.02

0.034

0.1

0. 5

1

2

0.1 0.3 0. 5 0.7 Standard parameters

307.6 161.2 – – –

356.8 240.9 147.8 122.9 43.18

405.9 306 156.8 132.4 43.78

406.6 331.2 243.1 165.5 43.96

406.9 331.8 253.8 171.9 44.04

407.3 332.3 254.3 172 44.07

407.3 332.4 254.3 172 44.07

407.3 332.5 254.3 172 44.07

407.3 332.5 254.4 172 44.07

407.3 332.5 254.4 172 44.07

simulation using the standard parameter set predicts a time average mucus velocity of u0 = 44.07 μm/s where u0 is calculated as

u0 =

Fig. 8. Variation of non-dimensional mean mucus velocity as against mucus relaxation time at various viscosity ratio β .

1 Tc (H2 − H1 )

Tc H2 uM dydt

(31)

0 H1

H2 − H1 in Eq. (31), is the depth of mucus layer and Tc is the consecutive beat cycle. As remarked earlier, using the standard parameter set our simulation predicts a mean mucus velocity of 44.07 μm/s. ICRP [53] reported a wide range of values depending upon disease, ambient conditions and other factors. For healthy subjects, values of 70 and 92 μm/s for tracheal transport, and 40 μm/s for bronchial transport were reported. The experimental investigation performed by Matsui et al. [41] showed a mean mucus transport of 39.2 μm/s. Smith et al. [43] showed a mean mucus velocity of 38.3 μm/s. Also, numerical simulation of Lee et al. [11] and Jayathilake et al. [28] predicts a mean

M.H. Sedaghat et al. / Mathematical Biosciences 272 (2016) 44–53

51

Table 2 Mean mucus velocity (μm/s) for various mucus viscosities and viscosity ratios for λ = 0.034 s.

β /ηM (Pa.s)

0.002

0.005

0.01

0.02

0.0482

0.06

0.08

0. 1

0.2

0 0.1 0. 3 0.5 0.7 Standard parameters

56.21 53.89 49.08 43.95 33.91 44

86.50 81.70 71.91 61.69 50.38 44

132.2 123.2 104.9 86.49 68.17 44

219.5 202.15 167.4 132.2 95.97 44.03

443.3 407.4 322.6 254.3 172 44.07

525.3 483.8 395.9 302.3 202.4 44.07

646.5 598.9 495.9 376.6 252.3 44.07

744.1 695.6 583.5 440 301 44.07

934.8 920.9 856 511.5 382.1 44.07

Fig. 9. Non-dimensional mean mucus velocity as a function of mucus viscosity at various viscosity ratio β .

mucus velocity of 44.38 μm/s using the standard parameter set. Evidently, the value predicted in this study is in a reasonable agreement with results reported in the literature. The velocity field and the beating cilia at two different times for the standard parameter set are plotted in Fig. 5. This figure shows a high forward velocity around the cilium at its effective stroke where the largest mucus velocity occurs close to the mucus-PCL interface. Fig. 6 illustrates enlarged location of the mucus – PCL interface at t = 8Tc /13. The temporal configuration of the membrane is affected by the mucus and PCL velocities and also the stiffness constant of the membrane. Fig. 7 demonstrates the non-dimensional mean mucus velocity versus the viscosity ratio. This figure indicates that the viscosity ratio has a great effect on mucus transport. The mean mucus velocity decreases rather linearly as β increases. It is apparent that the mean mucus velocity doubles as the viscosity ratio reduces from standard condition (β = 0.98) to β = 0.9. As Fig. 7 also indicates the mucus velocity is about ten times the standard mucus velocity when β = 0 which means that by reducing the elastic part or by increasing the Newtonian part of mucus viscosity, mucus velocity can increase significantly. This is partly due to chain breakage in the system which leads to decrease in the molecular weight of viscoelastic fluids [54]. Because of the pronounced effect of the viscosity ratio on mucociliary transport in this study, the effects of other rheological properties such as relaxation time and mucus viscosity are also investigated at different values of β . The variation of mucus velocity as a function of the relaxation time (λ) at various viscosity ratios (β ) is summarized in Table 1. This table reveals that mean mucus velocity increases by increasing λ at given β . This is primarily due to keeping the mucus in a strained condition and hence causing some amount of plastic strain and as a result, the mucus becomes stiffer and then the PCL can move the mucus more easily. Smith et al. [43] has also discussed briefly on this finding. The last row of this Table indicates that, increasing the value of λ have no great effect on the mucus velocity at standard parameter values. To better observe the effects of relaxation time on mucus velocity, Fig. 8 shows the variation of

non-dimensional mucus velocity as a function of λ. This figure demonstrates that variation of mucus velocity is almost negligible at higher values of relaxation time. This could be due to the fact that at higher values of relaxation time, mucus is elastic enough and PCL can move it easily. This figure also shows that the variation of mucus velocity is more significant at lower values of β and λ. This is because both viscous and elastic components of mucus play an important role in the mucus flow. By increasing β or λ, the elastic effects of fluid overcomes the viscous effect and hence the variation of mucus velocity almost remains constant. Table 2 shows the variation of mucus velocity with respect to mucus viscosity at different values of viscosity ratio. This Table indicates that mucus viscosity at lower values of β has more significant effect on the mucus velocity. For instance, at β = 0, mucus velocity changes from 56.21 to 934.8 for ηM ranging from 0.002 to 0.2, while at β = 0.7, mucus velocity changes from 33.91 to 382.1 for the same range of ηM . Last row of this Table indicates that at the standard condition, the variation of ηM have no a pronounced effect on the mucus velocity. For better understanding of the discussion, the variation of nondimensional mucus velocity has been plotted in Fig. 9 as a function of β and ηM . This figure also indicates that in addition to β , ηM is another effective parameter which influences mucus velocity. For instance, at β = 0 when the value of ηM is 0.2, mucus velocity reaches about 21 times the standard mucus velocity u0 . This is possibly due to the fact that by increasing the mucus viscosity, the mucus become stiffer and the PCL can move it forward more easily. Increasing the mucus velocity by increasing the viscosity of mucus also reported by Smith et al. [43]. Fig. 9 also shows that by increasing β , the variation of mucus velocity decreases with respect to mucus viscosity. It means that at higher values of β (e.g. β = 0.5, 0.7) the mucus velocity almost remains constant as ηM increases (for ηM ≥ 0.1). In other words, at higher values of β and ηM , increasing mucus viscosity, have no a great effect on the mucus velocity. This is partly due to the fact that by increasing β , the elastic part of the mucus layer intensifies, and hence increase of mucus viscosity cannot overcome its elastic properties. Fig. 9 also demonstrates that the influence of

52

M.H. Sedaghat et al. / Mathematical Biosciences 272 (2016) 44–53

β on the mean mucus velocity intensifies by increasing the value of ηM . It is noteworthy to mention that in the present numerical investigations, it is assumed that the cilia would not penetrate into the mucus layer and hence increase the mucus velocity with mucus viscosity might be due to this particular assumption. However, in real biological systems, the tips of cilia enter the mucus layer during their effective stroke and increase in mucus viscosity may decrease the cilia motion and also the mucus velocity. 6. Conclusions In the present study, muco-ciliary transport process of human respiratory system is numerically simulated considering mucus as a viscoelastic fluid. The two layer model is solved by using immersed boundary-lattice Boltzmann method. The cilia beat is prescribed by using the beat pattern reported in Folfurd and Blake [37]. An OldroydB model is employed to study the detailed effects of mucus rheological properties such as the mucus viscosity ratio, relaxation time and viscosity on the mean mucus velocity. The results indicate that: 1. Viscosity ratio has an immense effect on muco-ciliary transport and when reducing the viscosity ratio, the mucus velocity increases linearly. For example, at β = 0, mucus velocity is about ten times the standard mucus velocity. In other words, the mucus velocity can be increased either by reducing the elastic component of the mucus viscosity or by increasing its Newtonian component. 2. Mucus viscosity plays an important role in the mucus velocity. This effect is more considerable at lower values of β and almost negligible at higher values of β when ηM ≥ 0.1. Our results also indicate that the influence of β on the mean mucus velocity intensifies at higher values of ηM . 3. The effect of relaxation time of mucus is not as significant as mucus viscosity and viscosity ratio on the mucus flow, but especially at lower values of β , it can be an important factor to change the mucus velocity. It should be noted however that by increasing mucus relaxation from its standard value, the increase of λ does not have no effect on the mucus velocity. References [1] A. Wanner, M. Salathé, T.G. O’Riordan, Mucociliary clearance in the airways, Am. J. Resp. Crit. Care. 154 (1996) 1868–1902. [2] E.P. Lillehoj, K.C. Kim, Airway mucus: its components and function, Arch. Pharm. Res. 25 (2002) 770–780. [3] D. Smith, E. Gaffney, J. Blake, Modelling mucociliary clearance, Respir. Physiol. Neurobiol. 163 (2008) 178–188. [4] S.K. Lai, Y.-Y. Wang, D. Wirtz, J. Hanes, Micro-and macrorheology of mucus, Adv. Drug. Deliv. Rev. 61 (2009) 86–100. [5] A. Dauptain, J. Favier, A. Bottaro, Hydrodynamics of ciliary propulsion, J. Fluid Struct. 24 (2008) 1156–1165. [6] M. Murase, The Dynamics of Cellular Motility, John Wiley, Chichester, 1992. [7] J. Blake, Hydrodynamic calculations on the movements of cilia and flagella I, J. Theor. Biol. 45 (1974) 183–203. [8] S. Gueron, K. Levit-Gurevich, Energetic considerations of ciliary beating and the advantage of metachronal coordination, in: Proceedings of the National Academy of Sciences, 96, 1999, pp. 12240–12245. [9] S.M. Mitran, Metachronal wave formation in a model of pulmonary cilia, Comput. Struct. 85 (2007) 763–774. [10] C. Barton, S. Raynor, Analytical investigation of cilia induced mucous flow, Bull. Math. Biophys 29 (1967) 419–428. [11] W. Lee, P. Jayathilake, Z. Tan, D. Le, H. Lee, B. Khoo, Muco-ciliary transport: effect of mucus viscosity, cilia beat frequency and cilia density, Comput. Fluids 49 (2011) 214–221. [12] C.S. Peskin, The immersed boundary method, Acta Numer 11 (2002) 479–517. [13] C.S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 (1977) 220–252. [14] D.-V. Le, B.C. Khoo, J. Peraire, An immersed interface method for viscous incompressible flows involving rigid and flexible boundaries, J. Comput. Phys. 220 (2006) 109–138. [15] D. Le, B. Khoo, K. Lim, An implicit-forcing immersed boundary method for simulating viscous flows in irregular domains, Comput. Methods Appl. Mech. 197 (2008) 2119–2130.

[16] Z. Tan, D.-V. Le, Z. Li, K. Lim, B.C. Khoo, An immersed interface method for solving incompressible viscous flows with piecewise constant viscosity across a moving elastic membrane, J. Comput. Phys. 227 (2008) 9955–9983. [17] Z. Tan, D.-V. Le, Z. Li, K. Lim, B.C. Khoo, An immersed interface method for solving incompressible viscous flows with piecewise constant viscosity across a moving elastic membrane, J. Comput. Phys. 227 (2008) 9955–9983. [18] Z. Wang, J. Fan, K. Cen, Immersed boundary method for the simulation of 2D viscous flow based on vorticity–velocity formulations, J. Comput. Phys. 228 (2009) 1504–1520. [19] Z. Wang, J. Fan, K. Luo, Combined multi-direct forcing and immersed boundary method for simulating flows with moving particles, Int. J. Multiph. Flow 34 (2008) 283–302. [20] A. Dauptain, J. Favier, A. Bottaro, Hydrodynamics of beating cilia, in: Proceedings of IUTAM Symposium on Unsteady Separated Flows and their Control, Corfu, Greece, 2007, pp. 18–22. [21] A. Dauptain, J. Favier, A. Bottaro, Hydrodynamics of ciliary propulsion, J. Fluid. Struct 24 (2008) 1156–1165. [22] R. Dillon, L. Fauci, X. Yang, Sperm motility and multiciliary beating: an integrative mechanical model, Comput. Math. Appl. 52 (2006) 749–758. [23] R.H. Dillon, L.J. Fauci, An integrative model of internal axoneme mechanics and external fluid dynamics in ciliary beating, J. Theor. Biol. 207 (2000) 415– 430. [24] R.H. Dillon, L.J. Fauci, C. Omoto, X. Yang, Fluid dynamic models of flagellar and ciliary beating, Ann. N. Y. Acad. Sci. 1101 (2007) 494–505. [25] X. Yang, R.H. Dillon, L.J. Fauci, An integrative computational model of multiciliary beating, Bull. Math. Biol. 70 (2008) 1192–1215. [26] D.M. McQueen, C.S. Peskin, Heart simulation by an immersed boundary method with formal second-order accuracy and reduced numerical viscosity, Mechanics for a New Mellennium, Springer, 2002, pp. 429–444. [27] P.G. Jayathilake, Z. Tan, D. Le, H. Lee, B. Khoo, Three-dimensional numerical simulations of human pulmonary cilia in the periciliary liquid layer by the immersed boundary method, Comput. Fluids 67 (2012) 130–137. [28] P.G. Jayathilake, D. Le, Z. Tan, H. Lee, B. Khoo, A numerical study of muco-ciliary transport under the condition of diseased cilia, Comput. Methods Biomech 18 (2014) 944–951. [29] R. Chatelin, P. Poncet, A hybrid grid-particle method for moving bodies in 3D Stokes flow with variable viscosity, SIAM J. Sci. Comput 35 (2013) B925– B949. [30] R. Chatelin, P. Poncet, Hybrid grid–particle methods and penalization: a Sherman–Morrison–Woodbury approach to compute 3D viscous flows using FFT, J. Comput. Phys. 269 (2014) 314–328. [31] R. Chatelin, P. Poncet, M. Tokman, Computational aspects of mucus propulsion by ciliated epithelium, in: Proceedings of 2nd European Microfludics Conference, Toulouse, 2010. [32] R. Chatelin, P. Poncet, Particle methods for 3D biological flows with variable density and viscosity, in: 6th ECCOMAS Conference, Vienna, Austria, 2012. [33] R. Chatelin, P. Poncet, A. Didier, M. Murris-Espin, D. Anne-Archard, M. Thiriet, Mucus and ciliated cells of human lung: splitting strategies for particle methods and 3D stokes flows, in: IUTAM Symposium on Particle Methods in Fluid Mechanics, Lingby, Denmark, 2012. [34] B. Mauroy, C. Fausser, D. Pelca, J. Merckx, P. Flaud, Toward the modeling of mucus draining from the human lung: role of the geometry of the airway tree, Phys. Biol 8 (2011) 056006. [35] S.M. Ross, A Wavy Wall Analytical Model Of Muco-Ciliary Pumping, Ph.D thesis, Johns Hopkins University, Baltimore, Md, 1971. [36] J. Blake, A model for the micro-structure in ciliated organisms, J. Fluid. Mech. 55 (1972) 1–23. [37] G. Fulford, J. Blake, Muco-ciliary transport in the lung, J. Theor. Biol 121 (1986) 381–402. [38] J. Blake, H. Winet, On the mechanics of muco-ciliary transport, Biorheology 17 (1980) 125. [39] N. Liron, S. Mochon, The discrete-cilia approach to propulsion of ciliated microorganisms, J. Fluid. Mech. 75 (1976) 593–607. [40] M. King, M. Agarwal, J. Shukla, A planar model for mucociliary transport: effect of mucus viscoelasticity, Biorheology 30 (1992) 49–61. [41] H. Matsui, S.H. Randell, S.W. Peretti, C.W. Davis, R.C. Boucher, Coordinated clearance of periciliary liquid and mucus from airway surfaces, J. Clin. Investig. 102 (1998) 1125. [42] E. Puchelle, A.-L. Herard, J.-M. Zahm, Airway Mucociliary Epithelium Injury and Repair, Cilia, Mucus, and Mucociliary Interactions, Dekker, New York, 1998, pp. 203–217. [43] D.J. Smith, E. Gaffney, J. Blake, A viscoelastic traction layer model of muco-ciliary transport, Bull. Math. Biol. 69 (2007) 289–327. [44] R.B. Bird, R.C. Armstrong, O. Hassager, C.F. Curtiss, Dynamics of Polymeric Liquids, Wiley, New York, 1977. [45] J. Funk, G. Hall, J. Crandall, W. Pilkey, Linear and quasi-linear viscoelastic characterization of ankle ligaments, J. Biomech. Eng. T ASME 122 (2000) 15– 22. [46] T. Sochi, Pore-Scale Modeling of Non-Newtonian Flow in Porous Media, Imperial College London, 2007 Department of Earth Science and Engineering. [47] M.E. Mackay, D.V. Boger, An explanation of the rheological properties of Boger fluids, J. Non-Newton. Fluid. 22 (1987) 235–243. [48] J. Wu, C. Shu, Particulate flow simulation via a boundary condition-enforced immersed boundary-Lattice Boltzmann scheme, Comm. Comp. Phys. 7 (2010) 793.

M.H. Sedaghat et al. / Mathematical Biosciences 272 (2016) 44–53 [49] M.H. Sedaghat, M.M. Shahmardan, M. Nazari, M. Norouzi, Immersed boundarylattice Boltzmann method for modeling non-Newtonian fluid flow around curved boundaries, Modares Mech. Eng. 14 (8) (2014) 146–156 (In Persian). [50] Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev. E 65 (2002) 046308. [51] M. Sanderson, M. Sleigh, Ciliary activity of cultured rabbit tracheal epithelium: beat pattern and metachrony, J. Cell. Sci. 47 (1981) 331–347.

53

[52] P.S. Pedersen, N.-H. Holstein-Rathlou, P.L. Larsen, K. Qvortrup, O. Frederiksen, Fluid absorption related to ion transport in human airway epithelial spheroids, Am. J. Physiol. Lung C. 277 (1999) L1096–L1103. [53] I.C.o.R. Protection, ICRP, ICRP Publication 66: Human Respiratory Tract Model for Radiological Protection, Elsevier Health Sciences, 1995. [54] M. Yamamoto, The visco-elastic properties of network structure II. Structural viscosity, J. Phys. Soc. Jpn 12 (1957) 1148–1158.