## Abstract

The mechanisms underlying the coordinated beating of cilia and flagella remain incompletely understood despite the fundamental importance of these organelles. The axoneme (the cytoskeletal structure of cilia and flagella) consists of microtubule doublets connected by passive and active elements. The motor protein dynein is known to drive active bending, but dynein activity must be regulated to generate oscillatory, propulsive waveforms. Mathematical models of flagellar motion generate quantitative predictions that can be analysed to test hypotheses concerning dynein regulation. One approach has been to seek periodic solutions to the linearized equations of motion. However, models may simultaneously exhibit both periodic and unstable modes. Here, we investigate the emergence and coexistence of unstable and periodic modes in three mathematical models of flagellar motion, each based on a different dynein regulation hypothesis: (i) sliding control; (ii) curvature control and (iii) control by interdoublet separation (the ‘geometric clutch’ (GC)). The unstable modes predicted by each model are used to critically evaluate the underlying hypothesis. In particular, models of flagella with ‘sliding-controlled’ dynein activity admit unstable modes with non-propulsive, retrograde (tip-to-base) propagation, sometimes at the same parameter values that lead to periodic, propulsive modes. In the presence of these retrograde unstable modes, stable or periodic modes have little influence. In contrast, unstable modes of the GC model exhibit switching at the base and propulsive base-to-tip propagation.

## 1. Introduction

Flagella and cilia undergo bending deformations under the action of dynein, a motor protein powered by ATP hydrolysis. To produce bending, dynein molecules form an array of cross-bridges between pairs of microtubule doublets that comprise the flagellar cytoskeleton (the axoneme), and exert forces that cause sliding of one doublet relative to the other. These active shear forces interact with passive structural elements (doublets, nexin links and radial spokes) to produce bending [1]. Dynein activity must be coordinated in order to produce oscillatory, propulsive waveforms.

The mechanism of dynein regulation has been an active field of investigation for many years. In a remarkable series of experimental and theoretical studies [2–16], Brokaw explored a number of potential feedback mechanisms. Hines & Blum [17] contributed a seminal paper in which a detailed continuum model of the flagellum was derived, including delayed curvature feedback. Later, Murase and co-workers [18–21] proposed the ‘excitable’ dynein concept, in which sliding beyond a specific threshold stimulates dynein activity. Julicher and co-workers [22–24] further developed the concept of sliding-controlled, collective dynein behaviour to explain flagellar oscillation. These authors postulated a positive feedback mechanism in which the force per dynein head decreases as sliding velocity increases, which allows more dynein to be recruited, thus increasing net shear force [22]. Lindemann [25–27] proposed the ‘geometric clutch’ (GC) model of dynein regulation, in which the spacing between doublets controls the level of dynein cross-linking. In the GC model, interdoublet spacing is affected by cumulative shear force and curvature, providing a plausible explanation for mechanical feedback. The details of flagellar synchronization and mechanics remain topics of active research.

The predictions of theoretical models of flagella mechanics and dynein regulation may be obtained by computer simulation or mathematical analysis. Typically, these models are expressed as partial differential equations (PDEs) for the flagellar angle as a function of time and axial position [17,22–24]. In computer simulation, time-marching algorithms are used to iteratively solve the discretized equations of motion [17,26]. Alternatively, closed-form solutions may be found for simplified versions of the equations. Solutions to linearized models, valid for small-amplitude oscillations, are composed of linear combinations of characteristic modes of oscillation (eigenfunctions). For example, in a mathematical model based on sliding-controlled dynein regulation [22,23], oscillatory (periodic) modes were found that closely matched experimental measurements of the flagellar waveform. However, in this prior study, only periodic (neutrally stable) modes were sought, leaving open the possibility of coexisting unstable (growing) or decaying modes.

This study explores the range of behaviour of recent mathematical models of flagellar motion. The solution methods from references [22–24] are generalized to find and characterize unstable modes. The paper briefly reviews flagellar mechanics, dynein regulation models (sliding-controlled, curvature-controlled and separation-controlled), and the associated eigenvalue problems. Solution methods for these eigenvalue problems are described briefly, and unstable and neutrally stable modes are found for each model. In particular, in recent sliding-controlled models, exponentially growing oscillatory solutions with retrograde (tip-to-base) propagation exist at the same physical parameter values shown to produce flagella-like, anterograde, periodic modes. In contrast, the least stable modes of the GC model exhibit switching at the base and propulsive base-to-tip propagation.

## 2. General equations of flagellar mechanics

### 2.1. Nonlinear equations of motion

Hines & Blum [17] derived the equations of motion for a flagellum in two-dimensions, modelling it as a slender elastic beam in viscous fluid, subjected to both active and passive internal shear forces (figure 1). These equations, summarized below, form the basis for the flagella models considered here [17,22–24].

The equations of force equilibrium, neglecting inertia, at any point along the flagellum are written in terms of the flagellar tangent angle, *ψ*, the net internal tangential and normal force components (*T*, *N*) and the external viscous force components per unit length (*q _{T}* and

*q*) [17,28]. 2.1and 2.2Similarly, the moment balance for any element along the flagellum may be written as 2.3where

_{N}*M*

_{B}is the moment owing to elastic bending,

*a*is the effective diameter and

*f*is the net interdoublet shear force (from distributed active dynein arms and passive elements such as radial spokes or nexin links). The velocity of any point along the flagellum can be written as

**=**

*v**v*+

_{N}**e**_{N}*v*. The spatial derivatives of the normal and tangential components of velocity are 2.4and 2.5Finally, constitutive properties are postulated for the beam and surrounding fluid. The flagella is modelled as a slender elastic beam with flexural rigidity

_{T}**e**_{T}*EI*2.6and the fluid at low Reynolds number is assumed to provide resistive force proportional to velocity [17] 2.7and 2.8The equilibrium, kinematic and constitutive equations can be combined to form two equations describing the motion of a slender elastic beam with internal shear moving in viscous fluid [17,23]. 2.9and 2.10where (·)

_{,a}=

*∂*(·)/

*∂a*.

### 2.2. Linearized equation and boundary conditions

To describe small-amplitude motion about a straight, equilibrium configuration, equation (2.10) may be linearized to obtain a much simpler equation [22,23,29] which can be written 2.11

Solutions to equation (2.11) must also satisfy appropriate boundary conditions. For example, if the flagellum is fixed at its proximal end (*s* = 0) and free at its distal end (s = *L*), solutions must satisfy

— (2.11) (i) Zero angle at base:

*ψ*(0,*t*) = 0— (2.11) (ii) Zero normal velocity at base:

*EIψ*_{,sss}(0,*t*) −*af*_{,s}(0,*t*) = 0— (2.11) (iii) Zero moment at distal end:

*EIψ*_{,s}(*L*,*t*) = 0— (2.11) (iv) Zero transverse force at distal end:

*EIψ*_{,ss}(*L*,*t*) −*af*(*L*,*t*) = 0

To solve the equation or perform stability analysis, specific models of dynein regulation may be used to express *f* in terms of *ψ* and system parameters.

## 3. Models of dynein regulation

Models of dynein regulation are equations that relate the interdoublet shear force *f*(*s*, *t*) to mechanical variables such as curvature or sliding velocity. Examples explored in this paper are the model of sliding-controlled regulation described in reference [22], the implementation of curvature-controlled feedback described by Hines & Blum [17], and the model of dynein control by interdoublet spacing: the ‘GC’ hypothesis, proposed originally by Lindemann [25–27].

### 3.1. Sliding-controlled dynein regulation

#### 3.1.1. Basic equations of sliding-controlled dynein regulation

Several recent studies [22–24] have suggested that interdoublet sliding provides the feedback necessary to produce sustained, propulsive flagellar oscillations. Interdoublet sliding displacement is related to bending by the kinematic relationship [17,22]
3.1where Δ_{0} = Δ(0, *t*) represents sliding permitted at the base of the flagellum.

In the hypothesized feedback mechanism, local interdoublet sliding reduces the load per dynein motor, leading to recruitment of more dyneins, and greater net force [22–24]. The detailed explanation of this hypothesis is given in reference [22] and is summarized briefly here. The active shear force is related to the incremental change in probability of dynein cross-linking, *δp*, and to the local sliding velocity, ∂Δ/∂*t*.
3.2

The parameters *ρ*, , and *f*′, which are all positive and real, have the following meaning: *ρ* is the linear density of dynein arms; is the maximum force per dynein arm, is the mean baseline attachment probability and *f*′ is the magnitude of the slope of the dynein force–velocity curve. Because the force per dynein head decreases with velocity, dynein cross-linking probability changes in response to sliding speed, as described by the equation
3.3

where *τ* is the characteristic time for this effect and *f*_{c} is a characteristic force for cross-link detachment. This relationship can lead to an increase in net shear force in response to increasing sliding rate, which is mathematically equivalent to negative friction or stiffness.

#### 3.1.2. Eigenvalue problem for sliding-controlled dynein regulation

To describe the response of the flagellum, characteristic modes of oscillation are sought. Generalizing the approach of references [22–24], separable solutions of the form
3.4are sought with *σ* = *α* + *iω* (*α* and *ω* are real). Each such solution that satisfies the equation of motion, and all boundary conditions, is a valid solution mode. If *α* > 0, the mode grows exponentially. If *M* such modes are found with exponents *σ*_{m} and shape, , then a solution can also be formed from any linear combination of these modes: . In general, for arbitrary initial conditions, the least stable mode will dominate the response.

Consistent with the assumed form for *ψ*, analogous expressions for force, sliding displacement and cross-linking probability, for example, are [22]
3.5

The sliding-controlled model of reference [22] can be expressed in terms of a complex mechanical impedance *χ*(*σ*), using equation (3.5) in equations (3.1)–(3.3):
3.6

The complex impedance, *χ*(*σ*), is defined in terms of dynein kinetics (following appendix B of [22]), by substituting the assumed solution form into the equation governing cross-linking probability (equation (3.3))
3.7

Again, following appendix B of reference [22], the expression for may then substituted into the equation for shear force (equation (3.2)) 3.8

The complex impedance, *χ*(*σ*), is now written compactly in terms of the characteristic exponent, *σ*
3.9
3.10
3.11where the derived parameters *A* and *B* are positive and real. Equation (3.9) corresponds to the ‘active’ part of equation B9 in [22]. For the special case of the neutrally stable, periodic response [22,23] the characteristic exponent *σ* = *iω*, and *χ*(*iω*) = *k* + *iω**λ*. The effects of the feedback-controlled dynein motors can be combined into dynamic stiffness (*k*) and friction (*λ*) coefficients
3.12Passive stiffness and friction may also be included, but these contributions are considered negligible by the authors of reference [22]; the values of *k* and *λ* are expected to both be negative.

Substitution of equations (3.5)–(3.6) into the linearized equation of motion equation (2.11), as in [22], leads to 3.13

Following [22], equation (3.13) is written in non-dimensional form as 3.14with non-dimensional parameters defined as in [22] 3.15The boundary conditions for the fixed-free case are also written in non-dimensional form [22]

— (3.16) (i) Zero angle at base:

— (3.16) (ii) Zero normal velocity at base:

— (3.16) (iii) Zero bending moment at distal end:

— (3.16) (iv) Zero transverse force at distal end:

where describes the interdoublet sliding at the base [22].

Two cases of the sliding-controlled model, described originally in references [22,23], are considered here.

Case 1: sliding at the proximal end is resisted by base stiffness, *k _{s}*, and friction,

*γ*

_{s}, leading to the expression [22] 3.17where 3.18Case 2: sliding at the proximal end of the flagellum is prohibited: .

### 3.2. Curvature-controlled dynein regulation

#### 3.2.1. Basic equations of curvature-controlled dynein regulation

Hines & Blum [17] implemented curvature control in a continuum model of the swimming flagellum of the form of equations (2.9)–(2.10). The shear force *f*(*s*, *t*) is defined by the expression
3.19where *S _{r}* is the shear force owing to deformation of passive components and

*S*

_{d}is the active dynein force. Local dynein shear force is dynamically regulated by curvature, according to the equation 3.20where

*m*

_{0}(N m) determines the effect of curvature feedback on dynein activity. Hines & Blum [17] suggest that the passive component of shear is related to interdoublet sliding displacement Δ(

*s*,

*t*) by the nonlinear relationship 3.21

#### 3.2.2. Eigenvalue problem for curvature-controlled dynein regulation

By direct analogy to the sliding-controlled model in §3.1.2, separable solution forms (equations (3.4)–(3.5)) are substituted into the equations for the curvature-controlled model of Hines & Blum [17] (equations (3.19)–(3.21) and (2.11)). The corresponding non-dimensional differential equation and boundary conditions for the fixed-free case are [17] 3.22

— (3.23) (i) Zero angle at base: .

— (3.23) (ii) Zero normal velocity at base: .

— (3.23) (iii) Zero bending moment at distal end: .

— (3.23) (iv) Zero transverse force at distal end:

As before, . The new non-dimensional variables are (the ratio of the dynein time constant to the viscoelastic time constant of the flagellum) and (the ratio of the characteristic dynein moment *m*_{0} to the moment required to bend the flagellum to curvature 1/*L*). Note that passive shear forces may be added to account for elastic or dissipative shear elements, but in the linearized version of the original model, the passive shear-restoring force (equation (3.21), to linear order) is zero.

### 3.3. Dynein regulation by interdoublet spacing: the geometric clutch hypothesis

#### 3.3.1. Basic equations of the geometric clutch model

The GC model was proposed by Lindemann [25–27] and implemented as a discretized computer simulation. A version of the model, expressed as a set of PDEs, was derived recently [28] and is briefly summarized here. In the GC model, the net shear force is *f* = *f*_{P} + *f*_{R}, where *f*_{P} and *f*_{R} are the net shear forces in the principal and reverse directions, exerted between pairs of doublets on opposite sides of the axoneme. Including shear forces owing to passive stiffness and friction, net shear force is thus expressed as [28]
3.24

The probabilities *p*_{0} and *p*_{1} are system parameters that represent the baseline and maximum probabilities of dynein cross-linking, and is the maximum dynein force per unit length. The variable *A* = *A*_{P} − *A*_{R} represents net dynein activity in the principal bend direction. The parameters *k _{T}* and

*b*represent stiffness and damping in shear. (Note the sign convention for shear force here is opposite that of reference [28], and consistent with references [22,23].)

_{T}In the GC model, dynein activity is coupled to global flagellar motion by tension and curvature [28]. The distributed shear force causes a difference in tension between doublets in each active pair: and . If the doublets are curved, these tension differences lead to a component of force ( or ) that separates the doublets or draws them together. After adding the two sides together and linearizing, only the baseline difference in tension, *S*_{0}, appears in the equation governing net dynein activity, *A*. The following expressions are obtained [28]
3.25and
3.26The time constant *τ*_{N} describes the local dynein kinetics. The resting ‘isometric’ difference in tension, *S*_{0}, provides a baseline level of coupling between curvature and dynein activity even when the flagellum is almost straight. The parameter *C _{S}* controls the magnitude of the coupling [28]. The stability of a straight flagellum is significantly affected when the baseline tension difference

*S*

_{0}> 0.

#### 3.3.2. Eigenvalue problem for the geometric clutch model

Solutions to the flagella equations with GC dynein regulation are sought in the separable form above (equations (3.5)), and substituted into the expressions for dynein activity (equations (3.24)–(3.26)) and flagellar motion (equation (2.11)). The resulting expressions are combined and simplified to obtain the following equation for 3.27with coefficients 3.28and 3.29In non-dimensional form, the equation becomes 3.30or its equivalent, 3.31where the new non-dimensional parameters are and .

The corresponding boundary conditions are

— (3.32) (i) Zero angle at base:

— (3.32) (ii) Zero normal velocity at base:

— (3.32) (iii) Zero bending moment at distal end:

— (3.32) (iv) Zero transverse force at distal end:

By comparing equation (3.31) with the analogous equations describing sliding-controlled models (equation (3.14)) [22–24] and curvature-controlled models (equation (3.22)) [4,10,17,30], it is apparent that the GC model includes feedback from both curvature (the term containing ) and shear deformation (the term containing ). Notably, both terms become more destabilizing when *c*_{1}(*σ*) increases [28].

An important feature of the curvature-feedback term in the GC model is that it is proportional to . The feedback is thus strongest at the proximal end, which encourages switching at the base and thus proximal-to-distal propagation. The factor of also complicates the solution of the eigenvalue problem, so that numerical methods (weighted residual or finite-element calculations, e.g.) are required to find the natural modes and frequencies of oscillation.

### 3.4. Nonlinear-restoring force

The equations of sections 3.1–3.3 describe linearized models of the forces between doublets. These models either omit elastic shear forces or permit elastic force to remain proportional to displacement. In fact, it is likely that forces from passive components of the axoneme will increase nonlinearly with displacement [23]. Such nonlinearity does not affect the linear models above, but prevents unstable modes from growing without limit. Accordingly, in simulations, an additional nonlinear-restoring shear force proportional to the cube of the shear displacement was added. In each case, the shear force is 3.33

where Δ(*s*, *t*) is defined by equation (3.1), and *f*(*s*, *t*) by equations (3.2), (3.19) or (3.24). The nonlinear-restoring force of equation (3.33) approximates the assumed stiffening behaviour of elastic elements [4,23]. When the linearized model is unstable, the value of *k*_{3} may determine the amplitude of simulated oscillations.

## 4. Solution of the eigenvalue problem

The spatial differential equations derived in sections 3.1–3.3 are summarized here.

Sliding-controlled: , (3.14).

Curvature-controlled: , (3.22).

GC: . (3.31).

The behaviour of each model, for a given set of boundary conditions, can be described in terms of the characteristic modes, or eigenfunctions. Exact expressions for the eigenfunctions in terms of exponential functions can be found for the sliding-controlled and curvature-controlled models, as described below (§4.1); the eigenvalues themselves must be found numerically. The GC model requires a numerical approach to find eigenvalues and eigenfunctions. Accordingly, the method of weighted residuals and the finite-element method, both of which may be applied to all models, are invoked (§4.2).

### 4.1. Eigenvalue analysis by substitution of an assumed exponential solution

#### 4.1.1. Solution of the eigenvalue problem for the sliding-controlled model

The eigenvalue problem for the sliding-controlled model (equation (3.14) together with the boundary conditions equation (3.16)) may be solved as follows [22]. The differential equation for the mode shape, equation (3.14) is satisfied by exponential solutions of the form 4.1

Substitution of the assumed form into the equation of motion leads to the characteristic polynomial 4.2

For a given value of , the general solution to equation (3.14) can be constructed using the four roots of equation (4.2) 4.3

Equations (3.17) and (3.18) can be used to eliminate from the boundary condition equations, leading to a matrix equation of the form below (a representative column of the 4 × 4 matrix is shown): 4.4

If sliding at the base is not permitted, the base compliance parameter *Γ* = 0.

Note that the values of *β*_{n} in the matrix above depend not only on the physical parameters of the model, but also on the characteristic exponent, or eigenvalue, . For a specific parameter set, the eigenvalue problem can thus be written compactly as
4.5

Non-trivial solutions are found by seeking values of that lead to 4.6

Eigenvalues are found at the zeros of (or minima of ). Mode shapes are found by seeking vectors *a*^{(m)} that span the corresponding null-space of the matrix, ; this is accomplished by singular value decomposition. The mode shape is then reconstructed by substituting values for and (coefficients and roots corresponding to the eigenvalue ) into equation (4.3). Results may be expressed in dimensional form using .

#### 4.1.2. Solution of the eigenvalue problem for the curvature-controlled model

Substituting the assumed exponential form (equation (4.1)) into the equation of the curvature-controlled model, the characteristic polynomial is found to be 4.7

By imposing boundary conditions, the eigenvalue problem is expressed as the matrix equation (4.8) (again a representative column of the 4 × 4 matrix is shown) 4.8As in the sliding-controlled model, characteristic exponents (eigenvalues) and mode shapes (eigenfunctions) are obtained by finding the zeros of the determinant of the matrix on the left side of equation (4.8) and reconstructing the mode shapes using equation (4.3).

#### 4.1.3. The eigenvalue problem for the geometric clutch model

It is not possible to use the approach above to find closed-form expressions for the eigenfunctions for the GC model, owing to the factor in the equation for . Instead, approximate solutions to the eigenvalue problem were sought using numerical methods. Such methods include finite-element analysis, or the method of weighted residuals [31]. The method of weighted residuals is well-suited to this problem, as the flagellar eigenfunctions can be constructed from the vibration modes of an Euler–Bernouilli beam with corresponding boundary conditions. The application of the method of weighted residuals to this problem is described in the electronic supplementary material, part A.

### 4.2. Finite-element eigenanalysis and time-marching simulation

To check the stability analysis and characterize the corresponding behaviour, solutions to the equations of motion were also found using the PDE modelling capability of a commercial finite-element (FE) software package (COMSOL v. 4.3a, COMSOL, Inc., Burlington, MA). The one-dimensional domain was discretized into 50 elements with quartic interpolation. Eigenvalue/eigenfunction calculations (300 maximum iterations, relative tolerance 1 × 10^{−6}) and time-marching simulations (backward differentiation formula, variable time step, relative tolerance 1 × 10^{−4}) were performed. Representative results were confirmed at finer spatial resolution and smaller tolerance values.

## 5. Results: unstable and neutrally stable modes of flagellar models

For each model, we seek unstable modes that lead to oscillation. Each such mode is characterized by frequency, shape and propagation direction, which determines propulsive effectiveness. The basic physical parameters (length, flexural rigidity, diameter, viscosity) and fixed-free boundary conditions are consistent between the models. In the first sliding-controlled case, interdoublet sliding is permitted at the base, for comparison with the results of reference [22]. The role of unstable modes in initiation and maintenance of flagellar oscillations is investigated by numerical simulations of fully nonlinear models.

### 5.1. Unstable and neutrally stable modes of sliding-controlled models

#### 5.1.1. Case 1: modes of the sliding-controlled model with finite sliding compliance at the base

*Modes from assumed exponential solution:* a sliding-controlled model with fixed-free boundary conditions and sliding permitted at the base [22,23] was analysed by the assumed solution approach described above. This case, and corresponding parameters (table 1), are chosen for comparison to the analogous example in reference [22]. Sliding at the base is permitted, opposed by finite positive stiffness (*k _{s}* = 94.8 × 10

^{−}^{3}N m

^{−1}) and friction (

*γ*

_{s}= 0.273 × 10

^{−3}N-s m

^{−1}) as in [22]. Other parameters were also chosen to match the mechanical impedance used in reference [22]. For example, these parameters lead to negative effective shear stiffness (

*k*= −1620 N m

^{−2}) and friction (

*λ*= −7.60 N-s m

^{−2}) in the flagellum for the 20.6 Hz mode discussed in [22].

Characteristic exponents in the complex plane are found at local minima of the determinant magnitude, |*D*(*σ*)|, shown in figure 2*a*,*b* (results are shown in terms of physical, rather than dimensionless, variables). Two notable observations can be made. (i) A periodic solution (complex conjugate eigenfunctions with purely imaginary characteristic exponents: *σ* = *iω*) exists. The frequency of this periodic mode corresponds precisely to the frequency (*ω*/2*π* = 20.6 Hz) of the periodic mode reported in reference [22] for these parameters. (ii) Multiple unstable modes coexist with this periodic mode. Three of these unstable modes have characteristic exponents with positive real parts (*Re*(*σ*) = *α* > 0). Because of these coexisting unstable modes, the physical significance of the neutrally stable 20.6 Hz mode is minimal.

*Results from the method of weighted residuals*: complementary results obtained by the method of weighted residuals are shown in figure 2*c.* Figure 2*c* shows the paths of the eigenvalues in the complex plane as the baseline probability of dynein attachment, , is increased from 0 to 0.04. The final values of the eigenvalues (red ‘x’ symbols in figure 2*c*) agree closely with the values obtained from the local minima of |*D*| (figure 2*a,b*).

Mode shapes corresponding to these eigenvalues are shown in figure 3. The shape of the periodic mode at 20.6 Hz obtained here matches closely the shape of the periodic mode at the same frequency in reference [22]. The unstable modes of this model under the same conditions have not been described in prior studies; they are expected to influence physical behaviour more than the coexisting periodic mode.

Figure 4*a* shows the frequency of the least stable, unstable mode, as function of the mean cross-linking probability, and flagellum length, *L*. The white region of the plot indicates parameter combinations for which unstable modes are absent, so the edge of the coloured region indicates the stability boundary. Solutions with zero imaginary part (purely real eigenvalues *α* > 0) correspond to non-oscillatory unstable modes.

Figure 4*b* displays the gradient of phase for the least stable mode, , at each parameter combination. If the mode exhibits proximal-to-distal (anterograde) propagation, the phase gradient will be negative. For example, the phase gradient of the 20.6 Hz mode is negative. However, the phase gradient of the least stable mode is positive for all parameter combinations in this model, indicating that retrograde (distal-to-proximal) propagation is present in this model with these parameters.

#### 5.1.2. Case 2: stability of the sliding-controlled model with no sliding at the base

This case provides an opportunity for more extensive comparison to results from prior work [23] (figure 5). Periodic modes were found by the assumed solution method (§4.1) at frequencies from 1 to 100 Hz, using the same parameter values as in reference [23]. The values of the complex impedance *χ* corresponding to these periodic solutions are shown in figure 5 (red markers); the frequency of each mode is plotted on the vertical axis of figure 5*b*. Figure 5*a* is directly comparable to fig. 3 in [23].

Also shown in figure 5 are the frequencies of all unstable modes found by the method of weighted residuals, at all values of the complex impedance, *χ*. The frequency predictions of the weighted residual method at the boundaries of different solution regimes agree closely with the frequencies of periodic modes obtained by the assumed solution approach. These predictions also match corresponding results in reference [23]. However, for many of the parameter combinations at which periodic solutions exist, coexisting unstable modes are also found.

Figure 6*a* shows the locations of eigenvalues of this model as the mean probability of cross-linking, , is varied (other parameters are as in table 2). The frequency and direction of the least stable mode are shown in figure 6*b*,*c* as functions of the length of the flagellum (*L*) and . The neutrally stable mode at 28 Hz, identified by eigenvalues on the imaginary axis in figure 6 and by the mode shape in figure 7, matches the corresponding example result from [23].

### 5.2. Unstable modes of the curvature-controlled model

Eigenvalues, *σ*, of the classic curvature-controlled, fixed-free flagellum model of Hines & Blum [17], using parameter values in table 3, are shown in figure 8. The locations of these eigenvalues in the complex plane are found at local minima of |*D*(*σ*)|, which is shown in figure 8*a* and expanded in figure 8*c*. Notably, the least stable solutions are oscillatory with the imaginary parts of the characteristic exponents corresponding to frequencies of flagellar motion. The paths of eigenvalues in the complex plane, obtained by the method of weighted residuals, are shown in figure 8*b* (expanded in figure 8*d*) with the final eigenvalues closely matching the minima of |*D*(*σ*)| in figure 8*a*,*c*. The mode shape corresponding to the unstable positive eigenvalue is shown in figure 9.

Oscillation frequencies and phase gradients for the least stable modes of the curvature-controlled model are shown in figure 10. As the curvature feedback parameter *m*_{0} is decreased the modes become stable, as is characteristic of the underlying Hopf bifurcation. All the unstable modes in this parameter range exhibit anterograde (proximal-to-distal) propagation.

### 5.3. Unstable modes of the geometric clutch model

Figure 11 illustrates the increasing baseline probability of dynein attachment, *p*_{0} on the eigenvalues of the GC model of the flagellum. A dynamic instability occurs when complex eigenvalues cross into the right half plane with non-zero imaginary part: *Re*(*σ*) > 0; *Im*(*σ*) ≠ 0. Increasing the baseline probability of cross-bridge attachment, *p*_{0}, encourages instability. The modes corresponding to the least stable eigenvalues at *p*_{0} = 0.05 (with other parameters as in table 4) are shown in figure 12.

### 5.4. Results of finite-element analysis and numerical simulation

Eigensolutions were also obtained by FE analysis of the linearized models using COMSOL PDE modelling software, as described above. Characteristic exponents and mode shapes found by FE eigenanalysis (not shown) correspond closely to those found by the method of weighted residuals, shown in figure 3, figure 7, figure 9 and figure 12.

Numerical simulations (i.e. time-marching) of the full nonlinear flagellum equations in the time domain with fixed-free boundary conditions were performed in COMSOL to investigate the behaviour of the dynein regulation models at finite amplitudes. Time-marching simulation allows exploration of nonlinear, transient and non-periodic behaviour. Simulations were performed for the sliding-controlled model (equations (2.9) and (2.10) with equations (3.1)–(3.3)), for the curvature-controlled model (equations (2.9) and (2.10) and equations (3.19)–(3.21)), and for the GC model (equations (2.9) and (2.10) and equations (3.24)–(3.26)). A nonlinear elastic shear force (equation (3.33)) was included in each model for simulation; the value of the nonlinear coefficient *k*_{3} was adjusted to control amplitude (i.e. to produce oscillations of comparable amplitude in each model).

Simulation of each system led to oscillatory solutions with moderate amplitude (figure 13). Each row of figure 13 shows (i) snapshots of the waveform colour-coded by time; (ii) time series of angle at the tip and tension at the base; (iii) ‘phase plots’ of flagellar angle at *s* = 3*L*/4 and *s* = *L* and (iv) the fundamental mode obtained by Fourier analysis. Simulations of the sliding-controlled model with sliding at the base (case 1) lead to oscillatory waves with retrograde propagation. The positive mean value of tension at the base (*T*_{0}) signifies backward propulsive force. The sliding-controlled model without sliding at the base (case 2) exhibits non-propulsive, retrograde waves similar to the least stable mode of the linearized version. The curvature-controlled model exhibits clear anterograde wave propagation that leads to a forward propulsive force (negative *T*_{0}). The GC model exhibits propulsive, anterograde propagating waves that resemble the single unstable mode.

Fourier analysis of the steady-state spatio-temporal pattern of *ψ*(*s*, *t*) from simulation was performed along the time dimension (*fft*; Matlab^{™}, The MathWorks, Natick, MA) to obtain the frequency and shape of the fundamental mode of oscillation (figure 13 and table 5). The relative contribution of the fundamental mode to the simulated response was measured by the ratio of its magnitude to the summed magnitudes of all the Fourier coefficients (table 5). The similarity of the fundamental mode from simulation to the unstable modes of the linearized model (obtained by weighted residuals eigenanalysis) was measured by the magnitude of the correlation coefficient between the two shapes (*corrcoef*; Matlab^{™}, The MathWorks; table 5).

Finally, the sensitivity of flagellar behaviour to model parameters can be determined from the changes in eigenvalues as the parameter is varied. Electronic supplementary material, figure S5 shows the effects of flagellar length, *L*, flexural rigidity, *EI* and resistive force coefficient, *c _{N}* on eigenvalues of each model.

## 6. Discussion and conclusion

Unstable and neutrally stable periodic modes were identified for the sliding-controlled [22–24], curvature-controlled [17], and doublet separation-controlled (GC) [25–27] models of flagellar motion. Sliding-controlled models are characterized by sliding impedances with negative effective stiffness and friction coefficients. Previous studies of sliding-controlled models [22–24] have focused on the existence of and characterization of the neutrally stable periodic solutions. The current analysis confirms the existence of periodic modes that closely resemble observed flagellar behaviour. However, the uniqueness of these solutions was not demonstrated in prior studies, leaving open the possibility of coexisting unstable modes. In this study, we find that at the same parameter values that give rise to neutrally stable periodic modes in the sliding-controlled model, multiple unstable modes exist which exhibit retrograde (distal-to-proximal) wave propagation. Such unstable modes would dominate the response of a physical system.

The fixed-free, curvature-controlled model, derived originally by Hines & Blum [17] also exhibits an unstable mode at the parameter values in table 3. This unstable mode was characterized by anterograde, propulsive bend propagation. The dynamically unstable mode found in the present study parallels the propulsive modes of oscillation found by time-marching simulation in reference [17]. The existence of propulsive waveforms does not imply that this particular curvature-controlled model of flagellar motion is completely satisfactory. Others, including the authors of the original study [17], have noted that this model of curvature control is not based on a specific biophysical mechanism, and the original model relies on parameter values (particularly flexural rigidity) that may not be accurate. However, the ability to model stable, propulsive, oscillatory solutions as well as transient behaviour confirms the value of this seminal model.

The GC model exhibits a dominant unstable mode at the parameter values used here. The mode was characterized by clear anterograde waveform propagation. Physically, oscillations of the GC model are the result of a switching mechanism that is strongest at the base of the flagellum [26,32]. Active shear in one doublet pair (for example the P side) induces a tension difference in the doublets on that side, which combined with curvature, produces a transverse force which pushes the doublets apart and eventually terminates the active shear. This transverse force corresponds to the ‘global transverse force’ described by Lindemann [25–27]. At the same time, the corresponding passive shear force on the opposite doublet pair produces a transverse force which pulls the doublets on the passive side together and initiates active shear on that side. The GC model combines a physically intuitive mechanism with predicted behaviour that resembles observed waveforms.

Why do waves propagate from base to tip in some models but not others? The effect of dynein regulation mechanism on propagation direction is discussed appendix D of reference [22]. To summarize the result: in models with curvature feedback, propagation direction is determined by the phase of the (complex) curvature coefficient (equation D7 of reference [22]). In models with sliding control, propagation direction is not determined by the equation of motion. Thus, boundary conditions strongly affect propagation direction in sliding-controlled models. We hypothesize that in fixed-free, sliding-controlled models, larger sliding amplitudes at the free end may encourage early switching at the tip and retrograde propagation.

Nonlinear versions of these models, which more closely approximate the physical situation, were explored by time-domain simulation. Even in the nonlinear regime, unstable modes of the linearized models still appear to have pronounced effects on observed behaviour (figure 13 and table 5). When unstable modes exist, damped (decaying) or neutrally stable modes have relatively little influence on the large-amplitude behaviour. In addition, when multiple, strongly unstable modes exist (as in the sliding-controlled model, case 1), the influence of each individual mode is less clear. The results of this study complement the recent observation that distinct nonlinear modes of deformation in flagellar models may arise at large amplitudes [33] leading to asymmetric waveforms.

The linearized models in this paper are expected to be quantitatively accurate (say within less than 10%) only for oscillations with angular amplitudes of approximately 0.3–0.4 radians, or less. Oscillation amplitudes are reported to be approximately 1 radian [22] in bull sperm and even larger, asymmetric oscillations are seen in flagella of *Chlamydomonas reinhardtii* [34]. While eigenanalysis of linearized models illuminates the initiation and maintenance of propagating waves, and provides reasonable estimates of frequency, it does not provide any information about oscillation amplitudes. Nonlinear modelling of flagellar motion is a compelling target of future work. The geometric nonlinearity associated with large-amplitude motion of the slender Euler–Bernouilli beam is well understood, but much less is known about the possible structural and material nonlinearities of the axoneme.

In conclusion, the unstable modes of linear models of flagellar motion illuminate the hypotheses that define each model of dynein regulation, and clarify the consequences of underlying assumptions and parameter choices. With reasonable parameter values, the fixed-free GC model exhibits unstable modes that are propulsive and propagate base-to-tip. In contrast, the existence of unstable retrograde (tip-to-base) modes in current sliding-controlled models, with both fixed and sliding boundary conditions at the base, appears to weaken the evidence for this hypothesis of dynein coordination. The development of models of dynein regulation and flagellar motion remains an active topic of research. New details of the mechanics of flagella are still being uncovered [35] as are the mechanisms of synchronization between flagella [36,37]. The stability properties and propagation directions of all modes should be considered in evaluating proposed mathematical models and their underlying hypotheses.

## Acknowledgements

Support from NSF grant no. CMMI-1265447 and the Children's Discovery Institute is gratefully acknowledged. The authors thank Anders Carlsson, Susan Dutcher and Ruth Okamoto for helpful discussions.

- Received February 11, 2015.
- Accepted March 12, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.