Robin Chatelin

Philippe Poncet

A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus

Keywords: Mucociliary clearance, Cilia, Mucus, 3D computations of non-homogeneous flows, Fluid-structure interaction in complex geometry

Mucociliary clearance is the natural flow of the mucus which covers and protects the lung from the outer world. Pathologies, like cystic fibrosis, highly change the biological parameters of the mucus flow leading to stagnation situations and pathogens proliferation. As the lung exhibits a complex dyadic structure, in-vivo experimental study of mucociliary clearance is almost impossible and numerical simulations can bring important knowledge about this biological flow. This article brings a detailed study of the biological parameters influence on the mucociliary clearance, in particular for pathological situations such as cystic fibrosis. Using recent suitable numerical methods, a non-homogeneous mucus flow (including non-linearities) can be simulated efficiently in 3D, allowing the identification of the meaningful parameters involved in this biological flow. Among these parameters, it is shown that the mucus viscosity, the stiffness transition between pericilliary fluid and mucus, the pericilliary fluid height as well as both cilia length and beating frequency have a great influence on the mucociliary transport.

Introduction

In the lung, bronchial walls are covered by a physiological fluid: the airway surface liquid (ASL). It plays a major role to protect bronchus from the outer word: inhaled pathogens, dust and pollution particles are trapped in this viscous fluid film which is renewed to prevent contaminations. The ASL is constantly flowing form distal airways to trachea thanks to cilia located on the bronchial wall: they are beating at high frequency (4-10 Hz) to propel the thin film (10-15µm [START_REF] Fahy | Airway mucus function and dysfunction[END_REF]) which is swallowed in the stomach. This natural flow is called the mucociliary clearance (MCC). When it fails, the cough tries to overcome the ASL accumulation to prevent pathogens proliferation and contaminations.

Essentially composed of water, the ASL was shown to be composed of two different layers [START_REF] Lucas | Principles underlying ciliary activity in the respiratory tract: Ii. a comparison of nasal clearance in man, monkey and other mammals[END_REF][START_REF] Yoneda | Mucous Blanket of Rat Bronchus[END_REF][START_REF] Sanderson | Ciliary activity of cultured rabbit tracheal epithelium: beat pattern and metachrony[END_REF]: close to the bronchial wall, the "pericilliary fluid" layer (PCL) has a low viscosity (similar to water) whereas close to the air interface, the "mucus" layer (ML) is much more viscous. This difference is due to proteins: the mucins, which tend to accumulate close to the air-mucus interface and increase the viscosity of the fluid.

The MCC is very difficult to study experimentally since in vivo measurements are almost impossible. This problems comes from the complex structure of the tracheobronchial tree and the lack of non invasive experimental devices (a recurrent problem in biomechanics). Hence, numerical simulations can bring important informations to better understand the mucus flow, particularly under pathological situations. 1 Corresponding author: robin.chatelin@enise.fr Since the work of Blake in the 1970's [START_REF] Blake | A model for the micro-structure in ciliated organisms[END_REF] several researches have brought new insights on mucociliary clearance by using numerical computations -see [START_REF] Smith | Modelling mucociliary clearance[END_REF] for a review. Since the date of this review (2008), several works presented computations using 2D [START_REF] Smith | A viscoelastic traction layer model of muco-ciliary transport[END_REF][START_REF] Lukens | Using Lagrangian coherent structures to analyze fluid mixing by cilia[END_REF][START_REF] Montenegro-Johnson | Physics of rheologically enhanced propulsion: Different strokes in generalized stokes[END_REF][START_REF] Lee | Muco-ciliary transport: Effect of mucus viscosity, cilia beat frequency and cilia density[END_REF][START_REF] Jayathilake | A numerical study of muco-ciliary transport under the condition of diseased cilia[END_REF][START_REF] Sedaghat | On the effect of mucus rheology on the muco-ciliary transport[END_REF] or 3D modeling [START_REF] Jayathilake | Three-dimensional numerical simulations of human pulmonary cilia in the periciliary liquid layer by the immersed boundary method[END_REF][START_REF] Simons | A fully three-dimensional model of the interaction of driven elastic filaments in a Stokes flow with applications to sperm motility[END_REF][START_REF] Smith | Discrete cilia modelling with singularity distributions: Application to the embryonic node and the airway surface liquid[END_REF][START_REF] Li | An immersed boundary-lattice Boltzmann method for single-and multi-component fluid flows[END_REF]. Among these publications, some works are assuming a onelayer ASL with constant viscosity [START_REF] Jayathilake | Three-dimensional numerical simulations of human pulmonary cilia in the periciliary liquid layer by the immersed boundary method[END_REF][START_REF] Simons | A fully three-dimensional model of the interaction of driven elastic filaments in a Stokes flow with applications to sperm motility[END_REF] whereas others focused on a two-layer modeling of the ASL with constant viscosity in each layers [START_REF] Lee | Muco-ciliary transport: Effect of mucus viscosity, cilia beat frequency and cilia density[END_REF][START_REF] Jayathilake | A numerical study of muco-ciliary transport under the condition of diseased cilia[END_REF][START_REF] Li | An immersed boundary-lattice Boltzmann method for single-and multi-component fluid flows[END_REF] or integrating a more complex rheology in the ML (rheological parameters remain constant in each layer) [START_REF] Mitran | Metachronal wave formation in a model of pulmonary cilia[END_REF][START_REF] Smith | A viscoelastic traction layer model of muco-ciliary transport[END_REF][START_REF] Craster | Surfactant transport on mucus films[END_REF][START_REF] Montenegro-Johnson | Physics of rheologically enhanced propulsion: Different strokes in generalized stokes[END_REF][START_REF] Sedaghat | On the effect of mucus rheology on the muco-ciliary transport[END_REF]. The main difficulty of this non-Newtonian modeling is the controversy of mucus behavior: is it shear-thinning [START_REF] Montenegro-Johnson | Physics of rheologically enhanced propulsion: Different strokes in generalized stokes[END_REF][START_REF] Craster | Surfactant transport on mucus films[END_REF] or viscoelastic [START_REF] Mitran | Metachronal wave formation in a model of pulmonary cilia[END_REF][START_REF] Sedaghat | On the effect of mucus rheology on the muco-ciliary transport[END_REF] ? This remains an open problem and we chose to assume a variable viscosity mucus, with Newtonian rheology, in the present work ; this assumption is also a good alternative to the two-layers ASL models (as explained in the modeling section).

From an algorithmic point of view various numerical methods have been presented in the literature to compute the mucus flow. The Stokeslets method [START_REF] Cortez | The method of regularized stokeslets in three dimensions: Analysis, validation, and application to helical swimming[END_REF] is one of the most popular since it suits very well for such a fluid-structure interaction problem involving a creeping flow [START_REF] Smith | Discrete cilia modelling with singularity distributions: Application to the embryonic node and the airway surface liquid[END_REF][START_REF] Smith | A viscoelastic traction layer model of muco-ciliary transport[END_REF][START_REF] Lukens | Using Lagrangian coherent structures to analyze fluid mixing by cilia[END_REF][START_REF] Montenegro-Johnson | Physics of rheologically enhanced propulsion: Different strokes in generalized stokes[END_REF]. Nevertheless the computational cost grows with the number of cilia which is prohibitive in 3D. Moreover, Stokeslets can only be used for constant viscosity fluids and other methods need to be used for more complex flows exhibiting non-linearities. A complete discretization of the fluid-structure interaction was successfully proposed [START_REF] Mitran | Metachronal wave formation in a model of pulmonary cilia[END_REF] using a coupling of the finite volume and finite element methods but the computational cost in 3D grows too quickly to envisage parametric studies. Finally the immersed boundary method [START_REF] Peskin | The immersed boundary method[END_REF] was recently used coupled to finite differences [START_REF] Lee | Muco-ciliary transport: Effect of mucus viscosity, cilia beat frequency and cilia density[END_REF][START_REF] Jayathilake | Three-dimensional numerical simulations of human pulmonary cilia in the periciliary liquid layer by the immersed boundary method[END_REF][START_REF] Jayathilake | A numerical study of muco-ciliary transport under the condition of diseased cilia[END_REF] in a parametric study context. Other works have recently coupled this method to Lattice Boltzmann techniques to naturally consider the two-layers structure of the ASL [START_REF] Li | An immersed boundary-lattice Boltzmann method for single-and multi-component fluid flows[END_REF][START_REF] Sedaghat | On the effect of mucus rheology on the muco-ciliary transport[END_REF].

The main problem to investigate numerically regarding the MCC is the large number of biological parameters in the modeling. Moreover strong variations can be observed in pathological situations such as cystic fibrosis or primary ciliary dyskinesia. Hence parametric studies are essential to understand the role of the biomechanical parameters and to estimate the influence of these pathologies on the MCC. To our knowledge, only two studies focused on this strategy [START_REF] Lee | Muco-ciliary transport: Effect of mucus viscosity, cilia beat frequency and cilia density[END_REF][START_REF] Jayathilake | A numerical study of muco-ciliary transport under the condition of diseased cilia[END_REF] using 2D computations with respectively a two-layers ASL and a constant viscosity ASL. The present work proposes an investigation of the parametric influence using both 3D computations and a complex description of the mucus. Pair interactions of parameters are also originally proposed. The results shows a good agreement with these previous works and the more complex modeling brings new insights on parametric influence under pathological situations. In particular, the variable viscosity enables to investigate of high viscosity ratio between ML and PCL as well as PCL reduction. This study is based on the use of innovative numerical methods recently developed for 3D non-homogeneous creeping flows [START_REF] Chatelin | A hybrid grid-particle method for moving bodies in 3D stokes flow with variable viscosity[END_REF][START_REF] Chatelin | Méthodes numériques pour l'écoulement de stokes 3d: fluides à viscosité variable en géométrie complexe mobile ; application aux fluides biologiques[END_REF][START_REF] Chatelin | Hybrid grid-particle methods and penalization: A Sherman-Morrison-Woodbury approach to compute 3D viscous flows using FFT[END_REF]. These methods aim at solving a non-linear coupled problem with efficient algorithms based on (i) a Lagrangian discretization of advection, (ii) an iterative projection to handle both non-homogeneity and incompressibility, (iii) the use of fast FFT-solvers. The overall computational cost is quasilinear with respect to the number of discretization points and does not depend on the number of cilia. It allows fast computations and suits very well for parametric studies.

This article brings novel insights on the influence of biomechanical parameters on the MCC, in particular for pathological situations such as cystic fibrosis. State-of-the-art numerical methods enables new investigations considering (i) a nonhomogeneous viscosity description of the ASL compatible with recent experimental advances [START_REF] Button | Periciliary Brush Promotes the Lung Health by Separating the Mucus Layer from Airway Epithelia[END_REF][START_REF] Ehre | Cystic fibrosis, An inherited disease affecting mucin-producing organs[END_REF], (ii) 3D computations of a non-linear biomechanical problem and (iii) a wide range of parameters values.

In the first part of the paper, modeling assumptions are detailed to present the governing equations and the numerical methods used for the computations. The second section presents the reference simulation and reference parameters. The two last sections detail the parametric study and the results.

Mathematical modeling and numerical discretizations

Mucus modeling

As a physiological fluid, the ASL is essentially composed of water and proteins (see figure 1 for a sketch of the flow). This fluid is often seen as a two-layers structure: it was originally proposed in [START_REF] Lucas | Principles underlying ciliary activity in the respiratory tract: Ii. a comparison of nasal clearance in man, monkey and other mammals[END_REF] and verified by later experimental works on rabbits [START_REF] Sanderson | Ciliary activity of cultured rabbit tracheal epithelium: beat pattern and metachrony[END_REF] and rats [START_REF] Yoneda | Mucous Blanket of Rat Bronchus[END_REF] trachea. Among these proteins, mucins are responsible of this viscosity variation. Mucins are polymers, released into the ASL from the goblet cells, where they hydrate and expand. Then the cross-linking of mucin polymer chains form the ML [START_REF] Ehre | Cystic fibrosis, An inherited disease affecting mucin-producing organs[END_REF] which is much more viscous. Another recent experimental work [START_REF] Button | Periciliary Brush Promotes the Lung Health by Separating the Mucus Layer from Airway Epithelia[END_REF] shows that the cilia beating helps the separation of both layers, avoiding a return of large muncins chains in the PCL. It also "demonstrate that the Gel-on-Liquid model of a two-layered airway surface is qualitatively incorrect". Hence, even if the two-layer structure of ASL was widely accepted in the literature, the existence of a clear physical PCL-ML interface remains uncertain. As both layers are essentially composed of water [START_REF] Smith | Modelling mucociliary clearance[END_REF] and as it depends on mucins maturation [START_REF] Button | Periciliary Brush Promotes the Lung Health by Separating the Mucus Layer from Airway Epithelia[END_REF][START_REF] Ehre | Cystic fibrosis, An inherited disease affecting mucin-producing organs[END_REF], the transition between PCL and ML it might be continuous.

That is why a continuous non-homogeneous viscosity approach is used in this work, instead of a bi-fluid model with constant fluid features in each layer. The viscosity is assumed to be a function of a mucins' ratio, quantifying mucins maturation, which continuously varies from 0 in the deeper PCL to 1 in the upper ML. In the following α denotes this mucins' ratio in the fluid. An ad-hoc constitutive law Φ is used to compute the viscosity η as a function of α: viscosity varies continuously from water's in the PCL to a higher viscosity in the ML (saturated in mucins). For the computations a three-parameter sigmoid function is used to initialize α, allowing to change the stiffness of the profile as well as the transition highness (see figure 2 and section 3.1). The mucins' ratio is then transported by the ASL flow: α is the solution to a convection-diffusion equation (see equation set ( 2)).

Cilia modeling

The cilia beating model used in the following numerical simulations was proposed in [START_REF] Chatelin | A hybrid grid-particle method for moving bodies in 3D stokes flow with variable viscosity[END_REF][START_REF] Chatelin | Méthodes numériques pour l'écoulement de stokes 3d: fluides à viscosité variable en géométrie complexe mobile ; application aux fluides biologiques[END_REF] and recently reused in [START_REF] Li | An immersed boundary-lattice Boltzmann method for single-and multi-component fluid flows[END_REF]. It reproduces very well experimental observations [START_REF] Sanderson | Ciliary activity of cultured rabbit tracheal epithelium: beat pattern and metachrony[END_REF] as well as state-of-the-art numerical beatings [START_REF] Fulford | Muco-ciliary transport in the lung[END_REF][START_REF] Mitran | Metachronal wave formation in a model of pulmonary cilia[END_REF]. The resulting asymmetric beating is decomposed in two parts: during "recovery stroke" cilia are bending close to the bronchial wall and beating backward (it lasts 2/3 of the beating period) whereas during "effective stroke" they are tight (cilia tips penetrate the ML) and beating forward. In this study the beating is planar and a small shift (10% of the beating period) is added between two adjacent cilia to model the methachronal synchronization. Both these parameters do not change during the following simulation.

Each cilium is modeled by a filament which has a diameter of 0.3µm. The filament centerline is a parametric curve whose position is given by solving the following 1D advection equation (which is very cheap to compute compared to the computation of 3D fluid mechanics equations):

∂L ∂t + v(t) ∂L ∂ζ = 0 ; P(0, t) = 0 ; L(0, t) = ∂ ζ P(0, t) = g(t) = (2 cos(2πt/T ), 0 , 1)
(1) where ζ is the curve parameter, P(ζ, t) is the curve position at time t, L = ∂ ζ P and v(t) = (1 + 8 cos 2 (π(t + 0.25T )/T ))/T , with T the beating period. To obtain a methachronal synchronization the shift is added to both g and v.

Once each filament position is computed at time t the characteristic function χ(t) of the volume C(t) occupied by cilia can be easily computed as well as the cilia velocity field u(t) (which is non-zero only inside cilia). Both these quantities are used in the following equations.

In this work, the fluid-structure interaction is considered to be 1-way: it means that the force applied by the fluid on the solid is not taken into account in the beating model. To our knowledge, the counter force exerted by the mucus on cilia and the associated beating adaptation was only studied experimentally [START_REF] Hill | Force generation and dynamics of individual cilia under external loading[END_REF] and the development of a coherent model for numerical computations remains an open problem.

Governing equations

Mucus is a viscous fluid which constitutes a protective blanket on the bronchial wall. Since the mucus film is very thin (10 -15µm [START_REF] Fahy | Airway mucus function and dysfunction[END_REF]), the Reynolds and Womersley numbers associated to the flow are very small. Under this hypothesis, inertia terms in the Navier-Stokes equations can be neglected and the conservation of momentum is governed by the incompressible Stokes equations.

To take into account the interaction between mucus and cilia, the Penalization technique [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] is used to force the velocity to be equal to the cilia velocity u(t) (deduced from the beating model) in C(t): the domain occupied by cilia at time t. This technique has been widely used to handle fluid-structure interactions in computational fluid dynamics. This technique is compatible with any discretization and does not require a precise knowledge of the fluid-solid interface, which saves costly computations in 3D. The principle is very simple: it consists in adding a term in the conservation of momentum equation, which is only "visible" in the immersed obstacles to guarantee the equality u = u in C(t).

Hence the equation set (2) has to be solved. Both first equations model the evolution of mucins and viscosity (described previously); both lasts govern the evolution of velocity and pressure (Stokes equations).

                     ∂α ∂t + u • ∇α -ν∆α = 0 in Ω, η = Φ(α)
in Ω,

-∇ • (2ηD(u)) + χ(t) ε (u -u) + ∇p = 0 in Ω, ∇ • u = 0 in Ω\C(t). (2)
In these equations, Ω is the computational domain, C(t) is the (closed) domain occupied by cilia at time t, so Ω \ C(t) is the fluid domain. In this coupled problem u is the velocity field, p is the pressure, D = (∇u + ∇u T )/2 the strain rate tensor and u(t) is the cilia velocity deduced from the beating model. Finally, χ(t) is the characteristic function of C(t) and ε << 1 is the Penalization parameter. In (2) the momentum equation is a quasi-static problem which implicitly depends on time through cilia velocity and positions as well as the viscosity. This Stokes problem is completed with a no slip boundary condition at the epithelium-PCL interface, a free slip boundary condition at the air-mucus interface and periodic boundary conditions in the other (proximal and transverse) directions. The advection-diffusion equation is completed with homogeneous Neumann boundary conditions at both epithelium-PCL and air-mucus interfaces and periodic boundary conditions in other directions. Moreover an initial conditions is provided for the mucins ratio α (details are given in section 3.1).

It was shown in [START_REF] Chatelin | Analysis of penalized variable viscosity 3d stokes equations coupled to diffusion and transport[END_REF] that the problem (2) admits a unique solution in adequate Sobolev spaces. Moreover, when the Penalization parameter tends toward zero, this solution weakly converges toward the solution of the physical problem: u = u in the solid, (u, p) is the solution to the incompressible Stokes problem in the fluid with convection-diffusion of α.

Numerical algorithms

To compute efficiently the non-linear convection-diffusion equation (first equation of ( 2)) a splitting strategy is applied to uncouple both physical phenomena. This technique enables to choose a suitable discretization for each phenomenon: a Lagrangian method to compute convection and a grid-based method to compute diffusion. Error estimates for the splitting method were introduced in [START_REF] Beale | Rates of convergence for viscous splitting of the navier-stokes equations[END_REF] and have been widely used in the literature. In this work the second order splitting introduced in [START_REF] Chatelin | Mucus and ciliated cells of human lung : Splitting strategies for particle methods and 3d stokes flows[END_REF] is used: it permits to solve a diffusion problem after each advection computation to take advantage of the diffusion regularizing effects at each sub-step. The Lagrangian discretization of advection has two interested features: (i) it has a linear computational cost with respect to the number of particles and (ii) it is not subject to a CFL stability condition (even using explicit time discretizations), so large time steps can be achieved avoiding unnecessary computations of the velocity u. In this context of mucus flow, advection of mucins dominates diffusion so the parameter ν is small enough to permit explicit finite difference computations of the heat equation. When the diffusion CFL is too restrictive (meaning time steps are too large), FFT-based solvers are chosen [START_REF] Swarztrauber | Efficient FORTRAN subprograms for the solution of elliptic partial differential equations[END_REF] to get efficient implicit computations (see below).

The flow velocity u is computed by solving the two last coupled equations of (2) on a regular Cartesian grid, which correspond to non-homogeneous viscosity Stokes equations. An iterative projection method was introduced in [START_REF] Chatelin | A hybrid grid-particle method for moving bodies in 3D stokes flow with variable viscosity[END_REF] to compute this problem. It ensures an accurate computation of (i) boundary conditions -an inherent problem of projection methods [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF], (ii) immersed solid velocity, (iii) invicid velocity and (iv) nonhomogeneous effects. At each step of this iterative algorithm, a penalized Poisson problem is solved using a precise and efficient Sherman-Morrison-Woodbury approach [START_REF] Chatelin | Hybrid grid-particle methods and penalization: A Sherman-Morrison-Woodbury approach to compute 3D viscous flows using FFT[END_REF]. Its solution is then corrected by a projector which is also the solution to a Poisson problem. All these Poisson problems are computed with a FTT-based solver [START_REF] Swarztrauber | Efficient FORTRAN subprograms for the solution of elliptic partial differential equations[END_REF] which ensures a quasi linear computational cost with respect to the number of discretized points: O(N log N). This is made possible by the a priori structure of the finite difference operator on the regular Cartesian grid.

Using the Lagrangian method to compute transport phenomena and fast-solver-based algorithms for the Stokes problem, the overall computational cost is quasi-linear with respect to the number of discretized points: O(N log N). Moreover this computational cost was demonstrated [START_REF] Chatelin | A hybrid grid-particle method for moving bodies in 3D stokes flow with variable viscosity[END_REF] to be identical using homogeneous or non-homogeneous viscosity and it is independent to the number of immersed cilia. This result is very interesting for 3D computations since the computational time is significantly reduced compared to matrix assembling based methods such as straightforward finite differences/elements/volumes methods. These numerical methods are particularly suitable for parametric studies which requires a large number of simulations.

Numerical simulations of mucociliary clearance

Viscosity profiles

To model the transition between both ASL layers a sigmoid profile is used -see figure 2. Its expression depends on three parameters: the viscosity ratio β between ML and PCL, the transition length δ and the transition stiffness γ. The following expressions are used to initialize the mucins' ratio (equation ( 3)) and to compute the viscosity as a function of α at any time t (equation ( 4)) .

α(x, y, z, t = 0) = arctan (γ(z -δ)) -arctan(-γδ) arctan (γ(H z -δ)) -arctan(-γδ) (3) η(x, y, z, t) = Φ(α(x, y, z, t)) = η PCL (1 + βα(x, y, z, t)) (4)
where β = η ML /η PCL is the viscosity ratio between ML and PCL.

Reference parameters and simulation

The reference parameters used for the simulations are gathered in the table 1 (including viscosity parameters introduced previously). On the boundaries a no slip condition is imposed on the bronchial wall and a free slip condition is imposed on the air-mucus interface (assuming a flat interface). Finally periodic boundary conditions are imposed in the x and y directions.

On figure 3 several snapshots of a simulation are presented. This simulation was performed using the reference parameters of table 1. It presents two isosurfaces of mucins' ratio (which is transported by the flow), as well as velocity magnitude on the computational box boundaries.

During the simulation the mean velocity of the ML in the proximal direction (x) is computed. It is then averaged over three beating cycles (in fact this number of cycles do not change significantly the results):

U = 1 3T H x H y (H z -δ) 3T 0 H z δ H y 0 H x 0 u x (x, y, z, t)dxdydzdt (5)
This scalar quantity is a very good indicator of the MCC efficiency and it is used in the following parametric study to compare the different parameter sets. In the following we refer to this ML velocity as "MCC efficiency". In particular the power that cilia exert to achieve the beating is not investigated. It would be an interesting investigation in future works, particularly when the counter force exerted by the fluid will be integrated in the modeling.

Parametric study: methodology and results

Methodology for the parametric study

For the parametric study one simulation parameter is selected and K simulations are computed using K different values of this parameter. The other parameters remain constant, equal to the value presented in the table 1. The same methodology is then used for two parameters. For each simulation the quantity U introduced in equation ( 5) is computed to quantify the MCC efficiency. Hence for each simulation only a scalar quantity is recovered to avoid unnecessary memory storage.

The computational efficiency of the numerical algorithms enable sequential computations for each simulation using a 128×128×64 refinement (so more than a million grid points are used). The parametric study is achieved using HPC resources by performing simultaneously a large number of independent sequential simulations.

Variation of the beating frequency

When the beating frequency varies, results present a linear variation of the mean mucus velocity (that is why the figure is not presented). This result, also reported in [START_REF] Lee | Muco-ciliary transport: Effect of mucus viscosity, cilia beat frequency and cilia density[END_REF], is a direct consequence of the time-dilatation of the quasi-static Stokes equations.

Variation of the cilia length

In this section the cilia length varies from 3µm to 12µm using an increment of 1µm (ASL height is 13µm and the transition length between PCL and ML remains constant: 8µm). The results are presented on the figure 4. The MCC is the most efficient when the cilia length is larger than the ASL thickness, in fact it is maximized when cilia are slightly longer (9µm) than this thickness. It means that cilia need to penetrate into the ML to achieve an efficient MCC. In fact the best configuration is achieved when cilia penetrate the ML only during the effective stroke.

Mucociliary transport is divided by 10 when cilia are the smallest. In [START_REF] Jayathilake | A numerical study of muco-ciliary transport under the condition of diseased cilia[END_REF] similar results were reported (in particular the velocity magnitude is similar: ∼ 10 -5 m.s -1 but smaller than this work, due to 3D effects), even if the decreasing at high cilia length is lower in the present computations.

Variations of the mucus layer height

Figure 5 presents the influence of ML height variations. It shows that the mucociliary transport is reduced when the ML is Comparison between linear and sigmoid profiles using parameters of table 1. too thin. These computations are performed with a constant PCL thickness. The figure also shows a reduction for very high ML thickness. A 30% maximum difference is observed on this figure, hence it appears that this parameter have a significant impact on the MCC. Nevertheless, even if the mean velocity slowly decreases when the thickness grows, the throughput keeps growing. This result tends to prove that MCC can handle mucus overloads (if and only if other parameters do not change). On the contrary, it shows that mucus throughput is highly reduced when the film is too thin which can cause dangerous contaminations of unprotected bronchus.

Viscosity variations: linear profile

In this section a linear viscosity profile is used for α. Instead of equation ( 3) the following expression is used to initialize mucins' ratio:

α(x, y, z, t = 0) = z/H z (6)
The relation ( 4) is still used to define the viscosity as a function of the mucins' ratio and the coefficient β varies in this section from 10 -2 (almost constant viscosity ASL) to 10 3 (mucus is 1000 times more viscous than PCL) using 17 increments. The main interest is to get only one parameter involved, which reduces the cost of the parametric study but simplifies the modeling. The figure 6 presents these results. When β is small, the MCC remains constant but it decreases as β grows. In fact the viscous dissipation is more important when cilia are in the upper part of the ASL (during effective stroke), hence the fluid velocity is less important and is divided by 1.7 assuming a ML 1000 times more viscous than PCL. This figure suggest that this parameter has an important impact on the MCC and justifies the lack of mucociliary transport under cystic fibrosis situations.

Viscosity variations: sigmoid profile

As three parameters are used to describe the sigmoid profile, three figures are shown to present pair interactions between these parameters. As previously the same quantity (U) is displayed to quantify the MCC efficiency. Figures 7,[START_REF] Chatelin | Hybrid grid-particle methods and penalization: A Sherman-Morrison-Woodbury approach to compute 3D viscous flows using FFT[END_REF]

βγ

In this section 17 values of β are tested (as in previous section) against 14 values of γ. It means that 238 simulations are performed to obtain the figure 7. This figure shows that the MCC can be divided by 2.55 between the extrema. The most favorable case is obtain when β is moderate [START_REF] Cortez | The method of regularized stokeslets in three dimensions: Analysis, validation, and application to helical swimming[END_REF][START_REF] Craster | Surfactant transport on mucus films[END_REF][START_REF] Ehre | Cystic fibrosis, An inherited disease affecting mucin-producing organs[END_REF][START_REF] Fahy | Airway mucus function and dysfunction[END_REF][START_REF] Fulford | Muco-ciliary transport in the lung[END_REF][START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF][START_REF] Hall-Stoodley | Bacterial biofilms: from the natural environment to infectious diseases[END_REF][START_REF] Hill | Force generation and dynamics of individual cilia under external loading[END_REF][START_REF] Jayathilake | Three-dimensional numerical simulations of human pulmonary cilia in the periciliary liquid layer by the immersed boundary method[END_REF][START_REF] Jayathilake | A numerical study of muco-ciliary transport under the condition of diseased cilia[END_REF][START_REF] Lee | Muco-ciliary transport: Effect of mucus viscosity, cilia beat frequency and cilia density[END_REF] and γ is large enough, meaning a stiff transition between ML and PCL.

On the opposite to the result obtained for the linear profile (see figure 6 for a comparison), one can notice that, using a sigmoid profile, the most efficient MCC is not obtained for the smallest β, but for a viscosity ratio β = 10. It means that the stiffness of the transition has a major role: it permit to handle the accumulation of cross-linked mucins in the ML without reducing the MCC efficiency. This result demonstrate that a (moderated) more viscous ML increases the MCC efficiency.

Finally, for a very large ML viscosity MCC is highly deteriorated even using a stiff transition. This deterioration is even more important with the sigmoid profile than with the linear one.

βδ

17 values are still used for the parameter β whereas 10 values are used for the normalized transition length δ ; so 170 computations are necessary to obtain the figure 8. The most efficient MCC is obtained for similar values of β (around 10) as for the previous figure and for a normalized transition length between 0.6 and 0.9. This result confirms the assumption formulated in [START_REF] Fulford | Muco-ciliary transport in the lung[END_REF]: cilia need to penetrate the ML during the effective stroke to maximize MCC. On this figure a similar factor is obtained between extrema. The minimal values are obtained for very viscous mucus (large β) as for previous computations.

γδ

On figure 9 the influence of both γ and δ is presented for respectively 14 and 10 values (as in previous computations). The most efficient MCC is still obtained for the same values as previously, meaning largest γ and δ ∼ 0.6 (meaning a penetration of cilia in the ML during the recovery stroke). Here a factor 7.5 is observed between maximum and minimum result. This minimum is obtained for small δ and large γ. This figure exhibits several zone of vertical isolines, which means that only one parameter (δ) has an influence on the MCC in this zone and the number of computations could be reduced. As strong variations are observed on this figure both parameters can be identified as important, but not in all the regions: where isolines are vertical only δ has a significant impact.

Discussion and future works

In the previous section the influence of several biological parameters on the MCC was investigated. It was shown that MCC depends linearly on the beating frequency, hence it is a good non-dimensional parameter for the system. Three lengths were identified to play a major role on the MCC: the cilia length, the ASL height and the transition length between PCL and ML. Finally, the viscosity of ML and the stiffness transition between PCL and ML are also relevant parameters. In particular, this study tends to show that for healthy conditions a stiff transition between PCL and ML tends to maximize MCC efficiency, but this becomes wrong when PCL height is decreased as well as for highly viscous mucus.

The parametric study showed that under pathological conditions corresponding to cystic fibrosis (increase of the ML viscosity and reduction of the PCL height) the MCC is dramatically decreased, which leads to mucus stagnation and causes secondary infections: pathogens proliferation in the mucus can seriously contaminate the lung. Nevertheless, results obtained for the ASL height tend to prove that cilia beating can handle mucus overload. A cross study with this parameter and ML viscosity could also be investigated in the future.

A limitation of the present model (which remains, to our knowledge, an open question) is precisely these secondary infections: they would have a feedback on the biological parameters, which is not taken into account in the modeling. For example a cystic fibrosis patient, who has an increased ML viscosity, will present pathogen proliferation, for example Pseudomonas aeruginosa or Burkholderia cepacia [START_REF] Hall-Stoodley | Bacterial biofilms: from the natural environment to infectious diseases[END_REF], in the ASL, which will change again the ASL viscosity and so on... This requires a more realistic modeling which would increase dramatically the computational cost, and would currently make impossible the parametric study.

In the future, Non-Newtonian effects will also be integrated in the modeling. Since there is a controversy on the constitutive model, additional measurements are required to choose an appropriated rheology (shear-thinning or viscoelastic). An α-dependency of this parameters could be considered to model a Newtonian PCL and a non-Newtonian ML. From a computational point of view this will not be a locking point: computations of the resulting non-linearity will be integrated in the projection iterations. Nevertheless, at least two parameters will be added in the model which will increase the complexity of the parametric study.

Four meaningful parameters have been identified in this parametric study: two characteristic lengths (assuming two ratios involving the three previous length parameters), the viscosity of ML and the stiffness transition between ML and PCL. This would become an harder problem assuming non-Newtonian mucus (adding at least two parameters) and other biological parameters not investigated in this work, for example cilia density and cilia beating patterns [START_REF] Jayathilake | Three-dimensional numerical simulations of human pulmonary cilia in the periciliary liquid layer by the immersed boundary method[END_REF][START_REF] Jayathilake | A numerical study of muco-ciliary transport under the condition of diseased cilia[END_REF]. This more complex study would be performed using global sensitivity analysis or other optimization techniques to select appropriated parameter sets and avoid unnecessary computations in small variation zones of the parameters' space. These techniques will substantially reduce the number of simulations.

Figure 1 :Figure 2 :

 12 Figure 1: Sketch of the mucus flow on the bronchial wall: the nonhomogeneous mucus is propelled by beating ciliated epithelium cells.

Figure 3 :Figure 4 :

 34 Figure 3: Snapshot of the simulation performed with the parameters of the table 1. In red and black isosurfaces of level 0.1 and 0.9 respectively, of mucins' ratio are represented. The colormap present the velocity magnitude on the domain boundary. Top left: initial state at t = 0s. Top right: snapshot at t = 0.14s. Bottom left : snapshot at t = 0.273s. Bottom right: snapshot at t = 0.482s.

Figure 5 :

 5 Figure 5: Influence of the ML height on the MCC (PCL height remains constant for each computation).

Figure 6 :

 6 Figure 6: Influence of the viscosity ratio β between ML and PCL viscosity.Comparison between linear and sigmoid profiles using parameters of table 1.

  and 9 present the results for βγ, βδ and δγ interactions respectively.Mean mucus velocity (m/s)

Figure 7 :

 7 Figure 7: Influence of the viscosity parameters β and γ on the mucociliary clearance.Mean mucus velocity (m/s)

Figure 8 :

 8 Figure 8: Influence of the viscosity parameters β and δ on the mucociliary clearance.Mean mucus velocity (m/s)

Figure 9 :

 9 Figure 9: Influence of the viscosity parameters β and γ on the mucociliary clearance.

Table 1 :

 1 Reference parameters for the numerical simulations (gathered from[START_REF] Fahy | Airway mucus function and dysfunction[END_REF][START_REF] Sanderson | Ciliary activity of cultured rabbit tracheal epithelium: beat pattern and metachrony[END_REF][START_REF] Smith | Modelling mucociliary clearance[END_REF])

	Name	Value
	Cilia length	8 µm
	Cilia diameter	0.3 µm
	Number of Cilia	25
	PCL viscosity (η PCL )	0.001 Pa.s
	Viscosity ratio (β = η ML /η PCL )	20
	Viscosity transition stiffness (γ)	10
	Transition length (δ)	7 µm
	ASL height (H z )	13 µm
	Beating frequency	10 Hz

Acknowledgments

This work was granted access to the HPC resources of P2CHPD (Pôle de Calcul Hautes Performances Dédiés, from Fédération Lyonnaise de Simulation et de Modélisation Numériques). This work was supported by the ANR Grant BioFiReaDy, under the contract number ANR-2010-JCJC-0113-01.

Conflicts of interest declaration

We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

We confirm that the manuscript has been read and approved by all named authors and that there are no other persons who satisfied the criteria for authorship but are not listed. We further confirm that the order of authors listed in the manuscript has been approved by all of us.

We confirm that we have given due consideration to the protection of intellectual property associated with this work and that there are no impediments to publication, including the timing of publication, with respect to intellectual property. In so doing we confirm that we have followed the regulations of our institutions concerning intellectual property.

We understand that the Corresponding Author is the sole contact for the Editorial process (including Editorial Manager and direct communications with the office). He is responsible for communicating with the other authors about progress, submissions of revisions and final approval of proofs. We confirm that we have provided a current, correct email address which is accessible by the Corresponding Author and which has been configured to accept email from the Journal of Biomechanics.

R. Chatelin and P. Poncet