If the axonemal dyneins were not anchored to the A-tubule, they would walk along the B-tubule towards the base of the axoneme, where the microtubule’s minus end is located. Because the dyneins are anchored to the A-tubule, however, the force instead leads to sliding between adjacent doublets. This inter-doublet sliding in turn causes bending because there are mechanical constraints at the base of the axoneme that resist sliding there. The spatio-temporal coordination of the activity of the dyneins, through a mechanism that is not well understood (see e.g., Sartori et al., 2016b and references therein), gives rise to an approximately sinusoidal bending wave that travels along the axoneme and propels the axoneme through the fluid.

To connect motor activity with bending waves requires an equation of motion. Such an equation was first derived by Kenneth Machin (Machin, 1958), see also (Bayly and Wilson, 2014). It balances active forces against elastic forces (which oppose bending) and hydrodynamic forces (which oppose movement through the fluid). Using this equation, Machin deduced that active forces must be generated all along the flagellum; if motors were only active at the base, like the cracking of a whip, the amplitude would decay rapidly due to the damping from the fluid, and propagating bending waves would not be observed. In other words, “flagellum” (meaning whip in Latin) is a misnomer. Machin’s discovery is especially remarkable because dynein had not yet been discovered (Gibbons and Rowe, 1965) and he did not know that bending is driven by the sliding of the (nearly) incompressible microtubules (Satir, 1968).

Because Machin did not know how the forces were generated, he derived his equation using a continuum, non-molecular approach that is difficult to relate to our current understanding of motor proteins (e.g., (Ishibashi et al., 2020)). In this work, we rederive Machin’s equation by analyzing the forces generated by single dyneins. We show that the equation follows from just three molecular properties of dynein: 1) each active dynein generates a small bending moment, 2) a pair of active dyneins located on opposites sides of the axoneme and at different distances from the base bends the axoneme, 3) a difference in sliding force between two adjacent dyneins on one side of the axoneme produces a normal force that opposes hydrodynamic drag. Summing up these elementary interactions allows us to derive Machin’s equation and therefore deduce where active dyneins are located to produce the observed flagellar bends.

In this paper, we will assume that beating is driven by the outer-arm dyneins, which are anchored every 24 nm along the A-tubule and have two or three force-generating motor domains (depending on the species). This is a simplification as axonemes contain several different classes of inner-arm dyneins in addition to the outer-arm dyneins (Bui et al., 2008). However, there is functional redundancy among the dyneins: mutational studies show that Chlamydomonas cilia are still motile (though they beat more slowly) in the absence of the outer-arm dyneins or when individual classes of inner-arm dyneins are absent (Brokaw and Kamiya, 1987). Thus, given this redundancy, our simplification is likely to reasonable, at least at the level of analysis here.

Because the direction that dynein bends an axoneme depends on the doublet to which it is anchored, we need to use a numbering system for the doublets. The doublets (and associated dyneins) are numbered in Figure 1B according to the convention for sperm (Afzelius, 1988). Doublet microtubule 1 (DM1) is defined as the doublet that lies on the line that bisects the central pair. On the opposite side to DM1 there is usually a bridge that connects doublets five and six, between which the outer-arm dyneins are missing. In the unicellular alga Chlamydomonas reinhardtii, the numbering differs (Hoops and Witman, 1983): the bridge is between Chlamydomonas doublets one and 2 (cDM1 and cDM2), and cDM5 is equivalent to DM1 in sperm. The absence of dyneins between DM5 and DM6 in sperm (and cDM1 and cDM2 in Chlamydomonas) and the presence of the bridge, which presumably impedes sliding, tends to keep the axonemal beat in the x − y plane (see Note 2 in the legend to Figure 1).

To analyze mathematically how the dyneins bend the axoneme, we need to define a coordinate system. In the right-handed coordinate system shown at the bottom of Figures 1A,B, the x -axis is parallel to the axis of the axoneme and the y -axis points into the reverse bend (defined as the bend which has the bridge on the inside Lin et al., 2014); a useful mnemonic is RBI—the Reverse bend has the Bridge on the Inside. An axoneme shape corresponds to a curve y ( x ) . The z -axis points out of the page in Figure 1A and upwards in Figure 1B. This coordinate system defines the sign convention for curvature ( ≅ d 2 y / d x 2 ) and bending moment (see below). The origin of the coordinate system is at the base of the axoneme ( x = 0 , y = 0 ) and the tip of a straight axoneme is at ( x = L , y = 0 ) , where L is the length.

To simplify the analysis, we focus only on the dyneins between doublets three and four (DM3 filled black) and those between eight and nine (DM8 filled gray) (Figure 1C); these are the main drivers for bending in the x - y plane, with the other dyneins generating moments that make smaller contributions to the bend (see Notes 1 and 2 in the legend to Figure 1). The arms extend approximately parallel to the y -axis and are spaced with a period δ = 24 nm along the length of the doublet as shown in Figure 1D, which is in the same orientation as Figure 1A. The dyneins walk towards the minus ends of the microtubules, which are located towards the cell body ( x = 0 ). Therefore, they generate minus-end-directed forces, indicated by the downward arrows. Following this sign convention, the dynein force is − F , where F is positive. Single outer-arm dyneins can generate forces up to 5 pN (Hirakawa et al., 2000). Because the inter-dynein spacing (24 nm) is very small compared to the length of the bends (typical wavelengths are λ ∼ 10  μm), we can define a force density per unit length, − f , where f = F / δ .

The next step is to calculate the bending moments generated by the dyneins. In a beating axoneme, dynein forces vary with position (both along the length and on different sides of the axoneme) and time. It is instructive, however, to start with a simple scenario in which only two dyneins are active and the activity does not change in time. The dyneins are shown in Figure 1E: one at the lower left (gray MTB) and the other at the upper right (black MTB). The lower dynein is anchored to the A-tubule of DM8 and interacts with the B-tubule of DM9 and generates a downward force with a moment arm extending to the left (which is positive in our coordinate system). The magnitude of the moment arm, a ≅ 30 nm, corresponds approximately to the distance between the doublets, though the precise length depends on the molecular structure of dynein and how it generates force. The DM8 dynein generates a small counterclockwise (positive) moment M l = a F l , where − F l is the dynein force on the left (see Note 1 in Figure 1 legend for the definition of the sign of the moment). Note that if the dynein immediately to the right of this dynein (on the opposite side of the midline) were also active and generating a downward force, then the moments would cancel and there would be no net moment (see Note 2 in Figure 1 legend for a more precise statement). This illustrates that a net bending moment requires an imbalance of forces across the axoneme. The dynein anchored to DM3 and interacting with DM4 with its black MTB (Figure 1E, upper right) also generates a downward force, but the moment arm extends to the right (i.e., the negative direction); this dynein generates a small clockwise (negative) moment M r = − a F r . Therefore, in between this pair of dyneins is a region where the moment, M = a F = a ( F l − F r ) , is positive (Figure 1F). This moment bends the intervening axoneme (Figure 1G).

The magnitude of the bend produced by the pair of moments can be calculated using the “beam” equation, M = − κ C , where κ is the flexural rigidity of the axoneme and C = d 2 y / d x 2 is the curvature (Howard 2001). This equation is a consequence of Euler–Bernoulli beam theory and serves as a definition of the flexural rigidity. In accordance with our sign convention, the curvature of the bend in Figure 1G is negative (the angle decreases as x -increases).

The beam equation is analogous to Hooke’s equation for the extension of a spring: F = − k x , where k is the stiffness and x is the extension. To create a bend, two equal and opposite moments are needed (Figure 1H); a single bending moment will cause an object to spin. This is analogous to stretching a spring: two equal and opposite forces are needed (Figure 1I): a single force will cause an object to translate.

If only one dynein is actively generating force, then to get a bend, the moment must be balanced at the base (or tip) of the axoneme, for example by restricting sliding at the basal body or transition zone. If the moment is not balanced at the base (or along the length), then the doublets will slide apart without bending, as observed when the basal restriction to sliding is digested away with proteases (Summers and Gibbons, 1971). In addition to a sliding constraint at the base, bending as shown in Figure 1G requires an additional constraint: DM8 and DM3 must bend together otherwise the DM8-DM9 pair would bend in one direction and DM3-DM4 pair would bend in the other. This constraint is supplied by the radial spokes and additional electrostatic interactions between the doublets that keep the spacing between doublets fixed and maintain the circular cross-section of the axoneme.

A key finding of this analysis is that the location of the bend differs from the location of the active dynein motors that cause it. Figures 1E–G show this: the bend occurs between the active dyneins, and the curvature outside the dyneins is zero.

How big is the bend? Given that the axoneme contains 20 microtubules (9 doublets plus two making up the central pair), we expect that flexural rigidity to be ≥20 times that of a single microtubule (it could be much larger if there is resistance to inter-doublet sliding, which occurs in the absence of ATP and the motors are in rigor, (Howard, 2001)). Therefore, we expect the flexural rigidity to be at least 500 × 10^–24 N∙m² or 500 pN μm² (using the flexural rigidity of a single microtubule in (Howard, 2001)). This flexural rigidity agrees with experimental measurements on intact axonemes (e.g., 800 pN μm² in (Xu et al., 2016)). According to the beam equation, therefore, a single pair of dyneins is expected to generate only a very slight bend with curvature 0.0002 μm⁻¹, corresponding to a radius of curvature of about 5,000 μm ( | C | = a F / κ = 0.03   μ m × 5   p N ÷ 800   p N ∙ μ m 2 ). Thus, the large bends observed in beating axonemes (radius of curvature on the order of 1–10 μm) must be generated by hundreds of dyneins within each wavelength.

To find a general relationship between the distribution of motor forces and the curvature of an axoneme, we derive the static version of Machin’s equation. We start with the beam equation: M ( x ) = − κ C ( x ) . To calculate the total bending moment, M , at position x , we need to add up all the moment densities along the length: M ( x ) = ∫ 0 x m ( x ′ ) d x ′ = ∫ 0 x a ∙ f ( x ′ ) d x ′ , where f = f l − f r is the differential force density across the axoneme ( f l = F l / δ , f r = F r / δ ) . If we assume that the amplitude of the beat is small, then the curvature is approximately C ( x ) ≅ d 2 y / d x 2 . We therefore obtain M ( x ) = ∫ 0 x a ∙ f ( x ′ ) d x ′ = − κ C = − κ d 2 y / d x 2 . Differentiation gives a f ( x ) = − κ d 3 y d x 3 ( x ) (1)

This equation relates the differential dynein force density to the change in curvature in a static (unmoving) axoneme. If the forces are distributed along the axoneme as f ( x ) ∝ sin ( 2 π x ⁄ λ ) as shown in Figure 2A (i.e., f r = f l at x = 0 and DM6-9 maximally active at x = λ /4), then the amplitude can be obtained by triple integration of Eq. (1): y ( x ) ∝ − cos ( 2 π x ⁄ λ ) . The curvature, C = d 2 y / d x 2 = cos ( 2 π x / λ ) , is therefore λ /4 (ninety degrees) out of phase with the motor activity and the principal bend (maximum positive curvature) occurs before the force peaks. The motor forces have maximum amplitude at each end of the bend. This is analogous to the case illustrated in Figures 1E–G. Furthermore, where the curvature is maximum (and minimum) there is no differential motor activity! Thus, motor activity and bending are not co-localized. This is a key finding.

As an axoneme swims, the movement of each increment of length along the flagellum is opposed by viscous forces from the surrounding fluid. In a moving flagellum, therefore, the motor forces must also balance hydrodynamic forces, which are normal to the axis of the axoneme. That the motor moments produce normal forces can be seen with the help of Figure 3. At the ends of each increment of length (length δ in this example), the bending moment generated by the dynein also generates a normal force a F / δ . If each dynein generates the same force, and therefore the same bending moment, then the net normal force is zero. However, if there is a gradient of dynein forces, then the net normal force is non-zero: it is a ∆ F / δ , where ∆ F = F 2 − F 1 . This force can balance the hydrodynamic force acting on the segment: − ξ n v n δ , where v n ≅ d y / d t is the normal velocity of the increment and ξ n is the normal drag coefficient per unit length (Friedrich et al., 2010). Noting that F / δ is the force density, f , and that the change in f , namely ∆ f , occurs over distance ∆ x = δ , we obtain: a ∂ f ∂ x ( x ) = − ξ n ∂ y ∂ t ( x ) (2)

In other words, gradients of motor forces generate normal forces. This equation holds in the limit that the elastic forces are small, in which case the motors only balance hydrodynamic forces. This equation can also be derived by balancing the motor force against all the moments generated by the hydrodynamic forces at locations ≥ x and differentiating (Howard 2001; Appendix 6.2); we have derived it this way to more directly indicate that pairs of dyneins generate the normal force.

In the limit where the hydrodynamic forces dominate, Eq. 2 specifies how a force density, which varies in space and time, determines the shape of the axoneme. For example, let the force be a traveling wave, f ( x , t ) ∝ sin [ ( 2 π ( x ⁄ λ − ν t ) ] where λ is the wavelength and ν is the frequency (in Hz). f ( x , 0 ) is plotted in Figure 2A. The speed is λ ν (traveling from base to tip). In this case, the amplitude y ( x , t ) ∝ sin [ ( 2 π ( x ⁄ λ − ν t ) ] is in phase with the force (Figure 2C). The curvature ( C = − sin [ ( 2 π ( x ⁄ λ − ν t ) ] is therefore 180° out of phase with the motor activity. The sign of the motor activity can be understood by realizing that when hydrodynamic damping dominates, the place with highest velocity needs to be balanced by the motor moments. This place is the straight region where x = λ / 2 (which will become the reverse bend). Thus, motors anchored to DM6-9 need to be active before the straight part ( x = λ / 4 ) and the motors anchored to DM1-4 need to be active after the straight part ( x = 3 λ / 4 ).

That the dyneins anchored to DM6-9 are active in the reverse bend is counter-intuitive. This is because, if only dynein DM8 were active and there was no basal sliding (i.e., the negative moment is at the base and not at DM3 in Figure 1G), then a principal bend, not a reverse bend, would be generated.

In the prior sections we have considered the limiting cases in which motor forces balance only elastic forces (negligible hydrodynamic forces) or they balance only hydrodynamic forces (negligible elastic forces). In general, the motors balance the sum of the hydrodynamic forces and the elastic forces: a ∂ f ∂ x = − ξ n ∂ y ∂ t − κ ∂ 4 y ∂ x 4 (3)

This is Machin’s equation (Eq. 16) in Machin 1958; Machin’s parameter B is the integrated bending moment, M = a ∫ f d x , in our notation). This equation is equivalent to Eq. 16 in Camalet and Jülicher (2000) and Equation (14) Sartori et al. (2016b) if f is replaced by − f (due to the different sign convention used here). If Eq. 3 is differentiated again with respect to x , and we substitute ψ for d y / d x ( ψ is tangent angle), and equate x with arc length s , then we obtain a somewhat more general equation (it holds for small angles ψ and not just small amplitudes y ) used in (Riedel-Kruse et al., 2007). Eq. 3 can also be written m m o t o r s + m d r a g + m e l a s t i c = 0 (4)

showing that the moment densities are balanced.

In summary, we have shown that the Machin’s equation can be derived by summing up the forces generated by pairs of dynein molecules.

Machin’s equation can be used to predict where the force is being generated in a beating axoneme. If the amplitude y ( x , t ) is known in space and time, then we can use Eq. 3 to deduce the force f ( x , t ) . Note that if we have a “motor model”, meaning that we know how the activity of the motors depends on the shape (i.e., f ( y , ∂ y / ∂ x ,   ∂ 2 y / ∂ x 2 … ) ) then Machin’s equation becomes a dynamical system that can integrated (with appropriate boundary conditions) to predict self-organized waveforms. This approach, which entails completing a feedback loop in which motor activity bends the axoneme and the bending of the axoneme feeds back on motor activity, has been used in several earlier studies (see (Sartori et al., 2016b) for references). In this work, we ask the simpler question: given a shape, what is the force profile?

The force profile is readily deduced in the special case where the shape of a beating cilium or flagellum resembles a sinusoid, as is often the case. This so-called travelling wave, with amplitude y ( x , t ) ∝ sin [ ( 2 π ( x ⁄ λ − ν t ) ] , is an approximation that holds in the limit that the axoneme is infinitely long, in which case the boundary conditions can be neglected. The wave travels form base to tip with velocity λ ν . Traveling waves afford a particularly simple relationship between the shape and the motor force: the motor activity is also sinusoidal (seen by substitution into Eq. (3)) with a simple phase shift: f ( x , t ) = sin ⁡ ⁡ [ 2 π ( x ⁄ λ − ν t ) + ϕ ] (5) where the phase is ϕ = arctan ( 1 / M a ) ; (6) M a = ν ξ n   λ 4 ( 2 π ) 3   κ (7)

i s the Machin number, which quantifies the ratio of the viscous forces to the elastic forces for a traveling wave (Geyer et al., 2022). Eq. (6) shows that the Machin number can also be defined from the phase between the force and the amplitude for a traveling wave. ϕ > 0 means that the force leads the amplitude in space and lags the curvature in space. Equivalently, it means that the force leads the curvature in time: the place where the force is high is where the curvature will become high, as noted in Sartori et al. (2016b).

When M a ≪ 1 (which occurs when the wavelength is small, the frequency low, the drag is small, the flexural rigidity is large), the phase is approximately π / 2 . In this case, the curvature leads the motor force in space by ∼ π / 2 (Figure 2B). This is illustrated in (Figure 4B). Equivalently, the curvature lags the motor force in time by ∼ π / 2 : this means that the place where the motors are maximally active is where the principal curvature is increasing and will become the principal bend. When M a ≫ 1 (which occurs when the wavelength is long, the frequency high, the drag is large, the flexural rigidity is small), the phase is shifted towards zero (Figure 4D): now the amplitude is in phase with the force and the curvature lags the motor force (in time) by a phase that approaches π . At an intermediate Machin number, M a = 1 , the curvature lags the force by, 3 π / 4 intermediate between π/2 and π (Figure 4C). The additional phase lag in the presence of drag makes sense because increased damping generally causes an increase in the temporal lag between a response and the force that produces it (for example, a damped spring). Figures 4B,C shows this increasing temporal lag of the curvature behind the force as an increasing spatial lag of the force behind the curvature: the maximum force is further and further behind the curvature (i.e., towards the tip) as M a increases (blue arrows).