Journal of The Royal Society Interface

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 [216], 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 [1821] proposed the ‘excitable’ dynein concept, in which sliding beyond a specific threshold stimulates dynein activity. Julicher and co-workers [2224] 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 [2527] 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,2224]. 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 [2224] 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,2224].

Figure 1.

Figure 1. Schematic of a prototypical flagellum showing tangent angle, ψ(s, t). (b) Magnified view of a section of flagellum showing the sliding displacement, Δ(s, t), internal shear force f(s, t) and effective diameter a. In some sliding-controlled models, sliding at the base Δ(0, t) is permitted even if transverse motion of the base is constrained.

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 (qT and qN) [17,28].

Display Formula
2.1
and
Display Formula
2.2
Similarly, the moment balance for any element along the flagellum may be written as
Display Formula
2.3
where MB 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 = vNeN + vTeT. The spatial derivatives of the normal and tangential components of velocity are
Display Formula
2.4
and
Display Formula
2.5
Finally, constitutive properties are postulated for the beam and surrounding fluid. The flagella is modelled as a slender elastic beam with flexural rigidity EI
Display Formula
2.6
and the fluid at low Reynolds number is assumed to provide resistive force proportional to velocity [17]
Display Formula
2.7
and
Display Formula
2.8
The 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].
Display Formula
2.9
and
Display Formula
2.10
where (·),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

Display Formula
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 [2527].

3.1. Sliding-controlled dynein regulation

3.1.1. Basic equations of sliding-controlled dynein regulation

Several recent studies [2224] 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]

Display Formula
3.1
where Δ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 [2224]. 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.

Display Formula
3.2

The parameters ρ, Inline Formula, Inline Formula and f′, which are all positive and real, have the following meaning: ρ is the linear density of dynein arms; Inline Formula is the maximum force per dynein arm, Inline Formula 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

Display Formula
3.3

where τ is the characteristic time for this effect and fc 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 [2224], separable solutions of the form

Display Formula
3.4
are sought with σ = α + (α 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, Inline Formula, then a solution can also be formed from any linear combination of these modes: Inline Formula. 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]

Display Formula
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):

Display Formula
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))

Display Formula
3.7

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

Display Formula
3.8

The complex impedance, χ(σ), is now written compactly in terms of the characteristic exponent, σ

Display Formula
3.9
Display Formula
3.10
Display Formula
3.11
where 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 σ = , and χ() = k + λ. The effects of the feedback-controlled dynein motors can be combined into dynamic stiffness (k) and friction (λ) coefficients
Display Formula
3.12
Passive 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

Display Formula
3.13

Following [22], equation (3.13) is written in non-dimensional form as

Display Formula
3.14
with non-dimensional parameters defined as in [22]
Display Formula
3.15
The boundary conditions for the fixed-free case are also written in non-dimensional form [22]

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

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

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

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

where Inline Formula 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, ks, and friction, γs, leading to the expression [22]

Display Formula
3.17
where
Display Formula
3.18
Case 2: sliding at the proximal end of the flagellum is prohibited: Inline Formula.

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

Display Formula
3.19
where Sr is the shear force owing to deformation of passive components and Sd is the active dynein force. Local dynein shear force is dynamically regulated by curvature, according to the equation
Display Formula
3.20
where m0 (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
Display Formula
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]

Display Formula
3.22

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

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

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

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

As before, Inline Formula. The new non-dimensional variables are Inline Formula (the ratio of the dynein time constant to the viscoelastic time constant of the flagellum) and Inline Formula (the ratio of the characteristic dynein moment m0 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 [2527] 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 = fP + fR, where fP and fR 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]

Display Formula
3.24

The probabilities p0 and p1 are system parameters that represent the baseline and maximum probabilities of dynein cross-linking, and Inline Formula is the maximum dynein force per unit length. The variable A = APAR represents net dynein activity in the principal bend direction. The parameters kT and bT 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].)

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: Inline Formula and Inline Formula. If the doublets are curved, these tension differences lead to a component of force (Inline Formula or Inline Formula) that separates the doublets or draws them together. After adding the two sides together and linearizing, only the baseline difference in tension, S0, appears in the equation governing net dynein activity, A. The following expressions are obtained [28]

Display Formula
3.25
and
Display Formula
3.26
The time constant τN describes the local dynein kinetics. The resting ‘isometric’ difference in tension, S0, provides a baseline level of coupling between curvature and dynein activity even when the flagellum is almost straight. The parameter CS controls the magnitude of the coupling [28]. The stability of a straight flagellum is significantly affected when the baseline tension difference S0 > 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 Inline Formula

Display Formula
3.27
with coefficients
Display Formula
3.28
and
Display Formula
3.29
In non-dimensional form, the equation becomes
Display Formula
3.30
or its equivalent,
Display Formula
3.31
where the new non-dimensional parameters are Inline Formula and Inline Formula.

The corresponding boundary conditions are

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

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

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

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

By comparing equation (3.31) with the analogous equations describing sliding-controlled models (equation (3.14)) [2224] 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 Inline Formula) and shear deformation (the term containing Inline Formula). Notably, both terms become more destabilizing when c1(σ) increases [28].

An important feature of the curvature-feedback term in the GC model is that it is proportional to Inline Formula. The feedback is thus strongest at the proximal end, which encourages switching at the base and thus proximal-to-distal propagation. The factor of Inline Formula 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

Display Formula
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 k3 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: Inline Formula, (3.14).

Curvature-controlled: Inline Formula, (3.22).

GC:  Inline Formula. (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

Display Formula
4.1

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

Display Formula
4.2

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

Display Formula
4.3

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

Display Formula
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, Inline Formula. For a specific parameter set, the eigenvalue problem can thus be written compactly as

Display Formula
4.5

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

Display Formula
4.6

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

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

Display Formula
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)

Display Formula
4.8
As in the sliding-controlled model, characteristic exponents (eigenvalues) and mode shapes (eigenfunctions) are obtained by finding the zeros of the determinant Inline Formula 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 Inline Formula in the equation for Inline Formula. 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 (ks = 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].

Table 1.Parameter values for sliding-controlled model with base sliding—case 1 [22].

parameter value units description
cN 0.0034 pN-s μm−2 normal resistive force coefficient
cT 0.0017 pN-s μm−2 tangent resistive force coefficient (ctcN/2)
a 0.185 µm effective diameter of flagellum
L 58.3 µm length of flagellum
EI 1700 pN-μm2 flexural rigidity of flagellum
τ 0.004 s dynein time constant
Inline Formula 0.03 1 mean probability of cross-linking
Inline Formula 3.8 pN dynein stall force
fc 2.0 pN characteristic force
f 1.8 pN-s μm−1 slope of dynein force–velocity curve
ks 94.8 × 10−3 pN μm−1 sliding stiffness at base
γs 0.273 × 10−3 pN-s μm−1 sliding friction at base
ρ 150 1 μm−1 density of dynein motors
k3 800 pN μm−4 nonlinear shear stiffness

Characteristic exponents in the complex plane are found at local minima of the determinant magnitude, |D(σ)|, shown in figure 2a,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: σ = ) 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.

Figure 2.

Figure 2. (a) Eigenvalues from the assumed exponential solution of the sliding-controlled flagellar model (equation (3.14)) with sliding permitted at the base (case 1; parameters in table 1). Images of the log magnitude of the determinant, ln|D(σ)| (arb. units), are shown as a function of σ = α + . Eigenvalues (characteristic exponents) are found at local minima (blue) of |D(σ)|. Unstable modes have α = Re(σ) > 0. (b) Expanded view of panel (a) shows the eigenvalue at σ = i2π · 20.6 corresponding to a 20.6 Hz periodic mode. (c) Eigenvalues from the weighted-residuals method: paths of eigenvalues σ = α + are shown in the complex plane as Inline Formula is varied (Inline Formula). Other parameters are as in table 1. The red ‘x’ symbols denote the eigenvalues at the final value Inline Formula. The eigenvalues in panel (c) closely match the minima of |D(σ)| in panel (a). (Online version in colour.)

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

Figure 3. Unstable and neutrally stable modes of the sliding-controlled model with sliding at the base (case 1; equation (3.14)). (a,c,e,g) The mode shape expressed in terms of tangent angle, Inline Formula. (b,d,f,h) The mode shape expressed in terms of displacement Inline Formula The real (solid line) and imaginary (dashed line) part of each mode is shown. Each mode corresponds to eigenvalue σ = a + : (a,b) α = 51.3 s−1, ω/2π = 3.3 Hz (unstable); (c,d) α = 22.8 s−1, ω/2π = 22.5 Hz (unstable); (e,f) α = 4.0 s−1, ω/2π = 30.8 Hz (unstable); (g,h) α = 0.0 s−1, ω/2π = 20.6 Hz (periodic). See the electronic supplementary mateial, movies S1–S2 for corresponding animations.

Figure 4a shows the frequency of the least stable, unstable mode, as function of the mean cross-linking probability, Inline Formula 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.

Figure 4. (a) Frequency ω/2π (Hz) of the least stable mode of the sliding-controlled flagellar model (case 1; equation (3.14)) as a function of flagellar length, L, and mean probability of cross-linking, Inline Formula. Other parameters are as in table 1. At each parameter combination (Inline Formulafrequency is obtained from the imaginary part of the eigenvalue σ = α + with largest real part (α). (b) Median phase gradient, Inline Formula, of the least stable mode. Anterograde (proximal–distal) propagation corresponds to a phase gradient <0, for ω > 0. For all parameter combinations shown here, the least stable mode exhibits phase gradient >0 and thus retrograde propagation.

Figure 4b displays the gradient of phase for the least stable mode, Inline Formula, 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 5b. Figure 5a is directly comparable to fig. 3 in [23].

Figure 5.

Figure 5. Comparison of weighted residual stability predictions of frequency, ω/2π (Hz), to the frequency of periodic solutions obtained from the assumed exponential solution for the sliding-controlled model with no base sliding (case 2, described in [23]). Results are shown as a function of the mechanical impedance, χ. The colour scale in panels (a,b), and the vertical axis in panel (b) display the frequency of the least stable mode. Modes with zero frequency are non-oscillatory. The plots show good agreement as periodic (neutrally stable) solutions arise at transitions between solution regimes.

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 6a shows the locations of eigenvalues of this model as the mean probability of cross-linking, Inline Formula, is varied (other parameters are as in table 2). The frequency and direction of the least stable mode are shown in figure 6b,c as functions of the length of the flagellum (L) and Inline Formula. 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].

Table 2.Parameter values changed for sliding-controlled model with no base sliding—case 2 [22,23]. (L, EI, a, cN are as in table 1.)

parameter value units description
Inline Formula 0.01 1 mean probability of cross-linking
Inline Formula 5.2 pN dynein stall force
fc 1.0 pN characteristic force
f 0.5 pN-s μm−1 slope of dynein force–velocity curve
k3 40 pN μm−4 nonlinear shear stiffness
Figure 6.

Figure 6. Eigenvalues of the sliding-controlled model from the weighted residual method: (a) paths of eigenvalues σ = α + in the complex plane as baseline probability of dynein attachment, Inline Formula, is varied (Inline Formula) in the sliding-controlled flagellar model (equation (3.14)) with no sliding at the base (case 2). Other parameters are as in table 2. The red ‘x’ symbols denote the eigenvalues at the final value: Inline Formula. (b) Frequency ω/2π (Hz) of the least stable mode of this model as a function of flagellar length and Inline Formula. Other parameters are as in table 2. At each parameter combination (Inline Formulafrequency is obtained from the imaginary part () of the eigenvalue σ = α + with largest real part (α). (c) Median phase gradient, Inline Formula, of the least stable mode. Anterograde (proximal–distal) propagation corresponds to a phase gradient <0, for ω > 0. For all parameter combinations shown here, the least stable mode exhibits phase gradient > 0 and thus retrograde propagation.

Figure 7.

Figure 7. The 28 Hz, neutrally stable, retrograde propagating mode of the sliding-controlled model with no sliding at the base (case 2, described in [23]). (a) The mode shape expressed in terms of tangent angle, Inline Formula. (b) The mode shape expressed in terms of displacement Inline Formula The real (solid line) and imaginary (dashed line) part of each mode is shown. This periodic mode corresponds to eigenvalue σ = α + : α = 0.0 s−1, ω/2π = 28 Hz. See electronic supplementary material, movie S3 for corresponding animation.

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 8a and expanded in figure 8c. 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 8b (expanded in figure 8d) with the final eigenvalues closely matching the minima of |D(σ)| in figure 8a,c. The mode shape corresponding to the unstable positive eigenvalue is shown in figure 9.

Table 3.Parameter values for curvature-controlled model [17]. (L, EI, a, cN are as in table 1.)

parameter value units description
τ 0.015 s dynein time constant
m0 1600 pN-μm characteristic bending moment
k3 150 pN µm−4 nonlinear shear stiffness
Figure 8.

Figure 8. (a,c) Eigenvalues from the assumed exponential solution of the curvature-controlled, fixed-free flagellar model (equation (3.22)): images of the log magnitude of the determinant (ln|D(σ)|) (arb. units) are shown as a function of σ = α + . Values of key parameters are defined in table 3. Eigenvalues (characteristic exponents) are found at local minima (blue) of |D(σ)|. Unstable modes have α = Re(σ) > 0. (b) Eigenvalues from the weighted residuals method: paths of eigenvalues σ = α + are shown in the complex plane as m0 is varied (0 < m0 < 1600 pN-µm). Other parameters are as in table 3. The red ‘x’ symbols denote the eigenvalues at the final value m0 = 1600 pN-µm. (c,d) Expanded views of panels (a,b) showing the least stable eigenvalues. The eigenvalues in panels (b,d) closely match the minima of |D(σ)| in panels (a,c). (Online version in colour.)

Figure 9.

Figure 9. The unstable mode of the curvature-controlled model with parameter values in table 3. (a) The mode shape expressed in terms of tangent angle, Inline Formula. (b) The mode shape expressed in terms of displacement Inline Formula The real (solid line) and imaginary (dashed line) part of the mode is shown. The mode corresponds to eigenvalue σ = α + : α = 5.49 s−1, ω/2π = 36.3 Hz (unstable). See electronic supplementary material, movie S4 for corresponding animation.

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 m0 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.

Figure 10.

Figure 10. (a) Frequency ω/2π (Hz) of the least stable mode of the fixed-free, curvature-controlled flagellar model (equation (3.22)) as a function of flagellar length and characteristic dynein bending moment, m0. Other parameters are as in table 3. At each parameter combination (m0, L), frequency is obtained from the imaginary part () of the eigenvalue σ = α + with largest real part (α). (b) Median phase gradient, Inline Formula, of the least stable mode. Anterograde (proximal–distal) propagation corresponds to a phase gradient <0, for ω > 0. For all parameter combinations shown, the least stable mode exhibits anterograde propagation.

5.3. Unstable modes of the geometric clutch model

Figure 11 illustrates the increasing baseline probability of dynein attachment, p0 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, p0, encourages instability. The modes corresponding to the least stable eigenvalues at p0 = 0.05 (with other parameters as in table 4) are shown in figure 12.

Table 4.Parameter values for the geometric clutch model [28]. (L, EI, a, cN are as in table 1.)

parameter value units description
τ 0.05 s dynein time constant
p0 0.05 1 baseline probability of cross-linking
p1 0.15 1 maximum probability of cross-linking
fr 2000 pN μm−1 maximum dynein force per unit length
kT 12.5 pN μm−2 passive shear stiffness
bT 0.25 pN-s μm−2 passive shear friction
Cs 0.50 μm pN-s−1 interdoublet force and dynein coupling factor
k3 30 pN μm−4 nonlinear shear stiffness
Figure 11.

Figure 11. Eigenvalues of the GC model from the weighted residual method: (a) paths of eigenvalues σ = α + in the complex plane as p0 is varied (0 < p0 < 0.2) in the GC flagellar model (equation (3.31)). Other parameters are as in table 3. The red ‘x’ symbols denote the eigenvalues at the final value p0 = 0.2. (b) Frequency ω/2π (Hz) of the least stable mode of the GC model as a function of flagellar length and baseline probability of dynein activation, p0. Other parameters are as in table 4. At each parameter combination (p0, L), frequency is obtained from the imaginary part () of the eigenvalue σ = α + with largest real part (α). (c) Median phase gradient, Inline Formula, of the least stable mode. Anterograde (proximal–distal) propagation corresponds to a phase gradient <0, for ω > 0. For all parameter combinations shown, the least stable mode exhibits anterograde propagation.

Figure 12.

Figure 12. The unstable modes of the GC model (equation (3.31)) with parameter values in table 4. (a,c,e) The mode shape expressed in terms of tangent angle, Inline Formula. (b,d,f) The mode shape expressed in terms of displacement, Inline Formula The real (solid line) and imaginary (dashed line) part of each mode is shown. See electronic supplementary material, movie S5 for corresponding animation. This unstable modes corresponds to eigenvalues σ = α + : (a,b) α = 38.5 s−1, ω/2π = 43.2 Hz. (c,d) α = 7.75 s−1, ω/2π = 32.6 Hz. (e,f) α = 1.20 s−1, ω/2π = 20.3 Hz.

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 k3 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 = 3L/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 (T0) 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 T0). The GC model exhibits propulsive, anterograde propagating waves that resemble the single unstable mode.

Figure 13.

Figure 13. Results from simulation of the full nonlinear equations of each model. (Column 1) Successive snapshots of the flagellar waveform from time-marching simulations; colour shows time increasing from blue (early) to red (later); (column 2, top) time series of angle ψ(L)at s = L; (column 2, bottom) flagellar tension T0 at the base s = 0. Note that a mean value of T0 < 0 indicates axial compression of the flagellum, corresponding to pushing on the base and propulsive behaviour. (Column 3) Plot of ψ(3L/4) versus ψ(L) (clockwise loop = anterograde propagation; counter-clockwise = retrograde); (column 4) fundamental mode from Fourier analysis of simulation. Table 5 lists the relative contribution of each such mode and its correlation with unstable modes of the linearized version. (row 1, a–d, SC 1): sliding-controlled model with sliding at the base (reference [22]; parameters in table 1); oscillation frequency: 27.3 Hz. (Row 2, e–h, SC 2): sliding-controlled model, with no sliding at the base (reference [23]; parameters in table 2); frequency 31.3 Hz. (Row 3, i–l, CC): curvature-controlled model (reference [17]; parameters in table 3); frequency 37.8 Hz. (Row 4, m–p, GC): GC model (reference [28]; parameters in table 4); frequency 46.2 Hz. All models used the same basic physical parameters: length L = 58 µm, flexural rigidity EI = 1700 pN-µm2, diameter a = 185 nm and resistive force coefficient cN = 0.0034 pN-s µm−2.

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).

Table 5.Comparison of fundamental oscillation modes from simulation and unstable modes from eigenanalysis of sliding-controlled, curvature-controlled, and geometric clutch models. Italicized values correspond to the least stable mode.

model simulation: fundamental frequency (Hz) simulation: relative amplitude of fundamental mode eigenanalysis: frequencies of unstable modes (Hz) correlation coefficients: simulation to eigenanalysis
sliding-controlled case 1 27.3 0.881 3.3 0.739
22.5 0.667
30.8 0.272
20.6 0.354
sliding-controlled case 2 31.3 0.972 28.0 0.974
curvature-controlled 37.8 0.929 36.3 0.876
geometric clutch 46.2 0.948 43.2 0.971
32.6 0.804
20.3 0.430

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, cN on eigenvalues of each model.

6. Discussion and conclusion

Unstable and neutrally stable periodic modes were identified for the sliding-controlled [2224], curvature-controlled [17], and doublet separation-controlled (GC) [2527] 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 [2224] 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 [2527]. 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.

Footnotes

References