Equations of interdoublet separation during flagella motion reveal mechanisms of wave propagation and instability

Abstract

The motion of flagella and cilia arises from the coordinated activity of dynein motor protein molecules arrayed along microtubule doublets that span the length of axoneme (the flagellar cytoskeleton). Dynein activity causes relative sliding between the doublets, which generates propulsive bending of the flagellum. The mechanism of dynein coordination remains incompletely understood, although it has been the focus of many studies, both theoretical and experimental. In one leading hypothesis, known as the geometric clutch (GC) model, local dynein activity is thought to be controlled by interdoublet separation. The GC model has been implemented as a numerical simulation in which the behavior of a discrete set of rigid links in viscous fluid, driven by active elements, was approximated using a simplified time-marching scheme. A continuum mechanical model and associated partial differential equations of the GC model have remained lacking. Such equations would provide insight into the underlying biophysics, enable mathematical analysis of the behavior, and facilitate rigorous comparison to other models. In this article, the equations of motion for the flagellum and its doublets are derived from mechanical equilibrium principles and simple constitutive models. These equations are analyzed to reveal mechanisms of wave propagation and instability in the GC model. With parameter values in the range expected for Chlamydomonas flagella, solutions to the fully nonlinear equations closely resemble observed waveforms. These results support the ability of the GC hypothesis to explain dynein coordination in flagella and provide a mathematical foundation for comparison to other leading models.

Figure 1

(a) Diagram of flagellum with shape defined by tangent angle, ψ(s,t). (b–d) Cross-section, electron micrograph, and functional schematic (rotated 180°). (e) The simplified model is based on two pairs of doublets driving flagellar bending. Side views of P and R doublet pairs show sliding displacement, u(s,t); effective diameter a; active (green); and passive (red) shear forces (f_TP(s,t) and f_TR(s,t)) and internal doublet tension (T_1P, T_2P, T_1R, and T_2R). Schematic cross-sectional views illustrate activity of P and R doublet pairs (note: distortion is exaggerated). P activity drives doublet 4 tipward relative to doublet 2; R activity drives doublet 9 tipward relative to doublet 7. The dashed line segment, formed by connecting doublets 3 and 8, is normal to the beat plane. To see this figure in color, go online.

Figure 2

Free-body diagram of a differential axial element of one doublet pair. Interdoublet components (dynein, nexin links, radial spokes) are represented as a single lumped element contributing net forces f_Tds and f_Nds and net moment af_Tds/2 on each doublet. To see this figure in color, go online.

Figure 3

Schematic diagrams of interdoublet separation. (a) Doublet spacing with zero mean curvature. Interdoublet spacing is a smooth function of axial position, s, modulated by doublet flexural modulus and active and passive axoneme components. Shapes and sizes are exaggerated for illustration. (b) Effect of curvature and doublet tension on interdoublet force. (Blue arrows) Resultants. To see this figure in color, go online.

Figure 4

Simplified model of the effect of interdoublet separation on the rate of cross-bridge attachment or detachment (Eqs. 31 and 32). When the doublets become sufficiently close (h > h_on), attachment probability increases at a characteristic rate, k₀. When h drops below a different threshold (h < h_off). the probability of attachment decreases at a rate that approaches k₀. To see this figure in color, go online.

Figure 5

(a) Null clines and vector field corresponding to the local dynamics of interdoublet separation: Eqs. 34 and 35. (b) Trajectories in the A − h plane rapidly approach a curve representing the quasi-equilibrium value of A for a given value of h: A_eq(h). (c) The system approaches the behavior of a particle on a curved surface with two stable equilibria separated by an unstable equilibrium. To see this figure in color, go online.

Figure 6

Trajectories in (a) dh/dt − h plane and (b) dA/dt − A plane rapidly approach one-dimensional curves as in Fig. 5. These curves resemble cubic polynomial functions, with three zeros representing two stable equilibria and one unstable (threshold) equilibrium. (c) The cubic polynomial analogy to local attachment dynamics: ∂A/∂t ≈ E(A) = C_AA(A – A_th)(1 – A). Derivatives are normalized by τ_N = b_N/k_N. To see this figure in color, go online.

Figure 7

Propagation of interdoublet separation. (a) The decrease in separation h(s,t) is computed from the original PDE of interdoublet separation (Eqs. 33 and 30), and plotted versus s/L at discrete times t = 0.025 n, with n = 1, 2, ...20. (b) Solutions of the simplified excitable system (Eq. 39) plotted versus s/L at the same discrete times. See Movie S1 and Movie S2 in the Supporting Material for corresponding animations. Note the domain is extended to s = 2/L to better visualize propagation. To see this figure in color, go online.

Figure 8

(a) Frequency (imaginary part of eigenvalue: ω = Im(σ)) as a function of baseline attachment probability p₀ and flexural modulus, EI. (b) Region in which the straight position is unstable (R(σ) > 0). At a given value of EI, a Hopf bifurcation occurs (eigenvalues of the GC model cross into the right half-plane) as the baseline probability of attachment is increased. (c and d) Analogous plots of Im(σ) and Re(σ) > 0 in the p₀−L plane. (e and f) The least stable mode of the linearized GC model (either Eq. 46 or 49) with p₀ = 0.10, p₁ = 0.30, C_S = 0.5 μm/pN-s; other parameters are as in Table 4. The real (solid) and imaginary (dashed) part of the mode is shown; the frequency ω = 310 rad/s (49.3 Hz). See Movie S3 and Movie S4 for corresponding animations. To see this figure in color, go online.

Figure 9

Results of time-marching simulation. (a) Asymmetric waveform with baseline probabilities of dynein attachment p_0P = 0.05 and p_0R = 0.01. (b) Symmetric waveform with baseline probabilities of dynein attachment p_0P = p_0R = 0.10. (c–h) Corresponding time series of tangent angle, attachment probability, and deflection at s = L/2. (c, e, and g) Correspondence to the waveform in panel a. (d, f, and h) Correspondence to the waveform in panel b. Except as noted, parameter values are in Table 4. See Movie S5 and Movie S6 for animations.