Generic aspects of axonemal beating

We study the dynamics of an elastic rod-like filament in two dimensions, driven by internally generated forces. This situation is motivated by cilia and flagella which contain an axoneme. These hair-like appendages of many cells are used for swimming and to stir surrounding fluids. Our approach characterizes the general physical mechanisms that govern the behaviour of axonemes and the properties of the bending waves generated by these structures. Starting from the dynamic equations of a filament pair in the presence of internal forces we use a perturbative approach to systematically calculate filament shapes and the tension profile. We show that periodic filament motion can be generated by a self-organization of elastic filaments and internal active elements, such as molecular motors, via a dynamic instability termed Hopf bifurcation. Close to this instability, the behaviour of the system is shown to be independent of many microscopic details of the active system and only depends on phenomenological parameters such as the bending rigidity, the external viscosity and the filament length. Using a two-state model for molecular motors as an active system, we calculate the selected oscillation frequency at the bifurcation point and show that a large frequency range is accessible by varying the axonemal length between 1 and 50 µm. We discuss the effects of the boundary conditions and externally applied forces on the axonemal wave forms and calculate the swimming velocity for the case of free boundary conditions.


Introduction
Many small organisms and cells swim in a viscous environment using the active motion of cilia and flagella. These hair-like appendages of the cell can undergo periodic motion and use hydrodynamic friction to induce cellular self-propulsion [1]. Of particular interest for the present work are those flagella and cilia which contain an axoneme, a characteristic structure which occurs in a large number of very different organisms and cells and which appeared early in evolution. It consists of a cylindrical arrangement of nine doublets of parallel microtubules, which are elastic rod-like protein filaments, and one pair of microtubules in the centre. In addition, it contains a large number of other proteins, such as nexin, which provide elastic links between the microtubule doublets, see figure 1. The axoneme is inherently active. A large number of dynein molecular motors are located in two rows between neighbouring microtubules. In the presence of ATP, which is a chemical fuel, these motors can generate forces and thus induce local displacements between adjacent microtubules [2]. Axonemal cilia and flagella thus represent rod-like elastic structures which move and bend as a result of stresses that are generated internally. They are used by single-celled organisms such as paramecium, which has a large number of cilia on its surface, and chlamydomonas, which uses two flagella to swim [1]. The tail of sperm consists of a single flagellum used for swimming. Cilia and flagella can generate periodic waving or beating patterns of motion. In the case of sperm, for example, a bending wave propagates from the head towards the tail. In a viscous environment, the surrounding fluid is set in motion and hydrodynamic forces act on the filament. For typical values of frequencies and length scales given by the size of a flagellum, the Reynolds number is small and the inertia terms in the fluid hydrodynamics can be neglected [3]. Therefore, only friction forces resulting from solvent viscosity contribute to propulsion. Under such conditions self-propulsion is possible if a wave propagates towards one end, a situation which breaks the time-reversal invariance of the sequence of deformations of the flagellum [4].
How can bending waves be generated by the action of molecular motors within the axoneme? Dyneins induce relative forces between neighbouring microtubules. As a result of these forces, the filaments have the tendency to slide with respect to each other. If such a sliding is permitted globally, the filaments simply separate but no bending occurs. Bending results if the global sliding is suppressed by rigidly connecting the filament pair in the region close to one of the two ends. In this situation, sliding is still possible locally, however, but only if the filaments undergo a bending deformation. This coupling of axonemal bending to local microtubule sliding has been demonstrated experimentally. In situations where the axoneme is cut at its basal end, filaments slide and in the presence of ATP separate without bending [5]. Small gold-beads specifically attached to microtubules in fully functioning flagella can be used to directly visualize the local relative sliding during beating [6].
The dynamics of axonemal bending and wave patterns have been addressed theoretically by several authors [7]- [17]. One can distinguish two principally different mechanisms to generate oscillatory deformation patterns of the axoneme. (i) Deterministic forcing: a chemical oscillator could regulate the dynein motors to be activated and deactivated periodically. In this case, the internal active system creates a dynamical force pattern, which drives the system in a deterministic way [8,9]. (ii) Self-organized beating: the axoneme oscillates spontaneously as a result of the interplay of force-generating elements and the elastic filaments. Machin [10] and Brokaw [11] suggested the possibility of self-organized beating assuming a force-generating system capable of oscillations. Patterns of beating have been studied using numerical simulations of a variety of different models [10]- [16]. The effect of a coupling of the motor activity to the sliding displacement as well as to curvature were discussed in order to find specific situations where simulated and observed bending waves could be matched qualitatively and quantitatively.
Our present work approaches the problem differently. Instead of looking for specific models we are interested in characterizing the general properties of the class of systems which we call internally driven filaments [17] and for which the axoneme is an example. Most important, we focus on an oscillating instability of the motor-filament system. In general, spontaneous oscillations occur via a so-called Hopf bifurcation [18], where an initially stable quiescent state becomes unstable and starts to oscillate. In the vicinity of this bifurcation, the behaviour of the system is governed by linear terms which are generic and do not depend on the microscopic details. We thus can characterize the general behaviours of the system which follow from very few assumptions. Furthermore, we can discuss what types of nonlinearities occur as the oscillation amplitude increases. Far from the bifurcation the behaviour is not generic and is therefore difficult to predict since it could depend on minor structural details and could vary from one type of cell to another.
There is evidence for the existence of a dynamic oscillatory instability in the axoneme. Demembranated flagella show a behaviour which depends on the ATP-concentration C ATP in the solution. For small C ATP , the flagellum is straight and does not move. If C ATP is increased, oscillatory motion sets in at a critical value of the ATP concentration and persists for larger concentrations [19]. This implies an instability of the initial straight state with respect to an oscillating mode. Recently, it has been demonstrated, using a simple model for molecular motors, that a large number of motors working against an elastic element can generate oscillations via a Hopf bifurcation by a generic mechanism [20,21]. The fact that motors operate in large numbers allows the system to generate noiseless, phase coherent oscillations using motor molecules as active elements which operate stochastically. This suggests that a system consisting only of molecular motors and semiflexible filaments can in general generate self-organized oscillations. This idea is supported by the facts that flagellar dyneins are capable of generating oscillatory motion [22] and that experiments suggest the existence of dynamic transitions in many-motor systems [23,24]. The patterns of the motion of cilia and flagella can be complex and embedded in three-dimensional space. Many examples of propagating bending waves, however, are planar and can thus be considered as confined to two dimensions [6]. In the following sections, we present a systematic study of the dynamics of elastic filaments subject to internal forces [17]. In section 2, we derive the dynamic equations. Similar equations have been proposed by [12,13]. Our derivation follows recent work on the dynamics of semiflexible filaments subject to external forces [25]- [29], which we generalize for internal forces. We show how the dynamic equations can be solved perturbatively and we determine the shapes of bending waves, discuss the velocity of swimming and calculate the tension profile along the flagellum. In section 3, we discuss the general properties of an active system that generates internal forces. The behaviour of self-organized bending waves and oscillations in the vicinity of a Hopf bifurcation of the active system coupled to the filament dynamics are discussed in section 4. We calculate the generic wave patterns at the bifurcation and determine the selected oscillation frequency using a two-state model for active elements. Finally, we propose experiments which could be performed to test predictions that follow from our work.

Internally driven filament
The cylindrical arrangement of microtubule doublets within the axoneme can be described effectively as an elastic rod. Deformations of this rod lead to local sliding displacements of neighbouring microtubules. Here, we consider planar deformations. In this case, the geometric coupling of bending and sliding can be captured by considering two elastic filaments (corresponding to two microtubule doublets) arranged in parallel with a constant separation a along the whole length of the rod. At one end, which corresponds to the basal end of an axoneme and which we call the 'head', the two filaments are rigidly attached and are not permitted to 24.5 slide with respect to each other [2]. Everywhere else, sliding is possible, see figure 2. The configurations of the system are described by the shape of the filament pair given by the position of the neutral line r(s) as a function of the arclength s, where r is a point in two-dimensional space. The shapes of the two filaments are then given by where n, with n 2 = 1, is the filament normal. The local geometry along the filament pair is characterized by the relationṡ where t denotes the normalized tangent vector and C =ṫ·n is the local curvature. Throughout this paper the overdots denote derivatives with respect to s, i.e.ṙ ≡ ∂ s r ≡ ∂r/∂s.

Bending and sliding of a filament pair
The bending of the filament pair and local sliding displacements are coupled by a geometric constraint. The sliding displacement at position s along the neutral line is defined to be the difference of the total arclengths along the two filaments up to the points r 1 (s) and r 2 (s) which face each other along the neutral line. From equations (1) and (4) follows thatṙ 1,2 = (1 ± aC/2)t and thus Therefore, the sliding displacement is given by the integrated curvature along the filament.

Enthalpy functional
The bending elasticity of filaments (microtubules) characterizes the energetics of the filament pair. In addition to filament bending, we also have to take into account the internal stresses due to active and passive elements. We introduce a coarse-grained description, defining the force per unit length f (s) acting at position s in opposite directions on the two microtubules, see figure 2 (the internal forces must balance). This internal force density corresponds to a shear stress within the flagellum which tends to slide the two filaments with respect to each other. A static configuration of a filament pair of length L can thus be characterized by the enthalpy functional Here, κ denotes the total bending rigidity of the filaments. The incompressibility of the system is taken into account by the Lagrange multiplier function Λ(s) which is used to enforce the 24.6 constraintṙ 2 = 1 which ensures that s is the arclength [25]. The internal force density f couples to the sliding displacement ∆ as described by the contribution ds f∆. The variation of this term under small deformations of the filament shape represents the work performed by the internal stresses. Using equation (7) leads, after a partial integration, to where is the force density integrated to the tail. From (9) it follows that if the internal stresses or F are imposed, G is minimized for a filament curvature C = C 0 where C 0 (s) = aF (s)/κ is a local spontaneous curvature. The internal forces therefore induce filament bending. In order to derive the filament dynamics, we determine the variation δG with respect to variations δr. Details of this calculation are given in appendix A. As a result, we find δG δr Here, plays the role of the physical tension. This becomes apparent since from equation (11) it follows that where F ext (L) is the external force applied at the end which satisfies (20). Therefore, τ is the tangent component of the integrated forces acting on the filament.

Dynamic equations
We assume, for simplicity, that the hydrodynamics of the surrounding fluid can be described by two local friction coefficients ξ and ξ ⊥ for tangential and normal motion, respectively. The equations of motion in this case are given by [26,27] For the following, it is useful to introduce a coordinate system r = (X, Y ) and the angle ψ between the tangent t = (cos ψ, sin ψ) and the X-axis which satisfies C =ψ. We find with equation (11) that Noting that ∂ tṙ = n∂ t ψ, we obtain an equation of motion for ψ(s) alone:  Figure 3. Resting frame (X, Y ) and frame (x, y) moving with the filament at velocityv. The tangent t and the normal n are indicated, the angle between the X, x-axis and t is denoted ψ.
The tension τ is determined by the constraint of incompressibility ∂ tṙ 2 = 2t·∂ tṙ = 0. This condition and (15) leads to a differential equation for the tension profile: Equations (16) and (17) determine the filament dynamics [12]. The filament shape follows from where r(0, t) can be obtained from (15) evaluated at s = 0.

Boundary conditions
The filament motion depends on the imposed boundary conditions. The variation δG has contributions at the boundaries which can be interpreted as externally applied forces F ext and torques T ext acting at the ends, see appendix A. At the head with s = 0, Similarly, at the tail for s = L, 24.8 If constraints on the positions and/or angles are imposed at the ends, forces and torques have to be applied to satisfy these constraints. In the following, we will discuss different boundary conditions as specified in table 1. Case A is the situation of a filament with a clamped head, i.e. both the tangent and the position at the head are fixed, the tail is free. Case B is a filament with fixed head, i.e. the tangent at the head can vary. The situation of a swimming sperm corresponds to case C where the friction force of a viscous load attached at the head is taken into account, the head is otherwise free. As an example of a situation with an external force F ext (L) applied at the tail we consider case D. For simplicity, we assume a force parallel to the (fixed) tangent at the head.

Small deformations
In the absence of internal forces f , the filament relaxes passively to a straight rod withψ = 0, i.e. ψ is a constant. Without loss of generality, we choose ψ = 0 in this state, which implies that the filament is parallel to the X-axis. Internal stresses f (s) induce deformations of this conformation. We can perform a systematic expansion of the filament dynamics in powers of the dimensionless stress amplitude , where f (s, t) = f 1 (s, t) and f 1 characterizes the stress distribution along the filament. The dynamic equations (16) and (17) can be solved perturbatively by writing which allows us to determine the coefficients ψ n (s, t) and τ n (s, t) as described in appendix B. The coefficients ψ 2 and τ 1 vanish by symmetry. Indeed, under the transformation f → −f , we have ψ → −ψ and τ → τ . In the cases considered here, the coefficient τ 0 = σ is a constant, which vanishes for boundary conditions A-C. In case D, it is equal to the X-component of the external force applied at the end. To second order in , the filament shape is thus characterized by the behaviour of ψ 1 (s, t) which obeys In order to discuss the filament motion in space, we define the average velocity of swimminḡ v = lim which is independent of s. This velocity is different from zero only in the case of boundary condition C. We choose the X-axis parallel tov and introduce a coordinate system (x, y) which moves with the filament wherev = ±|v| and the sign depends on whether the motion is towards the positive or negative X-direction, see figure 3. It is useful to introduce the transverse and longitudinal deformations h and u, respectively, which satisfy and which vanish for = 0. The quantities h, u andv can be calculated perturbatively in , see appendix B. To second order in , transverse motion satisfies The boundary conditions for h(s) are given in table 2.

Swimming
The velocity of swimmingv can be calculated perturbatively in . We consider, for simplicity, periodic internal stresses with frequency ω given by † f (s, t) =f (s) iωt +f * (s) e −iωt (28) which after a long time leads to periodic filament motion h(s, t) =h(s) e iωt +h * (s) e −iωt (29) whereh satisfies, according to equation (26), Writingh = H e −iφ where H and φ denote the amplitude and phase, respectively, the filament motion can be expressed as which represents bending waves for which v p = ω/φ can be interpreted as the local wave propagation velocity. For a freely oscillating filament with a viscous load of friction coefficient ζ attached at the head (case C), we find an average propulsion velocityv is the velocity for ζ = 0. Equation (32), which is correct to second order in , confirms that motion is only possible if the filament friction is anisotropic ξ /ξ ⊥ = 1 and if a wave is propagating, φ = 0, consistent with earlier work on swimming [3,4,30,31,32]. For a rod-like filament with ξ < ξ ⊥ , the swimming motion is opposite to the direction of wave propagation. † Note, that more general periodic stresses f = nf n e inωt can also be considered. Since the equation of motion (26) is linear, different modes superimpose linearly and we can without loss of generality restrict our discussion to a single mode n = 1.

Internal forces
The filament dynamics are driven by internal stresses f (s, t) generated by a large number of active and passive elements. The general aspects of the active properties of a large number of force-generating elements can be studied using the simple two-state model [20,21], reviewed in appendix C.

Active elements coupled to a filament pair
We assume a coarse-grained description where at any position s an independent active system containing many force-generating elements is located, which generates the sliding displacements ∆(s) as well as the internal force density f (s). Note that we neglect fluctuations which arise from the chemical activity of a finite number of force generators and we assume homogeneity of all properties along the filament length.
Since we will study periodic motion, we express forces and displacements by a Fourier expansion The force f (s, t) induced by the active system located at position s along the filaments does in general depend on the full history of the sliding motion ∆(s, t ) which the active elements have experienced (t ≤ t). For limit cycle motion, this can be expressed by a nonlinear relation of the Fourier coefficients according to the expansion [21] f n = F (1) Time-translation invariance requires that the coefficients F (m) n,k 1 ,...,km = 0 if k 1 +· · ·+k m = n [21]. The derivation of the coefficients F (m) for the case of a two-state model is given in appendix C. Note that (35) is general and characterizes a whole class of active nonlinear mechanical systems. It is not restricted to a two-state model.

Symmetry considerations
Consider the case where the active elements are connected to one of the filaments and generate the sliding of the second. Since the filaments are polar structures, filament sliding is typically induced towards one particular end of the filaments. As a consequence, such a system has a natural tendency to create bending deformations with one particular sign of the curvature. Exchanging the role of the two filaments reverses this sign. The system therefore has a broken symmetry with respect to bending deformations with positive and negative curvature.
This broken symmetry does not however reflect the symmetry of the axoneme, which is symmetric with respect to all microtubule doublets. Each of the nine microtubule doublets within the cylindrical arrangement of the axoneme play identical roles and interact with rigidly attached motors as well as with motors which slide along their surfaces. If all dynein motors generate a constant force, there is consequently no resulting bending deformation since all bending moments cancel by symmetry. In order to introduce this axonemal symmetry in the active filament pair, we assume that both filaments play identical roles, see symmetry is that the filament polarity becomes irrelevant. This symmetry does not usually exist in situations where molecular motors occur because of the filament polarity which defines a direction of motor motion. In the following, we assume a symmetric motor-filament pair, which implies effectively a symmetric force-generating system. For this symmetry, even numbered coefficients F (2n) vanish in the expansion given in equation (35).

Linear response function
The perturbative treatment for small filament deformations introduced in section 2 can be applied to the coupled motor-filament system using the expansion (35). Up to second order in , only the linear response coefficient of the active system plays a role. In the following, we use F which is the linear response function obtained for a two-state model, see appendix C. Here, K is an elastic modulus per motor, λ an internal friction coefficient per motor, k is the cross-bridge elasticity of a motor and Ω, with 0 < Ω < π 2 , plays the role of a control parameter, α is a characteristic ATP cycling rate. Higher-order terms F (2n+1) have to be taken into account if the third or higher order in is considered.

Self-organized beating via a Hopf bifurcation
A Hopf bifurcation is an oscillating instability of an initial non-oscillating state which occurs for a critical value Ω c of a control parameter. For Ω < Ω c , the system is passive and not moving, while for Ω > Ω c it exhibits spontaneous oscillations. As we demonstrate below, self-organized oscillations of the driven filament pair are a natural consequence of its physical properties. The control parameter Ω in our description is related to the ratio of two chemical rates of the ATP hydrolysis cycle. In an experimental situation, it could be varied, for example, by changing the ATP concentration. If the molecular motors are regulated by some ion concentration such as, for example, Ca 2+ , this concentration could also play the role of the control parameter.

Generic aspects
For oscillations in the vicinity of a Hopf bifurcation, the filament deformations are small and we can thus use (26). The force amplitudef can be obtained from the expansioñ in the amplitude of sliding displacements. In addition, this amplitude is related to the filament shape: With equation (30), we find that spontaneously oscillating modesh(s) at frequency ω near the bifurcation are solutions to A B C Figure 6. Filament shapes in the (x, y) plane that correspond to the modes displayed in figure 5. Snapshots taken at different times to illustrate the bending waves for boundary conditions A-C. The arrows indicate the direction of wave propagation.
Note that this equation is general and its structure does not depend on the specific model chosen for the force-generating elements. Boundary conditions corresponding to cases A-D follow from table 2 by replacing h →h, ∂ t h → iωh and f →f = aχ(∂ sh − ∂ sh (0)). Equation (39) has the structure of an eigenvalue problem, with χ playing the role of a complex eigenvalue. As discussed in appendix D, it is convenient to introduce the dimensionless tensionσ ≡ σL 2 /κ, frequencyω ≡ ωξ ⊥ L 4 /κ, linear response functionχ = L 2 a 2 χ/κ and s = s/L. For every pair ofσ andω there exists, for given boundary conditions, an infinite set of non-trivial solutionsh n (s) to (39) ifχ is equal to a corresponding dimensionless eigenvalueχ n , n = 1, 2, . . . . Therefore, the complex eigenvalues of equation (39) can be expressed in the form We order these values according to their amplitude |χ n | < |χ n+1 |. The existence of discrete eigenvalues reveals that only particular values of the linear response function of the active material are consistent with being in the vicinity of a Hopf bifurcation. If the actual value of χ(Ω, ω) of the system differs from one of the eigenvalues χ n , the system is either not oscillating, or it oscillates with large amplitude for which higher-order terms in cannot be neglected. An important consequence of this observation is that the modesh n of filament beating close to a Hopf bifurcation can be calculated without knowledge of the properties of the active elements. Not even knowledge of the linear response coefficient is required since it follows as an eigenvalue at the bifurcation. Therefore, filament motion close to a Hopf bifurcation has generic or universal features which do not depend on the details of the problem, such as the structural complexity within the axoneme. The motion is given by one of the modesh n which are functions of s/L that only depend on two dimensionless parametersσ andω. Figure 5 displays the amplitude H and the gradient of the phaseφ of the fundamental modeh 1 = H e −iφ . Note, that we use arbitrary units for H since the amplitude ofh(s) is not determined by (39). Corresponding filament shapes represented according to (25) in the (x, y)

Selection of eigenmodes and frequency
In order to determine which of the modesh n is selected and at what frequency ω the system oscillates at the bifurcation point, the functional dependence of the linear response function χ(Ω, ω) on frequency ω and the control parameter Ω is crucial. Without taking into account aspects of the more microscopic processes within the active material, the selected frequency cannot, in general, be predicted. The linear response function calculated by the generic twostate model is a natural choice to discuss the selected frequency for a situation involving many molecular motors. We therefore use the expression for χ given by (36) in order to determine the selected frequency. For Ω = 0, χ = K + iλω, which is a passive viscoelastic response. In this case, no spontaneous motion is possible, which can be seen by the fact that all eigenvalues χ n have negative real and imaginary parts which cannot be matched by the linear response of a passive system. If Ω is increased, we are interested in a critical point Ω = Ω c for which the straight filament configuration becomes unstable. This happens as soon as the linear response function χ matches, for a particular frequency ω = ω c , one of the complex eigenvalues, κ a 2 L 2χ n (σ,ω c ) = χ(Ω c , ω c ) whereω c = ξ ⊥ L 4 ω c /κ. Since χ is complex, both the real and the imaginary parts of (41) represent independent conditions. Therefore, (41) determines both the critical point Ω c as well as the selected frequency ω c . The selected mode n is the first mode (beginning from Ω = 0) to satisfy this condition. Since |χ(Ω, ω)| increases with increasing Ω, the instability occurs almost exclusively for n = 1 since |χ 1 | has been defined to have the smallest value. For Ω > Ω c , but close to the transition, the system oscillates spontaneously with frequency ω c . The shapes of filament beating are characterized by the corresponding modeh n (s) whose amplitude is determined by nonlinear terms. For the case of a continuous transition, the deformation amplitude increases as |h n | ∼ (Ω − Ω c ) 1/2 .

Axonemal vibrations for different lengths
Choosing parameter values which correspond to the axonemal structure, we estimate the bending rigidity κ 4 × 10 −22 N m 2 which is the rigidity of about 20 microtubules. Furthermore, we use a 20 nm, a motor density ρ 5 × 10 8 m −1 , friction per unit length λ 1 kg (m s) −1 , a rate constant α 10 3 s −1 , a cross-bridge elasticity k 10 −3 N m −1 and a perpendicular friction ξ ⊥ 10 −3 kg (m s) −1 , which is of the order of the viscosity of water. A rough estimate of the elastic modulus per unit length associated with filament sliding can be obtained by comparing the number of dynein heads to the number of nexin links within the axoneme. This suggests that K is relatively small, K < ∼ kρ/10. The selected frequency ω c at the bifurcation point for case B is shown in figure 8 as a function of the axoneme length for K = 0. For small lengths, the oscillation frequency is large and increases as L decreases. In this high-frequency regime, which occurs for L < ∼ 10 µm, the system vibrates in a mode with no apparent wave propagation. For L > ∼ 10 µm the frequency is only weakly L-dependent and the system propagates bending waves of a wavelength shorter than the filament length. This length dependence of the frequency is qualitatively consistent with experimental results [33].
The limit of small lengths corresponds to smallω and can be studied analytically. For σ = 0, the eigenvalue χ 1 is given to linear order inω byχ 1 −π 2 /4 − iγω. The coefficient γ depends on the boundary conditions, but not on any model parameters. For a clamped head (case A), γ 0.184, for a fixed head (case B), γ 0.008, see appendix D. The criterion for a Hopf bifurcation for small L is

24.16
Together with equation (36), we find the critical frequency where L K ≡ (π 2 κ/4Ka 2 ) 1/2 and L λ ≡ (λa 2 /γξ ⊥ ) 1/2 are two characteristic lengths. For K and λ small, L K L L λ and the critical frequency behaves as ω c ∼ 1/L 2 . The critical value of the control parameter for small L is given by For axoneme lengths below a characteristic value, L < L min , the condition Ω c ≤ π 2 is violated. Therefore, oscillations exist only for L > L min . If (K + λα)/ρk 1, L min (κ/4a 2 ρk) 1/2 1 µm for our choice of the parameters. The maximal frequency of oscillations occurs for L = L min . Within the range of lengths between 1 and 50 µm, we find frequencies ω c /2π which vary between tens of hertz and about 20 kHz. The parameters chosen in figure 8 lead to frequencies above 150 Hz. Lower frequencies are obtained for larger values of the friction, λ, or smaller chemical rate α.

Effect of external forces applied at the tail
In the presence of an external force applied at the endσ = 0 and the filament pair is under tension. The eigenvaluesχ n depend onσ, therefore the tension affects the shape and frequency of an unstable mode. Inspection of (39) suggests that the tension σ can be interpreted as an additive contribution to the linear response coefficient χ → χ + σ/a 2 . Taking into account the boundary conditions complicates the situation, theσ dependence of the eigenvalues is in general non-trivial. However, for a filament with a clamped head and a force applied at the tail parallel to the tangent at the head (case D) the effect of tension can be easily studied. In this particular case, the tension leads to the same contribution to χ both in (39) and in the boundary conditions; the eigenvalues thus depend linearly on σ: Since the tension corresponds to a change of the real part of χ, its presence has the same effect on filament motion as an increase of the elastic modulus K per unit length in equation (36) by σ/a 2 . Therefore, we can include the effects of tension on the critical frequency in (43) by replacing K → K + σ/a 2 , which reveals that in those regimes where the frequency depends on K, it will increase for increasing tension. The critical value of the control parameter Ω c is a function of the tension; with equation (44), we find For our choice of parameters, ρ Ka 2 200 pN which indicates that forces of the order of 10 2 pN should have a significant effect on the bifurcation. Furthermore, equation (46) indicates that the tension σ can play the role of a second control parameter for the bifurcation. Consider a tensionless filament oscillating close to the bifurcation for Ω > Ω c (0). If a tension is applied, the critical value Ω c (σ) increases until the system reaches for a critical tension σ c a bifurcation point with Ω c (σ c ) = Ω and the system stops oscillating. 24.17

Discussion
We have studied a class of physical systems which we call internally driven filaments [17]. These systems are elastic rod-like filaments which contain active force-generating elements that generate bending deformations and periodic motion. Realizations of these systems exist in living cells in the form of cilia and flagella which contain microtubules and dynein molecular motors. We show that these systems, in general, can undergo a Hopf bifurcation for a critical value of a control parameter, which could in the case of biological systems be, for example, the ATP concentration. At this bifurcation an inert state becomes unstable with respect to the oscillatory behaviour in the form of beating motion and propagating bending waves. Close to the Hopf bifurcation the behaviour is completely determined by linear terms, which are generic and do not depend on the structural details of the system. Therefore, even though cilia and flagella are complex and consist of about one hundred different types of proteins, oscillatory motion and bending waves can be understood by a simple coarse-grained description based on effective material properties such as the bending rigidity, the oscillation frequency and the viscosity of the surrounding fluid.
Our results have been obtained using several simplifying assumptions. We have neglected hydrodynamic interactions between different parts of the filament which lead to logarithmic corrections. It is straightforward to generalize our approach to include the effects of hydrodynamic interaction; however the qualitative results remain the same. A more important simplification is the restriction to a two-dimensional system and to filament configurations which are planar. This is motivated by the fact that the observed bending waves of flagella are planar in many cases [6]. It has been speculated that structural features in the axoneme are responsible for the two-dimensional nature of the beating. For example, there is evidence for permanent bonds between the two microtubule doublets in some flagella which would prohibit interdoublet sliding and could confine bending deformations in one plane. Similarly, there could be connections linking microtubule doublets to the central pair via radial spokes; see, for example, [34] for a review. In order to address the question under what conditions beating is planar and to study nonplanar modes of motion, three-dimensional generalizations of our model could be used. Such generalizations introduce additional phenomena. In particular, torsional deformations become relevant. The local sliding displacements depend, in the three-dimensional case, both on the bending and torsional deformations, and the full torsional dynamics have to be accounted for in the dynamic equations.
In this work we have performed a systematic expansion for small deformations. Using this approach, the filament dynamics can be studied in some limits analytically. This expansion allows us to distinguish and clearly characterize linear and nonlinear terms. For filament beating with larger amplitude, nonlinearities become relevant. Such nonlinearities arise due to nonlinear geometric terms, non-Hookian corrections to the elasticity and as nonlinearities in the force-generation process of active elements. They could give rise to new types of behaviour such as additional dynamic instabilities. In contrast to the linear terms, the selection of the dominant nonlinearity can depend on structural details. Therefore, large amplitude motion is much more difficult to describe and to analyse since a knowledge of the dominating nonlinearities is necessary. However, if no subsequent instabilities occur, the principal effect of the nonlinear terms is to determine the amplitude of the propagating waves. Therefore, in many cases propagating waves with larger amplitudes should be well approximated by the linearly unstable modes.
Experiments have so far focused on the nonlinear regime with large amplitude motion, which is the regime of ordinary swimming flagella. Our results show that it would be important to perform experiments near the bifurcation point where simple predictions of our approach could be tested. Such experiments could be, for example, performed by using demembranated flagella, and ATP concentrations not far from the critical level where beating sets in. If the observed waves do not correspond to the modes obtained here, the physical mechanism underlying the oscillations must be different, for example a coupling of active elements to filament curvature might be involved. Our predictions could be tested experimentally by micromanipulation experiments where external forces are applied to individual beating flagella. Such experiments have recently become possible [35]. In particular, the influence of boundary conditions would be interesting to study since we find that boundary conditions and external forces should have a significant influence on the bending waves.
The frequency selected by the system at the Hopf bifurcation depends on the characteristic time scale of the chemical cycle of the active elements. Our calculation of the selected frequency as a function of the axonemal length shows that the physical mechanism that generates these oscillations can cover a large frequency range from several tens of hertz up to about 20 kHz. This finding is important since it shows that the axonemal structure can be used to build very naturally high-frequency mechanical oscillators using ordinary molecular motors with characteristic time scales of the order of milliseconds. Mechanical oscillators are of particular importance in the context of auditory cells, which amplify and detect faint sounds at one particular frequency. There is much evidence that these cells operate at a Hopf bifurcation of a dynamic mechanical system. The precise nature of this oscillator has, so far, been evasive and there is evidence that it could be located in the hair bundles of auditory cells [36]. A close look at the hair bundles of non-mammalian vertebrates reveals the existence of a single axonemal structure, called the kinocilium [37]. Its role is not understood, however our work shows that it is well suited to play the role of an active oscillating element, being able to cover the full frequency range of auditory cells [38].

Acknowledgments
We thank A Ajdari, M Bornens, H Delacroix, T Duke, R Everaers, K Kruse, A Maggs, A Parmeggiani, M Piel, J Prost and C Wiggins for useful discussions.

Appendix A. Functional derivative of the enthalpy
The variation δG of the enthalpy G given by (8)

24.19
The functional derivative δG/δr given by (11) can be read off from the integrand, the boundary terms provide expressions for forces and torques at the ends given by (19) and (20).

Appendix B. Small deformations
Inserting the expansions (21) in the differential equation (17) for the tension profile leads at each order in to a separate equation. With ψ 2 = 0 and τ 1 = 0 by symmetry, up to second order, we findτ From the boundary conditions, it follows that τ 0 = σ. Repeating the same procedure for the dynamic equation (16) using τ 0 = σ, we obtain The transverse and longitudinal displacements h and u and the velocityv are obtained perturbatively as Using equations (21) and (18) we find u 1 (s) = u 1 (0), h 2 (s) = h 2 (0) and The dynamics of h(s) and u(s) to second order in are therefore determined by ψ 1 (s, t) and the motion ∂ t r for s = 0. For the latter, we find from (15) thatv 1 = 0, ∂ t u 1 (0) = 0, ∂ t h 2 (0) = 0 and For the lateral motion, we obtain with equations (B2) and (B4), Using equation (B6), we obtain which is the result given in equation (26). In order to calculate τ 2 (s) andv 2 we integrate (B1). Together with (B9), we obtaiṅ

Appendix C. Two-state model for molecular motors
We review the two-state model for the generic behaviours of coupled motor molecules [20]. Each motor has two different chemical states, a strongly bound state, 1, and a weakly bound state, 2. The interaction between a motor and a filament in both states is characterized by energy landscapes W 1 (x) and W 2 (x), where x denotes the position of a motor along the filament. The potentials have the filament symmetry: they are periodic with period l, W i (x) = W i (x + l) and are spatially asymmetric, W i (x) = W i (−x). In the presence of ATP, the motors undergo transitions between states with transition rates ω 1 and ω 2 . Introducing the relative position ξ of a motor with respect to the potential period, where x = ξ + nl, 0 ≤ ξ < l and n is an integer, we define the probability P i (ξ, t) for a motor to be in state i at position ξ at time t, P 1 + P 2 is normalized within one period. The dynamic equations are [20] ∂ t P 1 + v∂ ξ P 1 = −ω 1 P 1 + ω 2 P 2 ∂ t P 2 + v∂ ξ P 2 = ω 1 P 1 − ω 2 P 2 .
The sliding velocity v = ∂ t ∆, where ∆ is the displacement, is determined by the force balance The coefficient λ describes the total friction per unit length in the system, ρ is the number density of motors along the filament and f is the force per unit length generated by the system. The elastic modulus per unit length K occurs in the presence of elastic elements such as nexins. For an incommensurate arrangement of motors as compared to the filament periodicity, P 1 + P 2 = 1/l and the equations can be expressed by P 1 alone ∂ t P 1 + v∂ ξ P 1 = −(ω 1 + ω 2 )P 1 + ω 2 /l f = λv + ρ l 0 P 1 ∂ ξ ∆W + K∆ (C3)
In order to determine the amplitudes A i the four boundary conditions have to be used. Since (39) is linear and homogeneous, the equations for the coefficients A i have the form where M ij is a (4 × 4) matrix which depends onχ,σ andω. Non-trivial solutions exist only for those valuesχ =χ n for which det M ij = 0. The corresponding eigenmode is the solution for A i of (D3). These eigenvalues and eigenfunctions can be determined by numerically obtaining solutions for det M ij = 0. In the limit of smallω, analytic expressions can be obtained by inserting the ansatzχ =χ (0) +χ (1)ω + O(ω 2 ) for the eigenvalue which leads to an expansion of det M ij in powers ofω 1/2 . This procedure leads forσ = 0 toχ (0) n = −[(n − 1)π + π/2] 2 independent of the boundary conditions. The linear coefficientχ