## 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

# Spontaneous oscillation and fluid–structure interaction of cilia

Contributed by Charles S. Peskin, March 13, 2018 (sent for review July 6, 2017; reviewed by Robert D. Guy and Philip C. Nelson)

## Significance

Motile cilia are hair-like structures that extend from the surfaces of eukaryotic cells and pump fluid. The motor protein dynein drives the motion of cilia, but it is unclear how the organized beat of a cilium emerges from the combined action of hundreds of dynein molecules. Here, we describe how the molecular architecture of a cilium constrains its possible motion, and we postulate a dynamical law that governs the dynein tensions. Based on this constraint and dynamical law, we construct a model cilium that beats spontaneously. Simulation of multiple copies of this model in fluid reveals that an array of cilia spontaneously breaks synchrony and adopts multiple phases and that such symmetry breaking is helpful in producing net fluid transport.

## Abstract

The exact mechanism to orchestrate the action of hundreds of dynein motor proteins to generate wave-like ciliary beating remains puzzling and has fascinated many scientists. We present a 3D model of a cilium and the simulation of its beating in a fluid environment. The model cilium obeys a simple geometric constraint that arises naturally from the microscopic structure of a real cilium. This constraint allows us to determine the whole 3D structure at any instant in terms of the configuration of a single space curve. The tensions of active links, which model the dynein motor proteins, follow a postulated dynamical law, and together with the passive elasticity of microtubules, this dynamical law is responsible for the ciliary motions. In particular, our postulated tension dynamics lead to the instability of a symmetrical steady state, in which the cilium is straight and its active links are under equal tensions. The result of this instability is a stable, wave-like, limit cycle oscillation. We have also investigated the fluid–structure interaction of cilia using the immersed boundary (IB) method. In this setting, we see not only coordination within a single cilium but also, coordinated motion, in which multiple cilia in an array organize their beating to pump fluid, in particular by breaking phase synchronization.

Motile cilia are hair-like organelles extending from the surfaces of eukaryotic cells and are characterized by their rhythmic wave-like motion. Most motile cilia have a “9 + 2” axoneme architecture, which refers to 9 peripheral doublet “AB” microtubules surrounding a central pair of singlet microtubules. Dynein motor proteins, responsible for force generation, extend from each A tubule of an outer doublet toward the B tubule of the adjacent doublet. Along each A tubule, the dynein motors are regularly spaced in groups with a period of 96 nm (1). The motor proteins generate relative sliding of the adjacent doublets by ATP-fueled cycles of attachment, translocation, and detachment along the B tubule of their target doublet. This active sliding interacts with the passive resistance from the other axonemal components, including the radial spokes and nexin links, and also with the basal body to produce the ciliary beat, yet it is still unclear how the organized motion of a cilium emerges spontaneously from the combined action of the hundreds of individual dynein motor proteins.

This question has been addressed in several different ways. The geometric factors that may regulate the dynein motor proteins include sliding displacement, curvature, and interdoublet distance as reviewed in ref. 2. In curvature-based models, the dynein motors are deactivated when the local curvature exceeds a specified threshold (3⇓–5). In these models, the direction in which the dynein motors move and the sign of curvature are opposite on two sides of the axoneme with respect to the beating plane. In sliding control models, linearity between the sliding displacements and the shear forces is assumed (2, 6⇓–8), and this has been combined with a two-state model of the dynein motor involving strongly and weakly bound states depending on the sliding displacement (8, 9). The geometric clutch model (10, 11) postulates that the interdoublet distance is responsible for the activation of the dynein motors. This distance depends on the transverse force, which develops across the structure in the beating plane when the structure is bent. The active bending and the passive bending of the cilium have an opposite effect on the transverse force, and this antagonism generates the ciliary beating. In another class of models, the cilium is treated as a slender body, and the dynein motors are not modeled in detail. Instead, an internal engine is postulated to generate the beat pattern of the cilium (12⇓–14).

Our model is based on a geometrical constraint that restricts the deformation of the cilium and on a postulated dynamical law for the tension generated by a dynein motor. The geometrical constraint allows us to retain important aspects of the microstructure of the cilium while reducing the dfs of the model to those of a single centerline space curve. The postulated dynamical law allows the tension generated by each dynein motor to evolve independently but in a manner that depends on the length of the link formed by that dynein motor and hence, on the deformation of the structure as a whole. The dynein motors in our model are, therefore, coupled to each other through the structure of the cilium, and this coupling is indirectly responsible for their coordination. The model cilium has an unstable symmetrical steady state, in which the cilium is straight and all motors are under equal tension. The instability of this steady state leads to the emergence of a globally stable limit cycle oscillation, which is the ciliary beat.

The interpretation of ciliary beating as a limit cycle oscillation around an unstable steady state was proposed in ref. 8. The authors presented a local drag model of a 2D cilium driven by two-state molecular motors (9). The instability of the steady state of the model cilium was shown to involve a Hopf bifurcation. Our model is 3D, and our postulated dynamical law governing the tension of an individual dynein motor is continuous rather than discrete.

In ciliary arrays, neighboring cilia interact with each other via surrounding fluid, and this leads to the emergence of phase synchronization and metachronal waves, in which a constant phase difference is maintained between adjacent cilia. The hydrodynamic interactions between cilia have been studied analytically and numerically with various models of cilia from minimal models capturing the essential features of ciliary beating (15⇓⇓⇓⇓–20) to more detailed models to simulate the beating shape of real cilia (12, 14, 21⇓⇓–24). In particular, metachronal coordination of arrays of cilia is believed to be connected to efficient pumping of fluid, which is the principal function of motile cilia. The efficiency of pumping fluid in a ciliary array has been studied by using magnetically actuated artificial cilia (25, 26), a model of cilia with active bending forces (23), or a model with an optimized beat shape for pumping efficiency (21). These studies address the pumping efficiency in relation to phase difference, wave vector, and interciliary spacing. In both studies of 1D (25, 26) and 2D ciliary arrays (21, 23), the metachronal wave pumps fluid more effectively than synchronized ciliary beating, and simplectic metachronal waves, which propagate in the same direction as the power stroke, are generally more efficient in fluid transport than antiplectic waves, which propagate in the opposite direction.

In contrast to previous works focusing on the emergence in phase difference in orthoplectic rows [running in the same or opposite direction as the ciliaty beat (27)], here, we present simulation results showing emergence of multiple phases (i.e., symmetry breaking) in a diaplectic row (which is perpendicular to the beating plane of the cilia) (27). In this simulation, we find that symmetry breaking is essential for effective fluid transport.

## Model

According to the 9 + 2 arrangement of the microtubules, we model the configuration of a cilium as 11 discretized space curves in

### Geometric Constraint.

We assume that all of the microtubules are inextensible, so that the distance between adjacent nodes along each microtubule remains constant. That is,**1**] and [**2**] are combined, however, it can be shown (*Supporting Information*) that the configuration of each cross-section is a rigid translation of the base cross-section of the model cilium. We assume, moreover, that the base cross-section of the cilium is planar and has the form of a fixed, regular nonagon. That is, for all **3** allows neither twist nor tilt of the different cross-sections of the cilium and implies that the configuration of the entire cilium is completely determined by that of its centerline. This enables us to retain the detailed structure of the cilium in our model and yet to describe its motion more simply in terms of the motion of the centerline. Since the base cross-section is a planar and regular nonagon, the whole structure can be thought of as a freely jointed chain of regular nonagons. Although the centerline is freely jointed, the nonagons themselves are constrained to move rigidly and by translation only.

In our model, the relative sliding of adjacent peripheral doublet microtubules happens automatically as the local orientation of the centerline changes. This consequence of our geometric constraints is derived in *Supporting Information*. It is shown there that, as the direction of the centerline changes, its rate of change induces a proportional relative sliding velocity of adjacent peripheral curves of the model. The relative sliding velocity varies sinusoidally around each cross-section, being zero in the direction in which the tangent to the centerline is changing and having the opposite sign on either side of that direction. All of these quantities and therefore, the relative sliding velocities can be different in the different cross-sections of the model cilium.

### Forces.

The motion of the cilium is driven by the active dynein motors coupled to the passive elasticity of the structure. We model the dynein motor proteins by active links connected between two points of the adjacent peripheral curves, *Supporting Information*.

The active links in our model generate tension governed by the following differential equation:*Supporting Information*). The level dependence of **6** and that the only information available to any active link about the state of the cilium as a whole at any given time is the instantaneous length of that particular link, **6** comes from a toy model, described in *Supporting Information*, in which a rigid rod hinged at its base is held in the vertical position by two guy wires under tension. The stability of the vertical equilibrium depends on the tension law of the guy wires. If the wires are under constant tension, the vertical equilibrium is unstable; if the tension is proportional to the length of the wire, the vertical equilibrium is neutrally stable, and if the tension is proportional to the square of the length or indeed, to any power of the length greater than one, the vertical equilibrium is stable. Our postulated dynamical law for the tension, [**6**], contains aspects of the first and third cases. On a fast timescale, Eq. **6** simply asserts that the tension is prescribed (by its instantaneous value, since Eq. **6** requires **S32**–S36, but we choose the power two for simplicity.) Indeed, when this equation is used to govern the guy wires of the toy model and when the parameters of that model are appropriately chosen, the result is a limit cycle oscillation about the vertical configuration.

Although the detailed mechanism by which dynein generates tension is unknown, the prevailing view (29, 31) is that this process involves cycles of attachment and detachment linked to ATP hydrolysis, like the cross-bridge cycle in muscle (32). We do not attempt to model the details of this process, however, and instead, we consider only the resulting tension, which we model in a phenomenological way, although we do give a possible mechanistic interpretation of Eq. **6** in *Supporting Information*. In particular, in our model, the whole group of dynein motors within a given 96-nm unit cell of an A tubule is represented by a single link that is permanently attached to the next B tubule in the next lower cross-section of the model (Fig. 1*A*).

### Cilia in Fluid.

We model the fluid–structure interaction of cilia by the immersed boundary (IB) method (33), with the following equations of motion:**7** and **8** are the incompressible time-dependent Stokes equations with fluid velocity *Supporting Information*. The term **S68**–S71). Inclusion of this term does not adversely affect the time step restriction for numerical stability, because our implementation is implicit (34) as described in Eqs. **S55**–S67. Eqs. **11** and **12** define the elastic energy of the model cilia. In these equations, **11**, **5** and also, the elastic energy of the 5-6 bridges), and the sum over *Supporting Information*, Eq. **6** is given a mechanistic interpretation in terms of the making and breaking of constant tension links. Eqs. **9** and **10** are interaction equations that describe the influence of the cilia on the fluid and that of the fluid on the cilia. In Eq. **9**, **S51**–S54). It defines the way in which the force is distributed over the neighborhood of the point **10**, the velocity of centerline node **9**. This ensures that

## Results

In this section, we report the results of two computer experiments. The first one involves a single model cilium immersed in fluid, and in the second one, we consider the fluid-mediated interaction of multiple cilia.

### Single Cilium in Fluid.

We put the cilium in a rectangular box of fluid with periodic boundary conditions, but the periodicity in the

In this study, we compare two different initial conditions. In both of them, the initial configuration of the cilium is straight and vertical, but in the first case, the initial tensions of the dynein links are all zero, whereas in the second case, the initial tensions are randomly chosen (Fig. 1 *C1* and *D1*).

Despite this difference in initial conditions, the long-term behavior is the same in both cases (except for phase). Both cilia settle into the same limit cycle oscillation, as shown in Fig. 1 *C3* and *D3*. This limit cycle shows a realistic ciliary beat with clearly distinct power and recovery strokes (Fig. 1*B*). At early times, however, we see a clear and qualitative distinction between the behaviors of the two model cilia. In the first case, in which the initial tensions are all zero, the system has to break symmetry to move, and this takes some time to occur. Thus, the tensions build up gradually while remaining nearly equal to each other within any given level of the cilium. This is shown in Fig. 1*C3* (up to *D3*), there is no need to break symmetry, because the initial tensions are random, and the cilium quickly and smoothly settles into its limit cycle oscillation.

The comparison of these two cases illustrates clearly the self-organizing character of the beat of the model cilium, in which hundreds of dynein motors coordinate successfully without any explicit coordination mechanism.

### Multiple Cilia in Fluid.

The computer experiment described in this subsection is similar to those of the previous subsection, but the dimensions of the fluid box are different (

The initial conditions here are the same as those of the first case above: zero initial dynein tension with the cilia straight and vertical. Thus, we expect and find an initial pause after which the cilia start to beat, and we observe that their early beating is synchronized. After a few beats, however, the synchrony spontaneously breaks, and the different cilia adopt different phases.

This symmetry breaking turns out to have an important effect on the ability of the ciliary array to pump fluid. Before the symmetry breaking, the cilia are unable to exhibit a clear distinction between their power and recovery strokes (Fig. 2*B*). To understand this, note that the asymmetrical beating pattern that pumps fluid requires the tips of the cilia to move down toward the cell boundary to which the cilia are attached (which is a no-slip planar boundary in our simulation) during part of the beat cycle, so that the tips can be closer to the cell boundary during the recovery stroke than during the power stroke. This downward motion is inhibited by fluid incompressibility and viscosity when the cilia are close together and synchronized. Symmetry breaking alleviates this effect.

## Summary and Conclusions

We have introduced a phenomenological model of a motile cilium, in which ciliary beating is driven by actively generated dynein tensions governed by a postulated dynamical law. Under two simple geometric constraints motivated by the inextensibility of microtubules and the centering of the central pair as maintained by the radial spokes of the cilium, we show that the entire configuration of the cilium can be expressed in terms of the configuration of its centerline.

A significant feature of our model is that the hundreds of dynein motors within a single cilium coordinate their tensions spontaneously to produce the ciliary beat as an emergent limit cycle oscillation. To emphasize this, we have shown that the same ciliary beat emerges regardless of initial conditions, in particular starting from zero initial tensions or from random initial tensions of the dynein motors. In the former case, there is a delay before beating starts, during which the dynein links of the model develop tension synchronously and the cilium does not move. The symmetry then breaks suddenly, and the beat is initiated. In the latter case, with random initial tensions, there is no need for symmetry breaking and hence, no delay before the cilium starts to beat.

We have also studied the fluid–structure interaction of multiple cilia by the IB method. In particular, we have considered the interactions that occur within a diaplectic row. In this setting, we see spontaneous symmetry breaking, in which the cilia adopt multiple phases, despite being started in the same initial state. Only after the symmetry breaks are the cilia able to pump fluid.

## Acknowledgments

We thank David McQueen for helpful advice concerning the computations in this paper and for providing the software that was used to visualize the results. J.H. was supported by the McCracken Fellowship and the Dean’s Dissertation Fellowship of the New York University Graduate School of Arts and Science.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: jihunhan{at}cims.nyu.edu or peskin{at}cims.nyu.edu.

Author contributions: J.H. and C.S.P. designed research, performed research, analyzed data, and wrote the paper.

Reviewers: R.D.G., University of California, Davis; and P.C.N., University of Pennsylvania.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1712042115/-/DCSupplemental.

Published under the PNAS license.

## References

- ↵
- ↵.
- Sartori P,
- Geyer VF,
- Scholich A,
- Jülicher F,
- Howard J

- ↵.
- Dillon RH,
- Fauci LJ,
- Omoto C

- ↵.
- Brokaw CJ

- ↵.
- Brokaw CJ

- ↵
- ↵
- ↵.
- Camalet S,
- Jülicher F

- ↵
- ↵
- ↵
- ↵.
- Gueron S,
- Levit-Gurevich K,
- Liron N,
- Blum JJ

- ↵
- ↵.
- Gueron S,
- Levit-Gurevich K

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Osterman N,
- Vilfan A

- ↵
- ↵.
- Elgeti J,
- Gompper G

- ↵
- ↵.
- Khaderi SN,
- den Toonder JMJ,
- Onck PR

- ↵
- ↵.
- Knight-Jones EW

- ↵.
- Afzelius B

- ↵.
- Lindemann CB,
- Lesich KA

- ↵.
- Han J

- ↵.
- Summers KE,
- Gibbons IR

- ↵
- ↵
- ↵.
- Mori Y,
- Peskin CS

- ↵

## Citation Manager Formats

### More Articles of This Classification

### Biological Sciences

### Biophysics and Computational Biology

### Related Content

- No related articles found.