## Abstract

Cilia and flagella are highly conserved organelles that beat rhythmically with propulsive, oscillatory waveforms. The mechanism that produces these autonomous oscillations remains a mystery. It is widely believed that dynein activity must be dynamically regulated (switched on and off, or modulated) on opposite sides of the axoneme to produce oscillations. A variety of regulation mechanisms have been proposed based on feedback from mechanical deformation to dynein force. In this paper, we show that a much simpler interaction between dynein and the passive components of the axoneme can produce coordinated, propulsive oscillations. Steady, distributed axial forces, acting in opposite directions on coupled beams in viscous fluid, lead to dynamic structural instability and oscillatory, wave-like motion. This ‘flutter’ instability is a dynamic analogue to the well-known static instability, buckling. Flutter also occurs in slender beams subjected to tangential axial loads, in aircraft wings exposed to steady air flow and in flexible pipes conveying fluid. By analysis of the flagellar equations of motion and simulation of structural models of flagella, we demonstrate that dynein does not need to switch direction or inactivate to produce autonomous, propulsive oscillations, but must simply pull steadily above a critical threshold force.

## 1. Introduction

Flagella and cilia are slender organelles that generate propulsive, oscillatory waveforms (figure 1*a*) to propel cells or move fluid. The common cytoskeletal structure of cilia and flagella, the ‘9 + 2’ axoneme, consists of nine outer doublet microtubules, a central pair, radial spokes and circumferential links (figure 1*b*). Dynein molecules form an array of cross-bridges between pairs of microtubule doublets and exert forces that cause sliding of one doublet relative to the other. These active shear forces combine with the forces from passive structural elements (doublets, nexin links and radial spokes) to produce bending [1].

Since dynein activity on opposite sides of the axoneme produces antagonistic bending forces, it is widely believed that dynein activity must be regulated (switched on and off, or modulated) in order to produce oscillatory waveforms [2–6]. This concept has been supported by a number of theoretical studies that explored the effect of feedback from mechanical deformation to dynein force. In a seminal paper [7], Brokaw showed that delayed feedback from curvature to dynein activity could lead to oscillations. Brokaw has continued to describe and evaluate potential feedback mechanisms in a comprehensive series of modelling papers [8–12]. Hines & Blum [13] developed a detailed mathematical model of flagellar motion in the form of partial differential equations (PDEs); numerical solutions to these equations exhibited sustained, wave-like oscillations when feedback from curvature was included. Murase [14] and Julicher and colleagues [5,15,16] have modelled sliding-controlled, collective dynein behaviour and demonstrated the existence of oscillatory modes in such models. In the ‘geometric clutch’ (GC) hypothesis of dynein regulation proposed by Lindemann [17–19], the spacing between each doublet pair is thought to control the level of dynein cross-linking. In the GC model, dynein activity increases when doublets approach one another and decreases when doublets separate.

Though these feedback models are appealing, experimental evidence that dynein regulation is required for bending oscillations remains circumstantial. For example, Lin *et al.* [20] has shown by cryo-electron microscopy (cryo-EM) that the dominant conformation of dynein arms (pre-powerstroke or post-powerstroke) differs between adjacent doublet pairs in bent sections of sea urchin flagella. Movassagh *et al.* [21] also observed that dynein activity was organized into clusters, though no clear relationship was seen between the observed patterns of activity and bending. By analysing waveforms of human sperm flagella, Gaffney *et al.* [22] and Smith *et al.* [23] obtained estimates of oscillatory dynein force from estimated viscous and elastic forces. On the other hand, Brokaw [3] has commented that the belief that dynein activity switches ‘on’ and ‘off’ derives mainly from the lack of alternative explanations for oscillation. Intriguingly, Geyer *et al.* [24] have shown show that gradually increasing dynein activity in reactivated axonemes leads both to (i) a steady increase in static mean curvature and (ii) an abrupt transition to oscillatory beating. These data suggest that simply increasing average dynein activity can lead to instability of the static equilibrium state of the axoneme.

In this study, we show that switching or modulation of dynein activity is not required to generate propulsive, oscillatory waveforms. In an elastic structure submerged in viscous fluid, like the axoneme, steady, distributed axial loading of the doublets (i.e. steady dynein activity) leads to a dynamic structural instability commonly known as flutter. Flutter is a well-known cause of oscillations in aircraft wings and panels exposed to steady flow [25,26], axially loaded beams [27–29] (figure 2*b,c*), and flexible pipes conveying fluid [32] (figure 2*d*). This dynamic instability can occur in flagella instead of the familiar static instability, buckling, because dynein forces remain tangent to the doublets as they deform. Oscillations arise because deformation of the structure itself changes the direction of these tangential loads.

The idea that steady dynein activity is sufficient to cause flagellar beating is supported by results from three different, but complementary, mathematical methods: (i) stability analysis of the linearized partial differential equations (PDEs) of doublet motion, using representative, realistic values for biophysical parameters including doublet tension, (ii) simulation of nonlinear equations of doublet motion, generalized from prior work [13,33], and (iii) simulation of three-dimensional (3D) finite-element (FE) models of flagella subjected to steady, distributed, axial inter-doublet forces. Each model and analysis, described in the sections below, predicts oscillatory, propulsive waves propagating from the base to tip of the flagellum.

## 2. Modelling and analysis of flagella

### 2.1. Axoneme structure and loading

Flagellar bending is caused by the activity of dynein arms between the outer microtubule doublets (figure 3; electronic supplementary material, figures S1 and S2). Dynein arms are permanently attached to the A-microtubule of doublet *N*; these arms transiently attach and ‘walk’ towards the base of the B-microtubule of doublet *N* + 1, exerting a tipward force on the higher numbered doublet. In combination with passive components, such as radial spokes, and nexin links, dynein forces also produce bending moments. If dynein activity occurs only between doublets on one side the entire flagellum will bend in response (figure 3*a*). The side that produces the ‘principal’, or *P* bend, is denoted the *P* side; dynein activity on the other (*R*) side produces the reverse (*R*) bend. If equal dynein activity occurs between doublet pairs on opposite sides of the flagellum (figure 3*b*), the net bending moment will be zero, and the equilibrium configuration of the flagellum will be straight (and slightly twisted). However, the equilibrium may not be stable.

When dynein activity is not uniform around the circumference (as in figure 3*a,b*) some doublets may be under significant tensile or compressive loads, with respect to critical values for beam instability (buckling or flutter). In flagella of *Chlamydomonas reinhardtii*, there are 19 dynein heavy chains per 96 nm repeat between each doublet pair, or over 195 arms µm^{−1}. Each arm may exert at least 2–4 pN [34,35] average force, so that coordinated dynein activity could lead to axial forces well over 1000 pN at the base of a flagellum. For comparison, the critical buckling load for a 12 µm long microtubule doublet with (order of magnitude) flexural rigidity *EI* = 100 pN µm^{2} is . Why do not axonemal doublets fail by buckling? One reason is that most doublets experience dynein forces in both directions (though the net dynein force will likely be non-zero). In addition, doublets under compression are connected by other axoneme components to doublets that are under tension. Nonetheless, axial loads are clearly significant, though they have not been considered in most models, the notable exception being the GC hypothesis [17–19].

To gain insight into the mechanism of instability, we first represent the 3D mechanics of the axoneme by a two-dimensional (2D) model (figure 3*c,d*) and a corresponding set of one-dimensional (1D) PDEs. Although 3D structural models are straightforward to simulate using FE methods, analysis of 1D PDEs is much more efficient and can provide greater insight into mechanisms. Figure 3*c* illustrates the accumulation of axial stress due to distributed dynein forces in a 2D, two-doublet model. Figure 3*d* shows the scenario in which the net bending moments are zero (counteracted by the opposing doublet pair), but the net axial forces are significant. The equations of motion for this situation were derived (electronic supplementary material) and analysed to determine the existence and nature of structural instability. To confirm and supplement the analysis of these 1D PDEs, the axoneme was also modelled as a 3D structure, using commercial FE software (COMSOL v. 5.1, COMSOL, Inc., Burlington, MA, USA), and stability and oscillations were investigated in this full 3D structural FE model.

### 2.2. Equations of motion: one-dimensional PDE model

#### 2.2.1. Nonlinear equations of motion for two coupled doublets

In the axoneme, doublets under compression are coupled to doublets in tension by the other components of the axoneme (dynein, radial spokes and nexin–dynein regulatory complexes). Equations governing the entire flagellum have been derived and analysed previously [5,13,15,16,33,36]. Using a similar approach to model individual doublet behaviour in a two-doublet model (electronic supplementary material, section A and figure S2), we obtain three equations for each doublet, describing the motion of a slender elastic beam moving in viscous fluid, subject to active and passive inter-doublet forces [13,15]:
2.1a
2.1b
2.1c
2.1d
2.1e
2.1fIn these equations . Key variables of the model are summarized in table 1. Parameters for *Chlamydomonas* flagella [37,38] (table 2) are used, unless otherwise noted. Most parameters have been measured experimentally (see references in table 2), but several are simply plausible estimates (denoted by superscript ‘a’ in table 2). In addition, boundary conditions (electronic supplementary material, section A and equations A.9(i–viii)) are needed to specify that each doublet is fixed at the base and free at its tip: zero angle, displacement and velocity at the base, and zero external moments and forces at the tip.

#### 2.2.2. Inter-doublet forces and moments

Passive components of the axoneme resist doublet separation (electronic supplementary material, figure S2), coupling the motion of the two doublets. These components are modelled as distributed elements with linear elastic and viscous coefficients, *k*_{N} and *b*_{N}, respectively, and a nonlinear elastic coefficient .
2.2In the tangential direction we include a *steady*, distributed dynein force per unit length, , as well as passive resistance to inter-doublet sliding. The dynein distribution function, *D*(*s*), allows the steady dynein activity to vary along the length. For longitudinally uniform dynein activity, *D*(*s*) = 1; to represent dynein activity that increases distally, . The moments due to equal, opposing dynein forces on the *P* and *R* sides cancel (electronic supplementary material, figure S1; figure 3*b*), but the moments due to passive shear resistance sum. The net inter-doublet moment per unit length is thus
2.3where the mean shear angle is . The coefficients *k*_{T} and *b*_{T} represent linear elastic and viscous resistance to sliding [5,33,40], and a nonlinear stiffness. The nonlinear terms in equations (2.2) and (2.3) represent assumed nonlinear, stiffening behaviour of elastic elements between doublets [8].

### 2.3. Stability analysis of the linearized one-dimensional PDE model

To describe small-amplitude motion about a straight, equilibrium configuration, equations (2.1*a*–*f*) may be linearized to obtain much simpler equations [5,15,42]. For a uniform, steady dynein force distribution , dropping terms that are nonlinear in the dependent variables, the baseline tension (or compression) in the two doublets is (electronic supplementary material, figure S3):

Doublet 1: (tension)

Doublet 2: (compression)

Equations (2.1*a*–*f*) become
2.4a
2.4b
2.4cwith corresponding fixed-free boundary conditions. Analogous linearized equations can be derived for distally increasing dynein activity (electronic supplementary material, figure S3).

Following [5,15,16] and [33,36] separable solutions of the form
2.5are sought, with (*α* and *ω* are real). If *α* > 0, the amplitude of the corresponding mode grows exponentially. These separable solutions are substituted into the linearized equations of motion (equations (2.4*a*–*c*)). After defining a characteristic time for the system, , and a normalized eigenvalue, , the resulting ordinary differential equations (ODEs) may be written in non-dimensional form:
2.6a
2.6bwhere the new non-dimensional parameters are , and .

These two coupled ODEs, together with the associated fixed-free boundary conditions, form an eigenvalue problem. Approximate solutions to the eigenvalue problem can be obtained using numerical methods such as FE analysis or the method of weighted residuals [36,43] (electronic supplementary material, section B).

### 2.4. Simulation of the nonlinear one-dimensional PDE model

Time-domain numerical simulations of the full nonlinear PDEs [37,38] for the two-doublet model (equations (2.1*a*–*f*)), with fixed-free boundary conditions, were performed using a commercial software package (COMSOL v. 5.1; COMSOL, Inc. ). Time-domain simulation can confirm the contributions of the unstable modes predicted by stability analysis, and also allows exploration of nonlinear, transient and non-periodic behaviour. Simulations were performed with the same baseline parameters used in the stability analysis, with various values of the steady, distributed dynein load *p*. The 1D domain was discretized into 201 elements with quartic polynomial interpolation, and a backward differentiation algorithm with variable time step and relative tolerance of 1 × 10^{−8} was used to solve the system. Representative results were checked with finer mesh and smaller relative tolerance.

### 2.5. Simulation of three-dimensional structural finite-element models

Three-dimensional FE structural models were constructed using the same commercial FE package (COMSOL v. 5.1, COMSOL, Inc.). The basic model consisted of six slender beams (outer doublets) in a hexagonal array connected to a central beam (central pair) by radially oriented beams (spokes). Radial beam elements coupled the doublets to each other and to the central pair, providing both normal and shear stiffness. Normal and tangential viscous force components, proportional to the normal and tangential velocity components, respectively, of each doublet, represented the effect of fluid at low Reynolds number. Dynein activity was modelled by steady, uniformly distributed axial forces tangent to each doublet, and corresponding steady, distributed bending moments to satisfy equilibrium conditions (electronic supplementary material, section A and figure S2).

The 3D FE model was composed entirely of linearly elastic, beam elements. Different values of flexural rigidity (*EI*) and axial stiffness (*EA*) were used for doublets (longitudinal members) and radial spokes. In the spokes, baseline values and were multiplied by quadratic functions of curvature or strain, to model assumed nonlinear stiffening (table 3). These nonlinear terms are chosen for convenience, to represent in simple form the behaviour of elastic components with limited range. Inertial effects were included, but inertial forces are extremely small relative to elastic and viscous forces (electronic supplementary material, section C and figure S8). Mass scaling [46] was used to speed computation; the ratio of kinetic energy to potential energy was checked to confirm that inertial effects could be neglected. The system was discretized into 260 beam elements, with cubic polynomial interpolation of the displacement field. Time-marching simulations (backward differentiation, variable time step, relative tolerance 1×10^{−6}) were performed. Representative results were checked with finer mesh and relative tolerance.

Parameters of the 3D FE structural model (table 3) were chosen to approximate the overall parameters of *Chlamydomonas reinhardtii* flagella [37,38]. Each of the seven individual doublets of the 3D model is more flexible (*EI* = 100 pN µm^{2}) than each of two doublets of the 1D PDE model (*EI* = 400 pN µm^{2}). As a result, critical dynein force densities differed between the two models by roughly a factor of 4.

## 3. Results: dynamic instability and oscillations

### 3.1. Instability and oscillatory modes of the one-dimensional PDEs

Modes and frequencies of oscillation were obtained by analysis of the linearized PDEs of two coupled doublets (equation (2.4*a*,*b*)) subject to steady, uniformly distributed dynein loading, using the method of weighted residuals (electronic supplementary material, section B). The eigenvalues of the system depend on the amplitude, *p*, of the distributed dynein force (figure 4*a,b*). Loss of stability is indicated by a positive real part of any eigenvalue; the finite imaginary part of that eigenvalue confirms the oscillatory nature of the instability, and specifies its frequency. The eigenfunction corresponding to each unstable eigenvalue (figure 4*c,d*; electronic supplementary material, movie M1) describes the shape of the growing, oscillatory waveform. All eigenfunctions of this model exhibited anterograde (base-to-tip) propagation.

The stability of the system can be established for multiple parameter combinations, as the stability analysis algorithm is quite efficient. As expected, increasing dynein force density and length are both associated with loss of stability (figure 4*e*). Increasing inter-doublet normal stiffness appears to increase the threshold for onset of oscillations (figure 4*f*). Neither inter-doublet viscous resistance nor shear stiffness had a significant effect on stability thresholds (electronic supplementary material, figure S4). Increasing dynein force density increases the frequency of oscillation (figure 4; electronic supplementary material, figure S4), whereas increasing the inter-doublet stiffness or viscosity tends to decrease frequency.

### 3.2. Propagating waves in time-domain simulations of the one-dimensional PDE model

Time-domain simulation of the system of 1D PDEs revealed oscillatory solutions (figure 5) consistent with modes predicted by stability analysis. Larger amplitude waves were obtained using the distally increasing, steady dynein force distribution. Near the stability threshold (*p* = 450 pN µm^{−1}), the solution consisted of small-amplitude waves that closely resembled the single unstable mode of the linearized model. At higher values of the steady dynein load (*p* = 600 and 800 pN µm^{−1}), the contributions of more than one unstable mode are apparent. All modes propagate from base to tip.

Above a critical value of axial force density (figure 6), the 3D structural FE model also exhibited oscillatory, wave-like behaviour propagating from base (fixed end) to tip (free end; figure 6*a,b*). The qualitative behaviour of the system, including onset of oscillation, and relative deformations of doublets in compression and tension are consistent with the predictions of the 1D PDE models. The threshold for onset of oscillations is lower in the 3D FE model because the flexural modulus of each doublet (100 pN µm^{2}) in the 3D FE model is lower than the flexural rigidity of each doublet (400 pN µm^{2}) in the two-doublet 1D PDE system. The 3D FE model also did not include a viscous inter-doublet component, so inter-doublet time constants are lower and oscillation frequencies are higher than those of the 1D PDE model.

## 4. Discussion and conclusion

Steady tangential dynein forces, combined with fluid–structural interactions, can initiate and maintain propulsive flagellar waveforms, without feedback to or regulation of dynein activity. Evidence for this mechanism is provided by (i) eigenanalysis of a simple linearized model, (ii) simulation of a set of nonlinear PDEs, and (iii) simulation of an FE structural model. The possibility that propulsive oscillations arise from a mechanical instability is distinct from prior models based on switching or modulation of dynein. These previous models include curvature-controlled dynein dynamics [7,13], and dynein regulation by inter-doublet forces (the GC model), which both produce propulsive, anterograde (base-to-tip) propagating waveforms [36]. Flagella models with dynein controlled by inter-doublet sliding can also produce oscillatory waveforms [5,15,16,36], though corresponding waveforms appear to be mostly retrograde. The possible role of nonlinear buckling instability in flagella was explored by Gadelha [47] in a model in which oscillations were driven by imposed, periodically varying dynein activity.

By contrast, the mechanism of fluid–structural flutter proposed here does not require dynein regulation or switching. Instead, the fluid–structural instability and resulting oscillations are driven by axial tension in individual doublets. The importance of doublet tensile and compressive loading was first raised by Lindemann [17–19] in the development of his (GC) hypothesis. However, as noted above, the GC hypothesis involves explicit mechanical feedback from transverse force to dynein activity. In the current model, dynein activity varies circumferentially (i.e. differs between doublets) and longitudinally, but does not vary with time.

Unlike classical systems that exhibit dynamic instability, inertial effects are not necessary for this instability and in fact inertia is negligible at the low Reynolds numbers involved in flagella motion. Viscosity, however, can help destabilize the flagellum by introducing phase shifts between doublet motion and force on the doublets, allowing work to be done on the flagellum during each cycle [32]. This instability occurs over a wide range of parameter values (figure 4*e,f*; electronic supplementary material, figure S4*a,b*).

The possibility of structural flutter is consistent with a number of prior experimental observations. If this model is correct, increasing steady dynein activity explains both the gradual increase in static curvature and abrupt onset of oscillations, which have been observed in reactivated axonemes as adenosine triphosphate concentration is increased [24]. The current model is also consistent with observations of intermittent buckling and re-attachment observed between doublets from which inter-doublet links are absent [48]. The effect of increased fluid viscosity on beat frequency [49–51] is captured naturally because oscillation frequency scales inversely with the major viscoelastic time constant of the system: . The decrease in beat frequency of outer dynein arm-deficient mutants (e.g. *oda2*) is consistent with the general decrease in frequency with force density (figure 4*e,f*). The paralysis of *pf13* mutants missing outer dynein arms and one species of inner dynein arm is consistent with sub-threshold dynein activity. The direction and timing of oscillations arising from fluid–structural flutter instability can be influenced by subtle cues, such as the orientation of the central pair [52], structural asymmetry of the axoneme [53] or hydrodynamic effects from neighbouring flagella [54].

The models of flagella mechanics used in this study are simplified in order to isolate and emphasize the basic mechanism of instability. These models have inevitable limitations. The effects of biased dynein activity, or inherent curvature, on flagellar stability and oscillations have not been explored. As a result, the predicted waveforms are symmetric. The interactions between doublets in the current model are not satisfactory for describing large amplitude behaviour; inter-doublet components are undoubtedly more complex than the simple spring–dashpot system used here. Neither the possible collective behaviour of dynein arms, nor their force–velocity relationships are modelled. The 3D structural FE models are composed completely of linear, elastic beam elements. While geometric nonlinearity due to large deformations is considered, possible material nonlinearities, sliding or loss of contact between components are not accounted for. A notable structural feature omitted from the current model is the twisted central pair of microtubules, which may enhance twisting or out-of-plane motion. Twisting is present, but minimal, in the current 3D FE models, but may be important in flagella beating [55].

Despite these caveats, this study shows that dynamic structural instability described here is *sufficient* to produce propagating, oscillatory, waveforms. In other words, dynein regulation is *not necessary* for flagella beating, although mechanical feedback to dynein activity may still affect waveform frequency, amplitude or symmetry. Estimates of propagating, oscillatory, dynein forces in prior studies of human sperm [23] support a role for such feedback in sperm flagella. However, this evidence remains indirect, and dependent on precise knowledge of fluid resistive force and flagellum elastic parameters.

In summary, propulsive flagellar beating can arise from steady, distributed dynein forces acting between coupled, flexible doublets in viscous fluid. This mechanism, analogous to flutter in airframes and flexible pipes, is simple, robust to parameter variations and consistent with experimental observations. The possible role of this dynamic instability in flagellar oscillations is a compelling target for future experimental and theoretical investigations.

## Authors' contributions

P.V.B. and S.K.D. conceived the model combining insights from mechanics (P.V.B.) and biology (S.K.D.). P.V.B. derived equations, and performed analyses and simulations. Both authors helped draft the manuscript and gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

Support from NSF grant no. CMMI-1265447 and CMMI-1633971 (to P.V.B.), NIH grant no. GM032842 (to S.K.D.) and the Children's Discovery Institute is gratefully acknowledged.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3512442.

- Received June 30, 2016.
- Accepted September 26, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.