New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Cilia internal mechanism and metachronal coordination as the result of hydrodynamical coupling
Abstract
We present a simple but realistic model for the internal bendgenerating mechanism of cilia, using parameters obtained from the analysis of data of the beat of a single cilium, and incorporate it into a recently developed dynamical model. Comparing the results to experimental data for twodimensional beats, we demonstrate that the model captures the essential features of the motion, including many properties that are not built in explicitly. The beat pattern and frequency change in response to increased viscosity and the presence of neighboring cilia in a realistic fashion. Using the model, we are able to investigate multicilia configurations such as rows of cilia and twodimensional arrays of cilia. When two adjacent model cilia start beating at different phase, they synchronize within two cycles, as observed in experiments in which two flagella beating out of phase are brought close together. Examination of various multicilia configurations shows that metachronal patterns (i.e., beats with a constant phase difference between neighboring cilia) evolve autonomously. This provides modeling evidence in support of the conjecture that metachronism may occur as a selforganized phenomenon due to hydrodynamical interactions between the cilia.
Ciliary motion and particularly the metachronism phenomenon have attracted a great deal of research effort both experimentally and theoretically. Metachronal coordination between cilia is a situation where cilia beat together with a constant phase difference between adjacent neighbors, their tips forming a moving wave pattern. The reason why arrays of cilia beat in a metachronal pattern is not fully understood. The work of Machemer (1), for example, shows that membrane voltage and calcium levels affect the direction of the metachronal wave and also the directions of the effective and recovery strokes of the ciliary beats. Some researchers speculate that the intriguing metachronism phenomenon is possibly the result of hydrodynamical coupling (e.g., see refs. 2–4). This work provides a theoretical model that supports this conjecture.
Gueron and Liron (ref. 5; GL hereafter) introduced an improved technique for describing the hydrodynamics of moving cilia, based on a refined slender body theory, which offers an alternative for the simplistic resistive force theory (known as the Gray and Hancock approximation) that relates drag forces and velocities. As a result, the accuracy and consistency of the model for cilia beating is markedly improved. More important, the GL equations provide a method for dynamical simulations of multicilia configurations that account for the effects of neighboring cilia and the effect of the surface from which the cilia emerge. This model was originally applied to a twodimensional setup and later extended to describe threedimensional beating (6).
The work of GL was mainly oriented toward developing the framework for dynamical modeling of multicilia configurations. For the internal mechanism of a cilium (hereafter referred to as the “engine”) they used an ad hoc equation for the active (normal) shear force generated inside a cilium. However, the GL model has the engine as a separate building block which can be replaced (and is replaced here) by any other engine. For simplicity, the GL engine included a builtin frequency term that controls the resulting beat frequency and predetermines the duration of the effective and the recovery strokes. Therefore, this model does not allow for realistic changes in the beat frequency in response to external load such as increased viscosity or external flow generated by neighboring cilia. The GL engine cannot, therefore, be used for investigating metachronism. Murase (7) proposed an engine which is analogous in its form to muscle sliding models. However, his model is restricted to smallamplitude motion, which is not the case for real cilia, and is based on the Gray and Hancock approximation, which does not allow for modeling multicilia configurations. This engine gives poor results when incorporated into the GL model.
We present a simple plausible functional form for the engine that is related to the dyenin arms and the radial spokes system (see refs. 8 and 9). This engine has a small number of parameters. To obtain realistic values for the parameters we use observations of ciliary motion to compute locations and velocities during the beat cycle, solve the GL equations to compute the drag forces, and then compute the internal shear forces. These forces are then used to determine the engine parameters. The computed engine is tested by incorporating it into the GL dynamical equations for twodimensional beats. The model produces realistic beats, and it reproduces experimental results such as exponential decrease in beat frequency with increased viscosity, selfsynchronization between two adjacent cilia, and frequency matching with the frequency of external flows. Finally, we investigate multicilia configurations and find that metachronal patterns evolve autonomously due to the hydrodynamical interaction between the cilia. This supports the conjecture that metachronism can be explained as the result of hydrodynamical coupling.
The Model Equations
The drag–velocity relation, combined with the twodimensional curve geometry, yields the following three integro partial differential equations, written in nondimensional form (see ref. 5): 1 2 3 We use boundary conditions that correspond to a cilium that is pinned at its basal end and free at its distal end, and initial conditions that correspond to an erect cilium: 4 In Eqs. 1–4 s is the arclength variable of the centerline of the cilium, measured from the basal end, t is time, and the functions depend on s and t (we omit writing this dependence explicitly when it is clear from the context). The subscripts s and t denote the space and time derivatives, respectively. The function α = α(s, t) measures the angle between the tangent to the cilium and the horizontal axis. α is related to κ = κ(s, t), the curvature of the cilium’s centerline, by the equation κ = α_{s}. F_{N} and F_{T} are the normal and tangential components of the shear force generated inside the cilium, respectively. S = S(s, t) is the active shear force due to the internal sliding filaments mechanism and the radial spokes system. The model equation representing S is derived below. S_{0} is a typical fixed magnitude of the internal shear force (used for scaling), L is the length of the cilium, and E_{b} is its elastic bending resistance. The expressions g_{T}, g_{N} involve integrals of the appropriate singular solutions of Stokes flow, C_{N} and C_{T} are the related GL coefficients, C_{TN} = C_{T}/C_{N}, and C_{NT} = C_{N}/C_{T}. These equations, derived in detail by GL, are outlined in the Appendix.‖
Because we deal here with a twodimensional problem, the angle (α) determines the centerline curve if the position and the orientation at s = 0 (i.e., the anchor of the cilium) are given.
The normal shear force F_{N} in Eq. 1 includes the contribution due to the elastic resistance of the cilium and the active shear force S(s, t). Clearly, S(s, t) determines the dynamics of the cilium in the following way: If the position of the cilium (α) is given at time t, one can calculate F_{N} from Eq. 1. F_{T} can then be obtained from Eq. 2, which relates the normal and the tangential components of the shear force (due to the inextensibility of the cilium). α is propagated in time by means of Eq. 3.
A LoadDependent Internal Engine
In this section we briefly describe our method for modeling the internal mechanism of cilia. As data, we use the observed (twodimensional) beat pattern of the cilium of Paramecium (diagrams and discussion in refs. 10 and 11). We measure the coordinates of equally spaced points along the cilium at different times during the beat cycle. The velocity of the cilium is calculated from these data (after proper smoothing). Using the calculated velocity, we solve the integral Eqs. 1–3 for the drag forces as unknowns, and then compute the active shear force S(s, t).
According to Sleigh and Barlow (9), the radial spokes system contributes to the total shear force at the bent region but not at the straight region. A reasonable intrinsic (and local) indicator of the place of the bend is the curvature κ(s, t) = α_{s}(s, t). Therefore, we fit the measured function S(s, t) to a simple combination of functions of the curvature, the arclength s and the deviation from the resting position: 5 Here, ω (a typical angular velocity), C̄_{N} = C_{N}L^{2}/S_{0}, and A_{1}, A_{2} are parameters used for scaling. ω, A_{1}, A_{2} have different values during the effective and the recovery strokes. To simulate the fast initiation of the bend at the beginning of the recovery stroke, we take A_{1} in the region 0 ≤ s < 0.1 to be twice as large as its value in the region 0.1 ≤ s ≤ 1. Altogether, the effective stroke is determined by four parameters (C̄_{N}ω, A_{1}, A_{2}, B) and the recovery stroke by five (A_{1} has two values). Note that C̄_{N} and ω appear only as the product C̄_{N}ω but are artificially separated to introduce a realistic frequency term (ω) and for consistency with the model equations developed by GL.
The first term in Eq. 5 represents the contribution of the dynein arms, and the second curvaturedependent term represents the contribution of shearresistant elements such as the radial spokes and nexins. Note that during the effective stroke the cilium is almost straight (κ ∼ 0), and thus only the dynein arms contribute to the motion.
Our engine has a similar functional form during the effective and the recovery stroke but with opposite signs (+ during the effective stroke and − during the recovery stroke). This is based on the evidence (10) that different groups of filaments are active at different stages of the beat cycle. To complete the definition of the model we must specify when the switching between the two phases of the motion occurs.
We speculate that the effective stroke continues as long as the filaments can slide in this direction, that is, until the effective shear σ between filaments reaches a certain maximal value (note that the effective shear is a geometric quantity, not to be confused with the active shear force. For details see ref. 8. Since the cilium is inextensible, we have ∂σ/∂s = ∂α/∂s, which implies σ(s, t) = α(s, t) + [σ(0, t) − α(0, t)]. It follows that during the effective stroke σ corresponds directly to the inclination angle of the cilium. Therefore, the criterion we choose for the switching between the effective and the recovery strokes is reaching a maximal inclination angle (denoted α_{right}) at the basal region.
During the recovery stroke, the maximal effective shear is attained when the cilium completes its straightening process. The criterion we choose for the switching between the recovery and the effective strokes is, therefore, returning to an erect position.
It is important to note that Eq. 5 represents a load (geometry)dependent engine, since the beat pattern and the frequency depend on the environment. External flow (e.g., due to neighboring cilia) and external load (e.g., changes in the viscosity of the surrounding fluid) affect α(s, t). Because the engine depends on the geometry and the switching depends on reaching some “final position,” changes in α(s, t) affect the beat frequency as well.
The parameter values used in our model are ω_{eff} = 11,000°/s, ω_{rec} = 2,290°/s, α_{left} = 130° and α_{right} = 20°. These are measured directly from data (beat cycles from ref. 10). During the effective stroke A_{1} = 0.26 and A_{2} = −0.17, and during the recovery stroke A_{1} = 2 for 0 ≤ s ≤ 0.1L and A_{1} = 1 for 0.1L ≤ s ≤ L, A_{2} = −2, and B_{eff} = 0 and B_{rec} = 2. These parameters were obtained by fitting the measured data.
Results
Fig. 1a displays snapshots of the beat cycle of a single cilium, having an engine described by Eq. 5 and dynamics computed by Eqs. 1–3. Clearly, the model produces a realistic beat pattern: during the effective stroke (dashed lines), the cilium moves almost as a straight rod, and during the recovery stroke (solid lines) a bend is generated at the basal end and propagates toward the distal end. Table 1 compares the model results to experimental data, demonstrating that the model cilium captures all essential features of the motion.
Machemer (1) observed that the beat frequency of the cilia of Paramecium decreases approximately linearly when plotted against the logarithm of viscosity. Brokaw (12) investigated this dependence with various types of flagella and found the same behavior. In the present model, both beat pattern and frequency change with increasing viscous load. In Fig. 2 we plot the calculated beat frequency as a function of the viscosity in the range μ_{water} ≤ μ ≤ 5μ_{water}, reproducing a linear dependence. Since this dependence is not explicitly assumed in the model, this suggests that the fundamental characteristics of the ciliary motility are captured by the present model. The ciliary beat frequencies (shown in Fig. 2) measured by Machemer in the range μ_{water} ≤ μ ≤ 5μ_{water} for a multicilia configuration in Paramecium are higher than what we obtain for the single isolated model cilium, but the slopes of the resulting lines are roughly similar. This difference is consistent with our results for multicilia configurations. As reported below, the beat frequency of the model cilia indeed increases in the presence of neighboring cilia. Note that since the cilia of Paramecium change the plane of beating when the viscosity is increased by more than about 5fold that of water (see ref. 1) and the present model does not include this possibility, we do not study the behavior of the model cilia outside this range.
To check the sensitivity of our model engine to external load, we subjected the model cilia to an external periodic flow, parallel to the cell’s surface and satisfying the noslip boundary condition. In the range from 29 to 55 Hz the cilium tends to match its beat frequency to the frequency of the external flow (data not shown).
Cilia interactions is another important aspect we study. When neighboring cilia beat close enough to each other, each one is subjected to the flow generated by the others. This influences their beat pattern and beat frequency. The model results indicate that the influence is negligible if the distance between the cilia exceeds two cilium lengths.
The simplest/smallest multicilia configuration consists of two cilia. Gray (13) reported on experiments with two flagella, beating initially at different frequencies and having different wavelengths. When these are put close enough to each other, they tend to quickly synchronize, and beat with the same frequency and in phase. An analogous phenomenon of selfsynchronization and phase locking between rings and lines of coupled oscillators is well known. It was investigated experimentally and theoretically (e.g., for onedimensional weakly coupled chains and rings of oscillators, see refs. 14 and 15), but here we have a complex realistic biological system. Fig. 3 displays snapshots of the beat cycles of two cilia separated by one cilium length, beginning their beat at opposite phases. As shown in the figure, they synchronize completely within two cycles. The final mutual beat frequency increases slightly from 29.5 Hz (for a single cilium) to 31 Hz. This is an indication of the “advantage” of some multicilia configurations. The presence of neighboring cilia (one in this case) “helps” them beat faster.
Metachronal coordination is characteristic of ciliary arrays. The cilia on same line (perpendicular to the direction of the effective stroke) are synchronized, and adjacent rows of cilia (parallel to the direction of the effective stroke) beat with a phase difference. Antiplectic metachronism is a situation in which the metachronal wave (defined by the tips of the beating cilia) propagates in a direction opposite to the direction of the effective stroke (e.g., in Paramecium, see refs. 1 and 10). In such configurations, each individual cilium is influenced by cilia located on the same line as well as by cilia located on neighboring rows. To ascertain whether the model would generate selforganized metachronism, we used multicilia configurations where the cilia start simultaneously with the same initial conditions. Fig. 4 displays snapshots of 10 cilia and a 100 cilia configurations. Both cases develop antiplectic metachronal waves, which is typical for the cilia of Paramecium. The beat frequency does not change significantly when the number of cilia exceeds 10 (42 Hz). This is because, as mentioned above, the spatial range of the interaction is limited. The cilia reach their steadystate beat patterns after a few cycles, and the results shown in Fig. 4 are already at steady state.
Liron (16) developed an efficient method for calculating the velocity induced at a point by an infinite line of Stokeslets, located at a fixed height above a flat surface (see also ref. 17). Changing the expression for g_{N} and g_{T} in Eqs. 1–3 to include the appropriate kernels (16) enables us to model the dynamics of an infinite line of synchronized cilia (the beat remains planar) and twodimensional arrays consisting of several rows each representing an infinite line of synchronized cilia. Fig. 1d shows the side view of an infinite line of synchronized cilia. It is interesting to observe that, although the beat pattern is changed compared with a single cilium (compare Fig. 1d with a), the resulting frequency (≈29 Hz) remains almost unchanged, unlike the case for a row of cilia. As for the singlerow configuration, metachronal phase differences emerge from the dynamics due to the hydrodynamical interactions.
This paper describes a framework for modeling the internal mechanism of cilia, which can be used for various cases based on experimental data. We used this approach successfully with data for Paramecium. Our results demonstrate that for closely packed ciliary arrays such as those found on ciliates, hydrodynamic interactions between neighboring cilia are sufficient to account for antiplectic metachrony. In principle, it is also possible to use our approach with appropriate data to investigate the formation of symplectic metachrony. Whether hydrodynamic effects alone are sufficient to cause diaplectic metachronal waves requires extension of the present model to include an “engine” that generates a realistic threedimensional beat. Such studies have been initiated. It should also be pointed out and that in airways, where the tips of the cilia penetrate into the lower surface of a thin layer of mucus during the effective stroke, hydrodynamic interactions between neighboring cilia may be insufficient by themselves to generate metachronal wave propagation. The properties of the mucus would need to be included in a suitable theoretical model.
Acknowledgments
This research was supported by the U.S.–Israel Binational Science Foundation (Grant 94242), by the Technion Vice President for Research Fund, and by the Fund for the Promotion of Research at the Technion. S.G. acknowledges the support of the Henri Gutwirth Fund for the Promotion of Research. K.L.G. acknowledges the support of the Technion Graduate School.
Appendix
The integral expressions for g_{N} and g_{T}, as derived by G.L. (5), are: 6 where G_{N} and G_{T} are the normal and tangential components of the vector G. 7 and 8 We use q = 1 μm (the parameter that determines the integration limits) and a = 0.1 μm (the radius of the cilium) as in ref. 5.
Here, U_{s}(r, r_{0}, φ) is the velocity induced at r by a Stokeslet with intensity φ, located at r_{0}. U_{si}(r, r_{0}, φ) is the velocity induced at r by a Stokeslet located at r_{0}, with intensity φ, and by its image system. U_{d}(r, r_{0}, φ) is the velocity induced at r by a doublet with intensity φ, located at r_{0}. U_{di}(r, r_{0}, φ) is the velocity induced at r by a doublet located at r_{0}, with intensity φ, and by its image system. V_{si}(r, r_{0}, φ) = U_{si}(r, r_{0}, φ) − U_{s}(r, r_{0}, φ). V_{di}(r, r_{0}, φ) = U_{di}(r, r_{0}, φ) − U_{d}(r, r_{0}, φ). Eq. 7 with the coefficients in 8 is an approximation up to the order O(), where ɛ = a/L measures the cilium’s slenderness (1/120 in our case). Details are given by G.L. (5).
Footnotes

↵† email: shay{at}math.technion.ac.il.

↵‡ email: konstant{at}math.technion.ac.il.

↵§ email: liron{at}math.technion.ac.il.

M. James Lighthill, University College London, London, United Kingdom

↵‖ The values used in the model are S_{0} = 10^{−12} N, L = 12 μm, E_{b} = 25⋅10^{−24} N⋅m^{2}, C_{N} = 0.0035 kg/m⋅s, and C_{T} = 0.0025 kg/m⋅s.
ABBREVIATION
 GL,
 Gueron and Liron
 Received December 11, 1996.
 Accepted March 24, 1997.
 Copyright © 1997, The National Academy of Sciences of the USA
References
 ↵
 Machemer H
 ↵

 Sleigh M A
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Sleigh M A
 ↵
 ↵
 Brokaw C J
 ↵
 Gray J
 ↵
 ↵
 ↵
 ↵