Effect of bunching of cilia and their interplay on muco-ciliary transport.

Cilia are hair-like organelles projecting from a eukaryotic cell, used either for locomotion or as sensors. Cilia commonly occur in patches. To take this into consideration, we represent cilia in multiple patches, instead of the conventional 'dense mat' representation. We focus on the combined action and interplay of these patches. The effects of varying the frequency, spacing and phase lag of the beating of one cilia bunch with respect to the beating of adjacent patches are studied. We model the Airway Surface Liquid (ASL) as a three-layer structure. The possibility of an optimum frequency of beating is noted and the change of mucous flow under different spacing and phase differences are observed.

Introduction

The study of cilia has shown escalating physiological significance in the past decade owing to the discovery of its role in the proper functioning of many mechanisms in the human body and other organisms. In a relatively short period of time, over ten diseases, such as primary ciliary dyskinesia (PCD), hydrocephalus, polycystic liver and kidney disease, some forms of retinal degeneration and Joubert Syndrome, have been linked to the malfunctioning of cilia. They have been classified as ‘ciliopathies’ (Badano [1]), (Adams et al. [2]). Some of the most vital ciliopathies are caused due to the lack of coordination in ciliary movement in the epithelium. Cilia exhibit a distinctive, complex asymmetric manner of movement. Each beat cycle comprises of an ‘effective stroke’ and a ‘recovery stroke’. In the former, the cilia are relatively straight, and their effect on the fluid is maximized. In the latter, their effect on the surrounding medium is minimized since the cilium has a more curved shape. Any disruption of this characteristic pattern can lead to physiological problems such as when the cilia on the epithelium show hyper/hypo motility or disoriented arrangement and beat direction (Van's Gravesande et al. [3]). The PCD is a classic example of such a condition. It involves the cilia on the respiratory tract and reproductive organs. PCD leads to poor muco-ciliary clearance and hence infection (Rayner et al. [4]), (Bush et al. [5]). This disease manifests itself during the embryologic phase of development and is a result of inefficient and unsynchronized ciliary movement.

Normally-functioning cilia determine the position of the internal organs during early embryological development. Therefore, individuals with PCD have a 50% chance of developing situs-inversus (El Zein et al. [6]) (major visceral organs are reversed or mirrored from their normal positions), accompanied by progressive sinusitis, and bronchiectasis. A sum total of these effects is called the ‘Kartergener Syndrome’ (Dhar et al. [7]).

Another severe condition associated with dysfunctional cilia is the Joubert Syndrome (Doherty et al. [8], Brancati et al. [9]) which is a rare genetic disorder that affects the cerebellum. This affects the balancing and coordinating control of the brain on the skeletal muscles. It is caused due to dysfunctional molecular structure of the cilia and adversely affects the numerous critical developmental signaling pathways essential to cellular development.

Also, Usher syndrome (USH), is a heterogeneous disease, which, in rare cases, can be associated with bronchiectasis, chronic sinusitis and reduced nasal mucociliary clearance, which is indicative of ciliary dyskinesia (Bonneau et al. [10]).

The diseases mentioned hitherto are genetically acquired. Several ‘acquired diseases’ are also consequent upon malfunction of cilia or may have ciliary dysfunction as a symptom. For example, when infected with various common cold viruses, the respiratory epithelium reacts by shedding individual cilia as well as entire ciliated cells (Rautiainen et al. [11]). Cilia may also become withdrawn into the cell body, as is seen in the nasal epithelium during coronavirus infection (Afzelius [12]). A substantial loss of cilia is likely to result in rhinorrhoea (a condition where the nasal cavity is filled with a significant amount of mucous fluid). A reduction to about half the normal number of ciliated cells has been recorded from the tracheae of hamsters that were deprived of vitamin A since birth (McDowell et al. [13]).

The cilia in the ependyma (the thin epithelial membrane lining the ventricular system of the brain and the spinal cord) of the brain, though not exposed to the same harsh environment as the tracheal cilia, can be affected by vitamin A deficiency (McDowell et al. [13]) in the same way as the respiratory cilia. They can also be infected by viruses such as the mumps virus or various types of influenza virus (Johnson [14], Tanaka et al. [15]). The first step of the infection is replication of the virus in close vicinity to the cilia.

This leads to a loss of most cilia. As a consequence hereof, there are no cilia to prevent particular matter from settling down on the ependymal surface (Stokroos et al. [16]).

Study of the physiology of numerous animals has shown that the synchronization of cilia occurs in patches. Patches are distinct regions of ciliary activity each separated from one-another, and spread over ciliary epithelium lining various organs of the body. These patches of cilia span some few micrometers (Eshel and Priel [17]). It appears that in experimental studies, coordination of cilia has not been seen to exceed distances greater than a few micrometers.

In some cases, the patches of cilia are physically separated by clefts or ridges. In Sanderson and Sleigh [18], patches of cilia seen on rabbit tracheal epithelium by SEM are separated by such clefts. The patches of ciliated surfaces are pressed together in these cleft regions. It is mentioned that the areas of activity and the direction of metachronic waves traveling across the areas are not randomly changing. Examination of high-speed cine-film showed that the activity on the epithelial surface is comprised of many small patches of regularly beating cilia. The extent of these patches depends on the type of control over ciliary beating. In frog palate, where patches are highly localized, it is under neural control. In rabbit epithelium, where the action of cilia is spontaneous, larger areas can operate in synchronization.

In other cases the patches of cilia are separated by groups of other cell types. In a publication by Blake [19] which describes various mechanisms of mucus clearance in the mammalian lung, it is stated that cilia certainly do not occur in a dense mat form as was previously thought and while ciliated cells cover approximately 30–60% of the epithelium, goblet cells occupy approximately 20–40% of the epithelium.

Also in a study by Antunes et al. [20] using ‘Scanning Electron Microscopy (SEM)’, the degree of ciliation of the murine respiratory system was observed. It was found that the ciliated cells occur only in scattered patches. However, the degree of ciliation varied from area to area; for example, while the mouse tracheal epithelium had 30–35% of the surface covered in ciliated epithelium, the mouse nasal septum had approximately 90% of the epithelial surface covered with cilia.

Thus, whether groups of ciliated cells are separated by clefts or simply by other cell types, the overall ciliary function is in patches.

Studies which assume the cilia to exist in a dense mat form are mostly mathematical. In such sources, attention is paid to the beating characteristics and function of individual cilia. In order to extend the effect of individual cilia over some area, the treatment of individual cilia is simply extended over some space. Indeed cilia are represented either as slender bodies or with a spatially continuous force distribution (Phenomenological modeling) (Smith et al. [21]). While the discrete cilia method is more mathematically precise, volume force models are generally used to include physical effects that are difficult to take into account in a more precise mathematical model even though such models require a number of simplifying approximations in their formulation.

In neither of the above approaches does one appreciate the reality of the situation i.e. that cilia in fact occur in patches. Therefore, we realize the need for one more factor to be analyzed, i.e. how the cilia, in different locations of the epithelium assist each-other in muco-ciliary clearance.

In an analytical work by Blake [19], cilia are considered to exist in an infinite array. The mucous flow rate is evaluated as a function of ciliary beat frequency (CBF), cilium length and cilia concentration. Though this provides a comprehensive study of how mucous flow is affected by important parameters of ciliary action, it does not fully capture the reality of cilia, i.e. that cilia are unevenly distributed over the epithelium and hence another parameter to be studied is how these tufts of cilia interact and cooperate with each-other. Similar models have been considered by Smith et al. [21] and Keller [22].

In this study, we try to understand this hitherto ignored aspect of ciliary action i.e. the collective behavior of the cilia patches. We do so by investigating the effects of varying the frequency, spacing and phase lag of the beating of one cilia bunch with respect to the beating of adjacent patches. Such an integrated study of all these factors is previously unseen.

Before proceeding with the analysis of our results, it is important to realize the implications of our studies. As mentioned earlier, an understanding of the biological mechanisms underlying ciliary motion plays a vital role in diagnosis and treatment of many fatal diseases. Extending the scope further, we find that the study of the synchronized movement of cilia has inspired to mimic them and develop ‘artificial cilia’. These artificial cilia, apart from application in MEMS and micro-robotics, have been used on lab-on-a-chip devices for varied applications including microfluidic pumps, acoustic detection, and heat transfer (Khaderi [23], Breidenich [24]).

Further, ‘grafting’ of ciliary epithelium has been shown to be possible. ‘Grafting’ of cilia into organ systems of those with vast areas of impaired ciliary function could be life-saving. In a study by Bootz, and Reuter [25], the aim was to establish whether it was possible to carry out such grafting in an animal model. They used purebred strains of Lewis blood group rats to avoid host-versus-graft rejections of mucous membrane transferred from a donor to a recipient animal. Their findings showed that respiratory epithelia heal when grafted onto a recipient site which has an adequate blood supply but do not lose their initial differentiation or functions. Such results are important clinically for reconstruction of the respiratory tracts.

Also, in experimental studies it has been shown that the frequency of ciliary beating plays a key role in altering the rate of transport of the fluid medium. Alteration of levels of ATP and ACH (acetylcholine) in the physiological environment has shown to drastically enhance ciliary beating, leading to an almost linear increase in ciliary tip forces (Zvi Teff [26]). In regard to inter-ciliary spacing, it has been studied by Grϋnbaum [27] in a computational hydrodynamic model of the ciliated tentacle arrays (using M.edulis as a case study), that spacing between cilia can alter the efficiency of clearance of the fluid medium in which cilia occur. Therefore a good understanding of the parameters looked at in our study can gain us valuable insight about the efficiency of muco-ciliary clearance.

Our approach is computational and uses a hydrodynamic model that calculates mucus flow rates induced by ciliary movement. We model a generic ciliary epithelium and this model could be used to represent ciliated epithelia in various organs and species.

At the beginning, we validate our computational model using published experimental and analytical studies. Our model of the ciliary epithelium is based on the three layered model of the ASL. Having proved that our model corresponds to reality, we then extend our model to incorporate the ciliary patches, since in real life cilia do not occur as a dense mat, but rather in patches. Finally, we discuss the observations made from our simulations. We show that the mucous clearance can indeed be greatly enhanced in certain conditions of inter-patch coordination. These observations will help us understand the ciliopathies and their possible management, in a whole new light.

Development of the model and validation

Physical model

The model is developed using the three layer model introduced in Smith et al. [28]. It considers a three-layered structure of the Airway Surface Liquid (ASL), including a bottom-most layer of Newtonian Periciliary fluid (PCL), an overlying traction layer and a top-most thick layer of mucous ( Fig. 1). The traction layer is essentially an extension of the linearly viscoelastic mucous layer, which defines the maximum extent to which the cilia-tips extend and hence the traction layer can be assumed to be responsible for transmission of cilia-tip forces to the mucous layer. The cilia length is considered to be around 6 μm, and most of this length is submerged in the PCL except during the momentary extension into the traction layer during the effective stroke. Thus, the PCL is taken to be 5.4 μm, and the traction layer about 0.6 μm. The overlying mucous layer is taken as 4 μm. The width of the domain was taken to be considerably larger than the height, in order to compensate for edge effects. A two dimensional model is created, with the effect of the cilia represented as propulsive body forces in the traction layer and resistive viscous forces in the PCL. A simple Maxwell model of viscoelasticity is used and is established in the traction and mucous layers.

Fig. 1

The three-layered structure of the Airway Surface Liquid (ASL), with a bottom-most layer of Newtonian Periciliary fluid (PCL), an overlying traction layer and a top-most thick layer of mucous. The PCL is taken to be 5.4 μm, and the traction layer about 0.6 μm. The overlying mucous layer is taken as 4 μm.

Mathematical model

Governing equations

The fundamental equations solved in all the three domains were the Navier–Stokes equations in two dimensions, along with the supplementary continuity equation. They are presented here for reference.

Continuity equation:

Momentum equations:

The propulsive force of the cilia is due to the beat cycle asymmetry. When the cilia are in their effective stroke, they extend into the traction layer, whisking forward the overlying mucous layer just like a sheet of fabric being propelled forward by numerous swaying hands holding it up from the bottom. On the other hand, during the relaxation stroke, the cilia lower themselves, to be completely submerged in the PCL, as if they were ducking their heads. Consequently, the mucous is not pulled backward during the relaxation stroke.

The dense mat of cilia resists the flow of PCL, analogous to a porous medium, where fluid flow is obstructed due to a very large pressure drop across the medium. Also, when the cilia oscillate, the fluid close to the surface of a cilium will move with similar velocity as that of the cilium. The net effect of these two factors is that the cilia reduce flow of the PCL, though they cause significant oscillations in the PCL. The resistive force, however, acts only in the PCL since the cilia hardly penetrate either of the other two layers.

Periciliary fluid resistive force

In the horizontal direction:

In the vertical direction:where and are the PCL resistive forces in the horizontal and vertical directions respectively and are the resistance coefficients of the PCL in the x (horizontal) and y (vertical) axes respectively; and are the cilia velocities in the x (horizontal) and y (vertical) axes respectively; and are the resulting PCL velocity at a given value of (x,y,t) coordinates and is the height of the PCL layer.

(Considering h to be the height between the base of the PCL to the base of the traction layer. The distance between h and L represents the total height of the traction layer and between L and H is the total height of the mucous layer), Fig. 1.

Ciliary velocities

The expressions for the ciliary velocities were formulated to incorporate the periodicity of the beat cycle, the linear increase of velocity with distance from the epithelium and the short duration of effective stroke (about 20% of the beat cycle). Also, the integral of ciliary velocity over time needs to be zero, since in each beat cycle, each cilium generates an ellipse with its tip (Ross and Corrsin [29]). Thus the most accurate formulation for ciliary velocity is given below, after incorporating the above mentioned specifications.

The expressions for the ciliary velocity used are

In the horizontal direction:

In the vertical direction:

Here,where is the wavelength (about 30 μm), the cilia beat frequency (60 rad/s), the fraction of the cilia beat stroke that is occupied by effective stroke.

Ciliary propulsive forces (in traction layer)

The propulsive forces are also formulated by a Fourier expansion with 15 terms. The expression for these forces indicates that they decrease monotonically to zero from the base to the top of the traction layer.

The expressions for ciliary propulsive forces are:

In the horizontal direction:

In the vertical direction: γ is a coefficient of resistance which depends on the geometry of the cilia and the gap between each cilia. There are two separate values of the coefficient of resistance which correspond to the horizontal and vertical directions. The values of the coefficients of resistance also depend on the dynamic viscosity of the fluid for which they are calculated. For both x and y directions, the γ values are of the order of 1012 in SI units. The different γ values for the PCL and of the traction layer are denoted by using the superscripts ‘P’ and ‘M1’ for the PCL and the traction layer respectively. All the force terms mentioned hitherto are in the units of N/m3.

The above form of input conditions and Fourier series coefficients were obtained from Smith et al. [28], with reference to which we validated our model. The expressions correspond exactly with that prescribed in this reference paper. In order to prove this point, we present the plots of the cilia velocities and forces generated on MATLAB and compare them with the plots in the reference paper.

In Fig. 2(a) we show the plot for the horizontal cilium velocity along a short length across the domain whereas Fig. 2(b) represents the vertical cilium velocity plot across the same distance along the domain. The red curves on each plot shows the curve as obtained from the reference paper (Smith et al. [28]) whereas the blue curve is the plot generated by MATLAB and represents the input to our model. The plot refers to the expressions given above for cilium velocities. For the plots, we have fixed the value of y (vertical coordinate) at 4 μm, which is in the interior of the PCL and we have fixed the value of t (time) at 4 s. Thus, the ciliary velocities have been plotted with respect to the horizontal coordinates.

Fig. 2

(Cilia velocities) The first plot (a) represents the horizontal cilium velocity whereas (b) represents the vertical cilium velocity. The solid line curve shows the curve as obtained from the reference paper [28] whereas the dashed line curve is the plot generated by MATLAB and represents the input to our model.

Fig. 3 shows the three-layered structure of the ciliary epithelium with an aim to clearly explain the position of the individual layers with respect to each-other. One can see the PCL, traction and mucous layers. We have shown the images at 4 s, 4.2s, 4.5 s, 4.7 s and 5 s of solution respectively. We find that the velocity field is positive at all points within the domain, indicating unidirectional transport of mucous and there is vortical motion in the PCL layer, indicative of ciliary action. We also present a quiver plot of the velocity field at 5 s in order to show clearly the circular motion in the PCL and the relatively uniform velocity in the mucous.

Fig. 3

The figures show the three-layered structure of the ciliary epithelium. Visible are the PCL, traction and mucous layers .We also show the applied boundary conditions. These features are shown clearly in the diagram, showing the enlarged view of our domain (a). In (b) we show the mesh used. Subsequently, snapshots were taken at 4 s, 4.2 s, 4.5 s, 4.7 s and 5 s of solution respectively from Figs. 3(c)–g. (h) shows the domain at 4 s of solution with an arrow plot to represent the velocity field. (b) The extra-fine mesh of the model, comprising 2444 triangular elements and 11,226 degrees of freedom. (c) A snapshot of the simulation model at 4 s solution time showing the three-layered structure of the ciliary epithelium. (d) A snapshot of the simulation model at 4.2 s solution time showing the three-layered structure of the ciliary epithelium. (e) A snapshot of the simulation model at 4.5 s solution time showing the three-layered structure of the ciliary epithelium. (f) A snapshot of the simulation model at 4.7 s solution time showing the three-layered structure of the ciliary epithelium. (g) A snapshot of the simulation model at 5 s solution time showing the three-layered structure of the ciliary epithelium. (h) The domain at 4 s of solution with the velocity field represented by an arrow field.

The quiver plot shows the uniform and steady mucous velocity through the cycle. This is how the mucous flows except for the time of penetration of the cilia into the traction layer. It also indicates that there is circulatory motion in the PCL. This circulatory motion caused by the interaction of the cilia motion and the pressure gradient. In fact, these circulatory motions probably assist the transfer of particles, such as pathogens from the PCL to the mucous layer for efficient removal (Matsui et al. [30]).

One of the most distinctive features of the mucous layer is its viscoelasticity. Apart from the characteristics of shear thinning, spinnability, adhesiveness etc., there are certain unique behaviors of the mucous which can be best explained by its viscoelastic nature. It is observed for example, that during the brief period of the beginning of the recovery stroke, the mucous reacts elastically (with minimum relative slip). Then during the recovery stroke, a period of sufficient length must again be granted for the cilium to release itself from the contact of the mucous (Silberberg [31]). This is because of the elasticity of the mucous. Thus, unlike most fluids which are characterized only by their viscosities, mucous exhibits complex visco-elastic behavior. It can be argued that non-Newtonian fluids can be exploited by biological systems to their advantage because viscoelastic stresses allow a tuning of kinematics of transport and locomotion in a manner that is not possible from Newtonian fluids (Lauga [32]). In the same reference, it is observed that the energetics of a fluid which has some elasticity is more favorable in many biological systems such as that of swimming spermatozoa.

Hence there is need for the mucous to be modeled with an elastic component, in addition to the fluid viscosity, which would allow it to deform and then recoil in response to penetration. A simple model for this purpose is the Maxwell Model. The Maxwell model is defined by the parameters of ‘single relaxation time’, λ and ‘steady flow viscosity’, μ.

The Maxwell model can be represented by a purely viscous damper and a purely elastic spring connected in series. According to this model, if the visco-elastic fluid is put under a constant strain, the stresses gradually relax. When a material is put under a constant stress, the strain has two components. First, an elastic component occurs instantaneously, corresponding to the spring, and relaxes immediately upon release of the stress. The second is a viscous component that grows with time as long as the stress is applied. The Maxwell model predicts that stress decays exponentially with time. Thus the spring should be visualized as representing the elastic or energetic component of the response, while the dashpot represents the conformational or entropic component.

If is the stress tensor for visco-elastic effect, is the strain tensor for visco-elastic effect and is the strain rate tensor for visco-elastic effect, the constitutive relations for the spring and dashpot are given as:and respectively, where k is spring constant and μ is steady flow viscosity (if all the time derivatives of a flow field are assumed to vanish, the flow is considered to be a steady flow and the viscosity in that state of the flow is called the steady flow viscosity).

If , then has the units of time and is a useful measure of the response time of the material's viscoelastic response.

In a series connection such as the Maxwell model, the stress on each element is the same and equal to the imposed stress, while the total strain is the sum of the strain in each element:andwhere ‘s’ and ‘d’ stand for the spring and dashpot respectively.

In seeking a single equation relating the stress to the strain, it is convenient to differentiate the strain equation and then write the spring and dashpot strain rates in terms of the stress:

Multiplying throughout with μ, and usingwe getSince we are dealing with the mucous layer, we put instead of μ.

Thus the final equation isHere, is the stress tensor for visco-elastic effect. the strain rate for visco-elastic effect. the 0.034 s (single relaxation time, which is the time required for the ciliary beat cycle to reduce to 1/e of its initial intensity). the 0.0482 Pa s.

The equations of the Maxwell model are linearized under the assumption that the penetration of the cilia into the mucous is small. Due to the simplification of the Maxwell model, which is permitted due to the consideration of only very short-lived deformations, this model does not consider effects such as shear thinning.

The mucous viscosity of 0.0482 Pa s was obtained through oscillatory testing. Two constants are found: , the storage modulus (associated with elasticity) and the loss modulus (associated with viscosity). Fung [33] describes how these constants can be related to the Maxwell constitutive equation. The strain on the mucus will be of the form .

For small oscillations the stress will respond sinusoidally, with a phase difference . Then . Noting that and we can substitute into Maxwell's constitutive equation to show that A set of values for and (measured over a wide range of will lead to a widely varying set of values for and ). One cannot fit the simple Maxwell model to a range of real experimental results. However, if one chooses a characteristic frequency of the system, 5–10 Hz, one can find appropriate values for and .

In the study by Lutz et al. [34], for canine tracheal mucus, values obtained were =1 Pa and =0.64491 Pa (at the frequency of 7.2497 Hz). These values correspond to the parameters λ=0.034 s−1 and =0.0482 Pa s. Since 0.001 Pa s, we have =48.2.

Numerical procedure

The numerical technique and implementation on COMSOL

Numerical software (COMSOL 3.4) is used to solve the model. COMSOL is a ‘Finite Element Analysis’ solver and simulation software for engineering purposes. The Maxwell model of viscoelasticity is established in the traction and mucous layers using the Partial Differential Equation (PDE) mode on the software. The relative tolerance was set to 0.01 and absolute tolerance as 0.001. An ‘extra fine’ mesh consisting of 3309 triangular elements and 27,645 degrees of freedom was created. Using the relation , the time step was fixed as 0.01 s with a maximum solution time of 5 s (approximately 50 cilia beat cycles).

Transient analysis mode was used. A time dependent linear solver, Direct PARADISO was used. The parallel sparse direct linear solver PARDISO works on general systems of the form A x=b. In order to improve sequential and parallel sparse numerical factorization performance, the solver algorithms are based on a Level-3 BLAS update, and they exploit pipelining parallelism with a combination of left-looking and right-looking super-node techniques. The Nested Dissection preordering method was used.

The consistent initialization of the DAE systems was Backward Euler and the error estimation strategy was to exclude algebraic terms. Also, the constraint handling method was elimination.

Boundary and interface conditions

Boundary conditions

The external boundaries were defined as open boundaries except for the bottom-most boundary and the top-most boundary (Fig. 3(a)). The epithelial layer spans a very large area within the body, lining various organs, and so mucous flow occurs over relatively large distances. However, our domain represents only a small part of the epithelium. Hence, is imperative to allow for free flow of fluid to and from regions adjacent to our domain. We thus use the open boundary condition at the sides of the domain. The bottom boundary was defined as a wall. It represents the fact that there is no flow through the epithelium and the no-slip condition at the epithelium-PCL interface. Also, due to the wall boundary condition, it is possible to impose the zero vertical velocity boundary condition at the epithelium-PCL interface. The top most boundary is a stress boundary. The stress boundary represents the interface between the mucus and the air and allows us to model the shearing effects of the air flow upon the mucus. In our model, however, we have used zero shear stresses in order to simplify the situation. All the interior boundaries had zero vertical velocity applied on them. This is because we assume in this model that there is no exchange of fluid between layers.

Subdomain conditions

The density of each of the layers was kept at 1000 kg/m3, which is the same as water. The absolute viscosity of the PCL was 0.001 Pa s. Since the viscosity of traction and mucus layers are given by the Maxwell model, we couple the Navier–Stokes module with the Partial Differential Equation (PDE) module. As mentioned earlier, the PDE mode is used to define the viscous stresses in the traction and mucous layer through the Maxwell model. The ciliary resistive forces are modeled as bulk forces working in the PCL layer and the propulsive forces of the cilia are modeled as bulk forces in the traction layer.

The reference paper(Smith et al. [28]) gives a complex mathematical description of this three layer model of the ASL. Therefore, in the later part of our study, when we study patches of cilia, the advantage of our numerical treatment over the analytical solution given in the reference is obvious. By using the numerical solution with COMSOL, we avoid the complex analytical treatment of the interfaces between the ciliated and unciliated regions. Since at these interfaces, it is difficult to prescribe any concrete boundary conditions, we describe them as continuity boundaries on COMSOL. Also, for the analytical solution of the simple domain as shown, the total number of variables required are 30 and also 30 boundary conditions (Smith et al. [28]). Understandably, it is of greater efficiency to use a numerical solution when and equal number of variables and boundary conditions are going to be required for each patch comprising the multiple patch model shown later, which is the main aim of our efforts.

We took the mesh with 2444 triangular elements and 11,226 degrees of freedom for our simulations. Fig. 3(b) shows the mesh used. In order to check the spatial grid independence, different mesh sizes were taken successively and for each, separate simulations were created. The maximum velocity for each simulation was noted and it was found that among the different mesh sizes, the value of this maximum velocity did not vary considerably, with a maximum deviation of about 1%.

Results and discussion

Validation

The validation of the model was carried out by comparing with the velocity profiles as shown in the reference paper Smith et al. [28] (which develops an analytical model of the three layer structure of the ASL).We were also able to match our values of mucous velocities with those found in the reference. The results in this reference are based on their analytical solution of the problem. However, their results have been proved against previous experimental work. Therefore, our validation is not only against the mathematical solution of Smith et al. [28] but also against various experimental studies.

The average value of mucus velocity given in Smith et al. [28] is 40 μm/s. Also, in Salathe et al. [35], a range of values between 67 μm/s and 333 μm/s are reported, greater accuracy being on the first value. In our simulation, the mucous velocity is of the order of microns/s up to tens of microns/s, varying slightly through the domain. Thus, the values lie sufficiently within the realistic range.

Fig. 4(a) and (b) represent the horizontal velocity at different levels in the fluid. Fig. 4(a) represents the velocity profile at the vertical level of the PCL-traction layer interface while 4(b) represents the velocity profile at the vertical level of the mucous surface.

Fig. 4

(a) and (b) These figures compare the horizontal velocity profiles at different levels from the epithelium above. We show the similarity between the profiles shown in the reference paper and that obtained from our simulation. (a) is for the vertical level of the interface between the PCL and traction layer while (b) is for the vertical level of the mucous surface. The solid line curves represent data from the reference whereas the dotted line curves represent the results obtained from the simulation. (c), (d) and (e) compare the horizontal mucous velocities at x = 7, 25 and 30 μm respectively. The solid line curves represent data from the reference whereas the dotted line curve represents the results obtained from the simulation.

We have also compared our velocity profiles to that given in the reference paper at different positions across the domain. We have shown three plots for the horizontal velocities to indicate that the velocity profiles obtained by our simulations are largely similar in nature to that in the reference paper (Fig. 4(c–e)).These profiles are also observed to be similar to that shown in Fulford and Blake [36] where a mathematical treatment of a two layer model of the mucociliary epithelium has been presented.

Behavior of cilia bunches

Description of the model consisting of ciliary bunches

The main objective of this study was to simulate the interplay of adjacent bunches of cilia, as opposed to a dense mat representation of the cilia. Therefore we observed the change of net mucous flow upon changing the frequency, spacing and phase lag of the beating of one cilia bunch with respect to the beating of adjacent bunches. In order to vary the frequency and phase lag of the ciliary movement, we altered those terms in the harmonic functions that define the ciliary velocity and propellary force.

The study involved placing individual units created in the first part of the study adjacent to each-other. The boundaries separating these units are continuity boundaries and thus allow flow from one unit to the other ( Fig. 5). We observe movement of the mucous to the positive direction in each of the domains (Fig. 5). In all our simulations, we have five such units. The first, third and fifth are un-ciliated, whereas the second and fourth are ciliated. The black line separating the adjacent units can be seen in the Fig. 5(b).Thus the second and fourth units of the model are the first and second ciliated domains respectively. Also, the third unit of the model is the second un-ciliated domain. In our studies, we focus on the first and second ciliated domains and the second un-ciliated domains (as this is the domain separating the two ciliated regions).When we vary the spacing between the two ciliated domains; we alter the length of the second un-ciliated domain. The study of the phase difference is conducted by altering the phase difference between the ciliary beat cycles in the first and second ciliated domains. Also, it is to be noted that it is the mucous layer in each unit is that of interest since we aim to study the factors altering mucous movement in order to gain deeper insight into the pathological conditions that may occur in the ciliary epithelium. We therefore study the movement of the mucous as a result of the ciliary motion in the PCL layer.

Fig. 5

(a) A sketch of the model comprising bunches of cilia. Two bunches of cilia are shown separated from each-other by un-ciliated domains. The nature of the boundary conditions are shown. (b) A snapshot of the domain at 4 s. As seen in (a), there are two ciliated domain separated by an un-ciliated domain and with two un-ciliated domains on each side. The surface plot shows the velocity field and the ciliated domains can be recognized by the greater values of velocity field within them, denoted by the red coloration. (c) A snapshot of the domain at 4.2 s. One can notice movement of the densely red regions of the velocity surface plot towards the right. (d) A snapshot of the domain at 4.5 s. One can notice movement of the densely red regions of the velocity surface plot further towards the right. (e) A snapshot of the domain at 5 s. One can notice movement of the densely red regions of the velocity surface plot further towards the right. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

First, in order to study mucous-motion due to one ciliated bunch (as a reference to see the effect of one cilia bunch on the surrounding bunches), a model with only one ciliated domain at the center, and un-ciliated domains on both sides was created ( Fig. 6). From this, it was observed that:

Fig. 6

(a) The above snapshot is taken from a simulation in which only the center unit was ciliated at solution time 4 s. The purpose of this simulation was to see the effect of an individual bunch on the neighboring units. (b) The above snapshot is taken from a simulation in which only the center unit was ciliated at solution time 4.5 s. (c) The above snapshot is taken from a simulation in which only the center unit was ciliated at solution time 5 s.

The horizontal component of velocity (x-velocity) shows low variation in magnitude throughout the domain at any given time. However, after regular intervals of time, there is a surge of high horizontal velocity in the mucus in the positive direction, indicating the effective stroke of the ciliary beating.

The vertical velocity component is some orders of magnitude greater at the ciliated domain than the other domains at any given time. The pattern of vertical velocities throughout the domain varies relatively less with time compared to the horizontal velocities.

The net velocity field is positive all through the domain, throughout the solution time. The magnitude is greater in the ciliated domain relative to the remaining domains. The ciliary action is demonstrated by the characteristic vortical region seen in the PCL, traveling within the ciliated domain periodically.

Observations made

Up till now, we have shown and explained our model, laying stress on the salient points. Now we show our observations of the effects of varying the three different parameters (frequency; phase difference and spacing as mentioned earlier) in various combinations. We look at their effect on the net mucous flow. The net mucous flow is estimated by calculating the subdomain integration in the top most layer in a given unit (i.e. the mucous layer).

The findings of our simulations are summarized as follows:

Effect of frequency

The ciliary beat frequency (CBF) variation is implemented by changing the frequency of the harmonic functions that describe the cilia velocity and ciliary forces in the traction layer.

The average mucous flux was calculated with both ciliated domains at the same frequency. The frequencies implemented were 10 rad/s up to 150 rad/s, with increments of 10 rad/s. The fluxes in the first ciliated domain were noted. The observations were made corresponding to each frequency, at 3 s till 3.2 s of solution time, with increments of 0.01 s. The 20 values thus obtained for each frequency were averaged and plotted. Fig. 7(a)–(c) shows the averaged value at each frequency. The ciliated domains were at 0° phase difference in Fig. 7(a), 60° phase difference in Fig. 7(b) and 120° phase difference in Fig. 7(c) respectively.

Fig. 7

The above figures show the average mucous flux at the first ciliated domain in the simulations with both ciliated domains at the same frequency of 10 till 150 rad/s. The fluxes (y axis) were found in the first ciliated domain at each frequency (x axis) for times 3 s till 3.2 s with a time step of 0.01 s. The 20 values thus obtained for each frequency were averaged and plotted. (a), (b) and (c) show the averaged value at each frequency at 0 phase degree phase difference, 60° phase difference and 120° phase difference respectively.

In each of the simulations performed, both the ciliated domains were in synchronization. They were both of equal size and were at the same frequency. Thus, the mucous fluxes in either of them were almost identical. Thus it was sufficient to note the values in either one of them.

It is seen from the plots (Fig. 7) that the average values of mucous flux (i.e. average values of subdomain integration of the velocity field) show peak values at approximately 10 rad/s (1.6 Hz). From 10 rad/s onward to higher frequencies, the average mucous flux decreases almost along a hyperbolic curve. This may be indicative of optimal mucous flows around this value of frequency. Up till10 rad/s, we did not observe any particular trend and due to the erratic nature of both the values of average mucous flux and the trend shown, we ignored this range of frequency as ‘noise’.

Effect of phase difference

The effect of phase difference is studied by placing the adjacent ciliated domains either in complete synchronization (0° phase difference) or completely out of synchronization (180° phase difference).

Implementing 180° phase difference causes drastic reduction in magnitude of the mucous flux in the ciliated domains, as noticed in the figures ( Fig. 8, Fig. 9, Fig. 10, Fig. 11, Fig. 12, Fig. 13).

Fig. 8

Dashed line: first un-ciliated domain; dashed line with plus signs: first ciliated domain; dash-dotted line: second un-ciliated domain; dotted line: second ciliated domain; solid line: third un-ciliated domain. (both ciliated domains at 30 rad/s, with 0 phase difference and equal spacing between all units).

Fig. 9

Dashed line: first un-ciliated domain; dashed line with plus signs: first ciliated domain; dash-dotted line: second un-ciliated domain; dotted line: second ciliated domain; solid line: third un-ciliated domain (both ciliated domains at 30 rad/s, with 180° phase difference and equal spacing between all units).

Fig. 10

Dashed line: first un-ciliated domain; dashed line with plus signs: first ciliated domain; dash-dotted line: second un-ciliated domain; dotted line: second ciliated domain; solid line: third un-ciliated domain. (both ciliated domains at 30 rad/s, with 0 phase difference and increased length of first ciliated domain).

Fig. 11

Dashed line: first un-ciliated domain; dashed line with plus signs: first ciliated domain; dash-dotted line: second un-ciliated domain; dotted line: second ciliated domain; solid line: third un-ciliated domain (both ciliated domains at 30 rad/s, with 180 phase difference and increased length of first ciliated domain).

Fig. 12

Dashed line: first un-ciliated domain; dashed line with plus signs: first ciliated domain; dash-dotted line: second un-ciliated domain; dotted line: second ciliated domain; solid line: third un-ciliated domain (both ciliated domains at 30 rad/s, with 0 phase difference and increased length of second un-ciliated domain).

Fig. 13

Dashed line: first un-ciliated domain; dashed line with plus signs: first ciliated domain; dash-dotted line: second un-ciliated domain; dotted line: second ciliated domain; solid line: third un-ciliated domain (both ciliated domains at 30 rad/s, with 180 phase difference and increased length of second un-ciliated domain).

For a particular frequency, when two adjacent ciliated domains are put at 180° phase difference from one another, the mucous flux for the ciliated domains are lower than the case where the ciliated domains are in synchronization.

Variation of inter-bunch spacing

The variation of inter-bunch spacing is implemented by changing the length of either the ciliated domain or un-ciliated domain. By changing the length of a ciliated domain, we alter the distance between two un-ciliated domains. Similarly, by changing the length of an un-ciliated domain, we alter the distance between two ciliated domains. In our study, we observed the effect of changing the distance between two unciliated domains by doubling the length of the first ciliated domain. In the simulations where the effect of increasing the un-ciliated domain was studied, we doubled the length of the second un-ciliated domain.

Increase in length of ciliated domain

Increasing the length of the ciliated domain increases the mucous flux in that domain drastically. This is true for both when the two ciliated regions are at 0° and 180° phase difference. The variation is observed to be almost linear (Fig. 10, Fig. 11)

The mucous flux for the second ciliated domain (whose length has not been altered) is noticeably less when there is 180° phase difference between the larger ciliated domain and itself. However, the mucous flux in the larger domain is not reduced in this situation, compared to when the two ciliated regions are in phase.

The mucous flux for the third un-ciliated domain (last unit from the left) also decreases considerably when 180° phase difference is applied.

Increase of un-ciliated domain

The length of the un-ciliated region between the first and second ciliated domains is doubled. When there is 0° phase difference between the ciliated domains on either side, the mucous flux for this un-ciliated region is larger than other un-ciliated domains. However, the increase in mucous flux due to increased length is not linear (unlike in the case of ciliated domains).

Placing the two ciliated domains at 180° phase difference reduces the mucous flux in this unciliated domain greatly. Even though its length is now double that of the other un-ciliated domains, the value of mucous flux here does not differ significantly from that of other the un-ciliated domains.

Again, the mucous flux for the third un-ciliated domain also decreases considerably when 180° phase difference is applied between the two ciliated domains (Fig. 12, Fig. 13).

The results stated above are shown diagrammatically in the Fig. 8, Fig. 9, Fig. 10, Fig. 11, Fig. 12, Fig. 13.The green line in the figures (Fig. 8, Fig. 9, Fig. 10, Fig. 11, Fig. 12, Fig. 13) represents the mucous flux (subdomain integration for the top-most layer i.e. mucous) in the first un-ciliated domain. Dark blue is for the mucous flux in the first ciliated domain; red for the second un-ciliated domain; magenta for the second ciliated domain and black for the third un-ciliated domain.

The Fig. 8, Fig. 9, represent the situation where all the units (both ciliated and un-ciliated) are of the same length. Fig. 10, Fig. 11, represent the situation where the length of the first ciliated domain is increased. Finally, Fig. 12, Fig. 13 represent the situation where the length of second un-ciliated domain is increased. Further, in Fig. 8 there is 0° phase difference between the two ciliated domains. In Fig. 9 there is 180° phase difference between the two ciliated domains. Similarly the pairs of Fig. 10, Fig. 11, Fig. 12, Fig. 13 show 0° phase difference and 180° phase difference respectively, for their corresponding situations of inter-bunch spacing. Also, in all these figures, the ciliated domains are beating at 30 rad/s.

Fig. 8, Fig. 9 (All domains of equal length, ciliary beating at 30 rad/s in both ciliated domains. Fig. 8 shows the situation where both ciliated domains are in phase whereas Fig. 9 represents the situation where they are out of phase).

We see that the dark blue and magenta lines almost coincide in Fig. 8. This indicates that at 0° phase difference, both the ciliated domains behave almost identically. However, the magenta line in Fig. 9 goes down drastically, revealing that at 180° phase difference, the mucous flux reduces in the ciliated domains. In both Fig. 8, Fig. 9, the un-ciliated domains show almost equal mucous fluxes with each other, and there is not much difference between the magnitudes in Fig. 8, Fig. 9. However, the black line i.e. mucous flux in the third un-ciliated domain reduces for 180° phase difference.

Fig. 10, Fig. 11 (Increased length of first ciliated domain, ciliary beating at 30 rad/s in both ciliated domains, Fig. 10 shows the situation where both ciliated domains are in phase whereas Fig. 11 represents the situation where they are out of phase).

The mucous flux in the first ciliated domain has increased immensely due to its increased length. In fact, it can be observed that the increase in mucous flow in this domain varies linearly with the change in length. When the length of this domain was doubled, the mucous flux doubled to 3×10−15 from approximately 1.4×10−15 in case of the original dimensions represented in Fig. 8, Fig. 9.When 180° phase difference is applied (as shown in Fig. 11), the mucous flux in the second ciliated domain reduces drastically, whereas that in the first ciliated domain is not much different from what we see in Fig. 10. The un-ciliated domains show the similar trend as we observed in Fig. 8, Fig. 9.

Fig. 12, Fig. 13 (Increased length of second un-ciliated domain, ciliary beating at 30 rad/s in both ciliated domains, Fig. 12 shows the situation where both ciliated domains are in phase whereas Fig. 13 represents the situation where they are out of phase).

We see that the two ciliated domains show almost identical mucous fluxes when in 0° phase difference (as evident from their over-lapping plots in Fig. 12). In Fig. 13,(where there is 180° phase difference between the ciliated regions), the magenta line goes down drastically. In Fig. 12, the red line indicating the mucous flux in the second un-ciliated domain is higher up than the lines for the other two un-ciliated domains. However, it reduces greatly in Fig. 13, where it becomes almost equal to the mucous flux values in the other two un-ciliated domains. In the previous plots (Fig. 8, Fig. 9, Fig. 10, Fig. 11), the mucous flux in the second un-ciliated domain, located between the two ciliated domains, did not vary significantly when the ciliated regions were put out of phase. Thus, with increasing distance between cilia bunches, it becomes increasingly necessary for the cilia to be in phase.

Insights obtained from our simulations

The simulations that represent the cilia to occur in patches are physiologically more accurate than those that represent the ciliary epithelium to be an infinite array of cilia, beating in synchrony. It was interesting to compare the average mucus velocity obtained from our ‘validation model’ (where cilia are modeled as a ‘dense mat’, as per the mathematical model by Smith et al.), and our new model with cilia represented in ‘patches’. We plotted the average mucous velocity, averaged over time, at different values of frequency (from 10 rad/s to 100 rad/s) for the validation model. We also plotted the average mucus velocity (averaged over both ciliated and un-ciliated patches) in our simulations with cilia in patches, averaged over time. In these simulations, both the patches were kept at the same frequency (which was varied from 10 rad/s to 100 rad/s). The obtained plots are shown in Fig. 14. We took separate sets of readings for when the patches were in 0° phase difference (represented by the blue line on Fig. 14), 60° phase difference (green line on Fig. 14), and 120 degree phase difference (black line on Fig. 14) from each-other, respectively. The set of readings obtained from the validation model (dense mat representation of cilia) is shown by the red line.

Fig. 14

Dashed line: single domain (dense mat); dash-dotted line: 0° phase difference between patches; dotted line: 60° phase difference between patches; solid line: 120° phase difference between patches. Comparison between mucus velocity obtained for ‘dense mat’ representation of cilia and ‘patched’ representation of cilia.

It can be seen that the values for all the plots are of the same order of magnitude and that the values obtained for our patched model do not vary appreciably from that of the validated ‘dense mat’ (or single domain) model. However, the values in the ‘dense mat’ model are about three times higher than that of the patched model. This difference, however, can be attributed to the obvious fact that in the patched domain, there are less cilia per unit length to propel the mucus. In fact, the reduction in mucous movement is proportional to reduction in ciliated space. Over a considerable large surface area however, such as the entire area of epithelial surface lining the respiratory system of an organism, this effect will not be prominent.

Therefore, the patched model corresponds well with the previously validated model, the ‘dense mat’ model of cilia. However, it is more realistic, since it incorporates the true physiology of the ciliary epithelium i.e. the fact that cilia exist in patches, separated by either clefts or by other types of cells (such as goblet cells and Clara cells). Due to this fact, the patched model also allows us to dig deeper in to the mechanism of mucus clearance. It allows us to analyze a factor which has hitherto not been stressed upon i.e. how the cilia, in different locations of the epithelium assist each-other in muco-ciliary clearance.

Another fact to be stressed upon is that cilia operate in metachronal waves and the dynamics of this kind of wave are complex. This adds complexity to the study of the resulting mucous motion. Furthermore, mucus is a complex visco-elastic fluid. Therefore, it is difficult to predict how the mucous flow will be affected upon varying the parameters that define inter-bunch interaction. There are several questions that need to be answered. Do all the patches beat in the same direction at all times? What happens if some of the patches stop beating due to disease or even during normal conditions? In those situations, do the remaining patches change their mode of beating in order to keep the mucous flow rate undisturbed? Through our simulations, we try to find the answers to some of these questions.

We notice from our results that the interplay of multiple patches on net mucous flow can be largely explained by the superposition principle. For example, the reduction of mucous flux when the adjacent cilia bunches are put out of synchronization could perhaps be due to the destructive interference between them.

One would expect that by increasing the length of a ciliated domain, the mucous flux in the ciliated region would increase. This is indeed observed in our simulations. It is further found that by doubling the length of the ciliated domain, the mucous flux in this region also almost doubles, indicating an almost linear variation of mucous flux with length of domain. Since by increasing the dimensions of a ciliated domain we are increasing the region of activity of the cilia, this is quite understandable. However, after increasing the length of the ciliated domain, the effect of putting the bunches out of phase seems to be largely subverted. This is evident from the fact that in Fig. 11, we see that the value of mucous flux in the first ciliated domain (whose length has been increased) is almost identical to that seen in Fig. 10, whereas the mucous flux in the second ciliated domain seen in figure is much less in Fig. 11 compared to Fig. 10. This could mean that larger patches have the advantage of being able to maintain their mucous flow irrespective of the activity within surrounding patches.

It is slightly more difficult to predict the effect of increasing the gap between the cilia patches. In our simulation, this was done by increasing the length of the unciliated region in between the two ciliated domains. In Ref. [30] it is mentioned that the region of ciliary activity is divided into distinct patches and the influence of each such patch extends only the short distance across the area and fades at the edges of the region. Owing to the increased distance between the ciliated regions therefore, we would expect the mucous flow within the ciliated regions to not alter significantly and that the flow in the intermediate unciliated region would reduce, since the effect of the ciliated bunches would not be as pronounced in this larger bulk of mucous. It is seen, however, that despite the increased distance between the bunches, the mucous flow in the unciliated region is not reduced. When the cilia bunches on either side of this unciliated domain are in synchrony, for example, the flux in this unciliated region is larger than that of the other unciliated regions. This indicates that even if sections of ciliary epithelium are impaired due to disease, the muco-ciliary clearance rates can remain conserved provided surrounding patches are functioning. Of course, this argument fails in case of large scale damage of ciliary epithelium. Though the mucous flow in this unciliated region is much greater than that of the other unciliated regions, the increase is not linear with increase in length, as was observed when the length of the ciliated domain was increased. More interestingly, however, when the ciliated patches are completely out of phase, the mucous flux in the intermediate unciliated region greatly reduces and becomes almost equal to that of the other regions without cilia. This is interesting because in the cases where all the domains were of equal length and also in the other case where the length of the first ciliated domain was increased, the flow in the intermediate unciliated region was not affected much when the bunches were put out of synchronization. It seems therefore, that when the ciliated regions are farther away from each-other, the effect of putting them out of synchronization becomes more pronounced.

Frequency is another important factor that alters mucous flux. It is known that the frequency of ciliary beating can be altered due to the influence of chemicals such as Ach. (Kawakami et al. [37]). In a study to investigate the effects of pharmaceutical agents on the ciliary beat frequency (on rabbit tracheal epithelium), it was found that benzalkonium chloride induced inhibition of CBF in a concentration dependent manner (Wang et al. [38]). On the other hand, potassium sorbate resulted in an increase in CBF (Wang et al. [38]). Also it has been found that the hormone Progesterone has the effect of reducing ciliary beat frequency in experimental studies on mouse fallopian tube (Bylander et al. [39]). Thus, it is clear that the ciliary beat frequency can be modulated by extraneous agents, which in turn alters the flow of mucus. This can have therapeutic application in treating various muco-ciliary disorders. Recognizing this fact, we made a detailed study of how modulation of CBF of patches of cilia affects mucus flow.

While testing the effects of frequency, we also put the two cilia bunches at different phase differences from each-other, in order to provide a greater range of trial conditions. For each value of phase lag applied, we put both the bunches at the same frequency. We applied frequencies between 10 rad/s and150 rad/s, with increments of10 rad/s. Interestingly, for all the plots we obtained, the peak value of mucous flux was seen at the frequency value of 10 rad/s (1.6 Hz). It is not easy to provide an explanation for this observation. However, the fact that we have obtained a probable optimum frequency is interesting. It is highly probable that the value of optimum frequency we have obtained is highly specific to our simulation, with the given arrangement and shape of cilia bunches. Nevertheless, it does make us aware that there is a possibility of something such as an optimum frequency of ciliary beating. This is an interesting aspect deserving greater attention in the future. The observations mentioned above are summarized in Table 1.

Table 1

Summary of observations.

Conditions applied	Effects observed
The two bunches are out of phase	Reduced mucous flow.
Increased length of ciliated region (0° phase lag with the other bunch)	More mucous flow in ciliated region
Increased length of ciliated region with adjacent one out of phase (180° phase lag)	Larger mucous flow in the one with increased length, decreased flow in the adjacent bunch
Increased length of unciliated region(0° phase lag)	Increased flow in this larger unciliated region
Increased length of unciliated region(180° phase lag)	Flow in this larger unciliated region greatly reduced, becomes almost equal to that of other unciliated regions
Frequency variation from 10 rad/s to 150 rad/s with increments of 10 rad/s	Maximum flow at 10 rad/s

Summary of the above is given as follows:

Our model is a two dimensional representation of the ciliary epithelium. Though the expressions for the ciliary velocity and ciliary forces do bring out the essential features, such as periodicity of beat cycle, the linear variation of ciliary velocity with distance from the epithelium and the resistance the cilia create for PCL motion, there are certain factors that could not be taken into account due to the two-dimensional nature of our model. For example, the bunches of cilia depicted in our model have regular geometry. In reality, the patches of cilia are of irregular shape. Also, in reality there is a spatial distribution of the patches. Thus, interaction of the ciliary bunches does not take place only along one dimension, and the analysis is more complex. Our simulation has been based upon a three-layer model presented in Smith et al. [28]. We have validated against this paper and shown that certain important parameters and trends observed in our model are adequately close to the values given by the reference paper.

Conclusion

We have carried out a computer simulation of the interplay of ciliary patches. The patches have been modeled based on the three layer model of the ciliary epithelium. The effects of varying the frequency, spacing and phase lag of the beating of one cilia bunch with respect to the beating of adjacent patches have been studied. The interplay of multiple patches and the resulting mucous flow is largely based on the superposition principle. By increasing the length of a ciliated domain, the mucous flux increases and there is almost linear variation of mucous flux with length of domain. The mucous flow in unciliated regions increases when ciliated regions beat in synchrony but the variation is not linear with increase in length. However, when the ciliated patches are completely out of phase, the mucous flux in the intermediate unciliated region greatly reduces. The peak value of mucous flux was seen at the frequency value of 10 rad/s (1.6 Hz).

Future studies based on our results will bring out more on the study of cilia and ciliopathies. Our effort is expected to be contributory towards the understanding of cilia action and especially the coordination of ciliated regions on an epithelial layer. It thus helps gain a deeper understanding of the pathology resulting due to malfunctioning of cilia and lack of coordination among them.

Conflict of interest

The authors hereby state that there is no conflict of interest.

Mesh size (Number of mesh elements)	Maximum velocity seen
1500	1.623×10−4
2444	1.635×10−4
4640	1.638×10−4