Physical Aspects of Axonemal Beating and Swimming

We discuss a two-dimensional model for the dynamics of axonemal deformations driven by internally generated forces of molecular motors. Our model consists of an elastic filament pair connected by active elements. We derive the dynamic equations for this system in presence of internal forces. In the limit of small deformations, a perturbative approach allows us to calculate filament shapes and the tension profile. We demonstrate that periodic filament motion can be generated via a self-organization of elastic filaments and molecular motors. Oscillatory motion and the propagation of bending waves can occur for an initially non-moving state via an instability termed Hopf bifurcation. Close to this instability, the behavior of the system is shown to be independent of microscopic details of the axoneme and the force-generating mechanism. The oscillation frequency however does depend on properties of the molecular motors. We calculate the oscillation frequency at the bifurcation point and show that a large frequency range is accessible by varying the axonemal length between 1 and 50$\mu$m. We calculate the velocity of swimming of a flagellum and discuss the effects of boundary conditions and externally applied forces on the axonemal oscillations.


I. INTRODUCTION
Many small organisms and cells swim in a viscous environment using the active motion of cilia and flagella. These are hair-like appendages of the cell which can undergo periodic motion and use hydrodynamic friction to induce cellular self-propulsion (Bray 1992). In this paper, we are interested in those flagella and cilia which contain force generating elements integrated along the whole length of the elastic filamentous structure. They represent rod-like elastic structures which move and bend as a result of internal stresses. Examples for these systems are paramecium which has a large number of cilia on its surface; sperm, which use a single flagellum to swim; and chlamydomonas which uses two flagella to swim (Bray 1992). Cilia also occur in very different situations. An example is the kinocilium which exists in many hair bundles of mechanosensitive cells and has the ability to beat periodically (Rüsch and Thurm 1990).
The common structural theme of cilia and flagella is the axoneme, a characteristic structure which occurs in a large number of very different organisms and cells and which appeared early in evolution. The axoneme consists of a cylindrical arrangement of 9 doublets of parallel microtubules and one pair of microtubules in the center.
In addition, it contains a large number of other proteins such as nexin which provide elastic links between microtubule doublets, see Fig. 1. The axoneme is inherently active. A large number of dynein molecular motors are located in two rows between neighboring microtubules and can induce forces and local displacements between adjacent microtubules (Albert et al. 1994).
Axonemal flagella can generate periodic waving or beating patterns of motion. In the case of sperm for example, a bending wave of the flagellum propagates from the head which contains the chromosomes 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 inertia terms in the fluid hydrodynamics can be neglected (Taylor 1951). Therefore, only friction forces resulting from solvent viscosity can contribute to propulsion. Under such conditions self-propulsion is possible if a wave propagates towards one end, a situation which breaks time-reversal invariance of the sequence of deformations of the flagellum (Purcell 1977). How can bending waves be generated by axonemal dyneins which act internally within the axoneme? Dynein molecular motors induce relative forces between parallel elastic filaments, which are microtubule doublets. As a result of these forces, the filaments have the tendency to slide with respect to each other. If such a slid-ing is permitted globally, the filaments simply separate but no bending occurs. Bending results if 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 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 (Warner and Mitchell 1981). Small gold-beads specifically attached to microtubules in fully functioning flagella can be used to directly visualize the local relative sliding during beating (Brokaw 1991).
The theoretical problem of oscillatory axonemal bending and wave patterns have been addressed by several authors. One can distinguish two principally different mechanisms to generate oscillatory forces within the axoneme: (i) Deterministic forcing: a chemical oscillator could regulate the dynein motors which are activated and deactivated periodically. In this case, the regulatory system defines a dynamical force-pattern along the axoneme and drives the system in a deterministic way (Sugino and Naitoh 1982). (ii) Self-organized beating: the axoneme oscillates spontaneously as a result of the interplay of force-generating elements and the elastic filaments (Machin 1963, Brokaw 1975, Brokaw 1985, Lindemann and Kanous 1995. In particular, Brokaw has studied thoroughly self-organized patterns of beating using numerical simulations of simple models (Brokaw 1985, Brokaw 1999. Such models are based on the bending elasticity of the flagellum and an assumption on the coupling of motor activity to the flagellar deformations. In the present work, we are mostly interested in selforganized beating. Our approach is conceptually different from other works as we focus on an oscillating instability of the motor-filament system. In general, spontaneous oscillations occur via a so-called Hopf-bifurcation where an initially stable quiescent state becomes unstable and starts to oscillate. There is evidence for the existence of such a dynamic instability in the axoneme. Demembranated flagella show a behavior which depends on the ATP-concentration C AT P in the solution. For small C ATP , the flagellum is straight and not moving. If C ATP is increased, oscillatory motion sets in at a critical value of the ATP concentration and persists for larger concentrations (Gibbons 1975). This implies an instability of the initial straight state with respect to a wave-like 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 . This suggests that a system consisting only of molecular motors and semiflexible filaments can in general undergo self-organized oscillations. This idea is supported by the facts that flagellar dyneins are capable of generating oscillatory motion (Shingyoji et al. 1998) and that experiments suggest the existence of dynamic transitions in many-motor systems (Riveline et al. 1998, Fujita andIshiwata 1998). Patterns of 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 (Brokaw 1991).
In the following sections, we present a systematic study of a simplified model for the axoneme which has been introduced recently (Camalet et al. 1999) and which captures the basic physical properties that are relevant for its dynamics. In Section II, we present a thorough analysis of the dynamics of flexible filaments driven by internal forces. This work is inspired by recent studies of the dynamics of semiflexible filaments subject to external forces (Goldstein and Langer 1995. We show how the dynamic equations can be solved perturbatively and we calculate the shapes of bending waves, the velocity of swimming and the tension profile along the flagellum. In Section III, we briefly review a simple two-state model for a large number of coupled molecular motors. This model is well suited to represent the dynein molecular motors which act within the axoneme. Self-organized bending waves and oscillations via a Hopf bifurcation of the coupled motor-filament system are studied in section IV. We calculate the wave-patterns close to a Hopf bifurcation and determine the frequencies selected by the system. Finally, we discuss the relevance of our simple model to real axonemal cilia and flagella and propose experiments which could be performed to test predictions that follow from our work.

II. A SIMPLE MODEL FOR AXONEMAL DYNAMICS
The cylindrical arrangement of microtubule doublets within the axoneme can be modeled effectively as an elastic rod. Deformations of this rod lead to local sliding displacements of neighboring 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 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 "head", the two filaments are rigidly attached and not permitted to slide with respect to each other. Everywhere else, sliding is possible, see Fig. 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 r 1 (s) = r(s) − a n(s)/2 r 2 (s) = r(s) + a n(s)/2 , where n with n 2 = 1 is the filament normal. The local geometry along the filament pair is characterized by the relationsṙ where t denotes the normalized tangent vector and C = t · n is the local curvature. Throughout this paper dots denote derivatives with respect to s, i.e.ṙ ≡ ∂ s r ≡ ∂r/∂s. Internal forces f (s) are exerted in opposite directions, tangential to the filaments. The sliding displacement ∆ at the tail is indicated.

A. Bending and sliding of a filament pair
The two filaments are assumed to be incompressible and rigidly attached to each other at one end where s = 0. Bending of the filament pair and local sliding displacements are then 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 Eqns.

B. 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 large number of passive and active elements (e.g. nexin links and dynein molecular motors) which give rise to relative forces between neighboring microtubule doublets. These internal forces can be characterized by a coarse-grained description defining the force per unit length f (s) acting at position s in opposite directions on the two microtubules, see Fig. 2 (internal forces must balance). This force density is assumed to arise as the sum of active and passive forces generated by a large number of proteins. 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 constraintṙ 2 = 1 which ensures that s is the arclength (Goldstein and Langer 1995). The internal force density f couples to the sliding displacement ∆ as described by the contribution dsf ∆. The variation of this term under small deformations of the filament shape represents the work performed by internal stresses. Using Eq. (7) leads after a partial integration to where is the force density integrated to the tail. From Eq.(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 Here, plays the role of the physical tension. This becomes apparent since from Eq. (11) it follows that where F ext (L) is the external force applied at the end which satisfies Eq. (20). Therefore, τ is the tangent component of the integrated forces acting on the filament.

C. Dynamic equations
We derive the dynamic equation with the simplifying assumption 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 (Wiggins and ) 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 Eq. (11) Noting that ∂ tṙ = n∂ t ψ, we obtain an equation of motion for ψ(s) alone: The tension τ is determined by the constraint of incompressibility ∂ tṙ 2 = 2t · ∂ tṙ = 0. This condition and Eq. (15) leads to a differential equation for the tension profile: Eqns. (16) and (17) determine the filament dynamics. The filament shape follows from where r(0, t) can be obtained from Eq. (15) evaluated at s = 0.

D. Boundary Conditions
The filament dynamics 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, 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 I: Case A is the situation of a filament with 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.

E. Small deformations
In the absence of internal forces f , the filament relaxes passively to a straight rod withψ = 0, i.e. ψ constant. Without loss of generality, we choose ψ = 0 in this state which implies that straight filament is parallel to the Xaxis. Internal stresses f (s) induce deformations of this straight conformation. For small internal stresses, we can perform a systematic expansion of the filament dynamics in powers of the stress amplitude.
We introduce a dimensionless parameter ǫ which scales the amplitudes of the internal stresses, f (s, t) = ǫf 1 (s, t), where f 1 is an arbitrary stress distribution. We can now solve the dynamic equations (16) and (17) perturbatively by writing which allows us to determine the coefficients ψ n (s, t) and τ n (s, t). This procedure is described in Appendix B. We find that ψ 2 and τ 1 always vanish and τ 0 = σ is a constant tension which is equal to the X-component of the external force applied at the end. Note, that for boundary conditions A-C, σ = 0. The filament shape is thus characterized by the behavior of ψ 1 (s, t) which obeys In order to discuss the filament motion in space, we define the average velocity of swimminḡ 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 motion is towards the positive or negative X-direction, see Fig.  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 Longitudinal displacements satisfy The boundary conditions for h(s) are given in table II 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 ψ.

F. Swimming
The velocity of swimmingv can be calculated perturbatively in ǫ. We consider for simplicity periodic internal stresses with frequency ω given by 1 which after long time leads to periodic filament motion whereh satisfies according to Eq. (26) Writingh = He −iφ where H and φ denote the amplitude and phase, respectively, filament motion can be expressed as which represents bending waves for which v p = ω/φ can be interpreted as the local wave propagation velocity. Examples of bending waves for self-organized beating discussed in section IV are displayed in Figs. 4 and 5. For a freely oscillating filament with a viscous load of friction coefficient ζ attached (case C), we find an average is the velocity for ζ = 0. Eq. (32) which is correct to second order in ǫ reveals that motion is only possible if the filament friction is anisotropic ξ /ξ ⊥ = 1 and if a wave is propagating,φ = 0. This is consistent with earlier work on swimming (Taylor 1951, Purcell 1977, Stone and Samuel 1996. For a rod-like filament with ξ < ξ ⊥ , swimming motion is opposite to the direction of wave propagation.

III. INTERNAL FORCES GENERATED BY MOLECULAR MOTORS
The filament dynamics is driven by internal stresses f (s, t) generated by a large number of dynein motors and passive elements. We describe the properties of a forcegenerating system using a simplified two-state model.

A. Two state model for molecular motors
We briefly review the two-state model for a large number of motor molecules attached to a filament which slides with respect to a second filament Prost 1995, Jülicher et al. 1997). Each motor is assumed to have two different chemical states, a strongly bound state 1 and a weakly bound state or detached state 2. The interaction between a motor and a filament in states 1 and 2 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 reflect the symmetry of the filaments: 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 which is the chemical fuel that drives the system, the motors are assumed to undergo transitions between states. The corresponding transition rates are denoted ω 1 for detachments and ω 2 for attachments. We introduce the relative position ξ of a motor with respect to the potential period where x = ξ + nl, 0 ≤ ξ < l and n is an integer. The probability to find a motor in state i at position ξ at time t is denoted P i (ξ, t), P 1 + P 2 is normalized within one period. The dynamic equations of the system are given by ∂ 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 x is determined by the forcebalance Here, the coefficient λ describes the total friction per unit length in the system. The number density of motors along the filament is denoted by ρ and f is the force per unit length generated by the system. The elastic modulus per unit length K occurs in presence of elastic elements such as nexins. If motors have an incommensurate arrangement compared to the filament periodicity, P 1 + P 2 = 1/l and the equations simplify and can be expressed by P 1 alone where ∆W = W 1 − W 2 .

B. Molecular motors coupled to a filament pair
The two state model for a large number of motors is well suited to represent the internal forces acting within the axoneme. We assume that at any position s an independent two-state model described by Eq. (35) is located which generates the internal force density f (s). The filament sliding displacement is identified with the motor displacement, ∆ ≡ x. Note, that here we neglect for simplicity fluctuations which arise from the chemical activity of a finite number of force generators and we assume homogeneity of all properties along the axonemal length. The energy source of the active system is the chemical activity characterized by the transition rates ω 1 and ω 2 , the generated forces induce bending deformations of the filaments.
The dynamic equations (35) represent a nonlinear active system which generates time-dependent forces. Since we will study periodic motion, we can express forces and displacements by a Fourier expansion The relation between force amplitudes f n and amplitudes of sliding displacements ∆ n can in general be expressed in powers of the amplitudes ∆ n (Jülicher and Prost 1997) Here, the summation over common indices is implied. The coefficients F (n) are calculated for the two-state model in Appendix C. However, Eq. (38) is more general and characterizes a whole class of active nonlinear mechanical systems.

C. Symmetry considerations
The interaction of the motors with the two filaments is asymmetric. Motors are rigidly connected to one of the filaments and slide along the second. In addition, the filaments are polar and filament sliding is induced by motors towards one particular end of the two filaments. As a consequence, the 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 of bending deformations. This broken symmetry does however not reflect the symmetry of the axoneme which is symmetric with respect to all microtubule doublets. Each of the nine microtubule doublet 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 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 our two-dimensional model, we assume that both filaments have motors which are rigidly attached and which slide on the second filament, see Fig. 6. The consequence of this symmetry is that the filament polarity becomes unimportant. Exchanging the role of the two filaments is equivalent to replacing the sliding displacement ∆ of motors by −∆. Therefore, in the context of the two-state model this symmetrization is equivalent to the requirement that energy landscapes W i (x) and transition rates ω i (x) are symmetric functions with respect to x → −x. As a consequence of this symmetry, there is no preferred direction of bending. In the following, we adopt this choice of a symmetric motor/filament pair. If we use the expansion given in Eq. (38), this symmetry implies that all even coefficients F (2n) = 0.

D. Linear response function
The perturbative treatment to study small filament deformations introduced in section II can be naturally extended to the coupled motor-filament situation using the expansion (38). Up to second order in ǫ, only the linear response coefficient of the active system plays a role which for the two-state model is given by F (1) nk = χδ nk with Here, we have introduced R ≡ ω 2 /αl and α = ω 1 + ω 2 . Higher order terms F (2n+1) have to be taken into account if the third or higher order in ǫ is considered. The linear response function χ as well as nonlinear coefficients can be calculated most easily for a simple choice of symmetric potential and transition rates ∆W (ξ) = U cos(2πξ/l) ω 1 (ξ) = β − β cos(2πξ/l) where α and β are ξ-independent rate constants. For this convenient choice where k ≡ U/l 2 is the cross-bridge elasticity of the motors. We have introduced the dimensionless parameter Ω = 2π 2 β/α with 0 < Ω < π 2 which plays the role of a control parameter of the motor-filament system, α is a characteristic ATP cycling rate.

IV. 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, selforganized oscillations of the driven filament pair are a natural consequence of its physical properties. The control parameter Ω is in our two-state model a ratio of chemical rates of the ATP hydrolysis cycle. In an experimental situation, it could be varied e.g. by changing the ATP concentration. If the molecular motors are regulated by some other ion concentration such as e.g. Ca 2+ , this concentration could also play the role of the control parameter.

A. Generic aspects
For oscillations in the vicinity of a Hopf-bifurcation, filament deformations are small and we can thus use Eq. (26) to describe the filament dynamics. 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 Eq. (30), we find that spontaneously oscillating modesh(s) at frequency ω are to linear order near the bifurcation solutions to 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 II by replacing h →h, ∂ t h → iωh and f →f = aχ(∂ sh − ∂ sh (0)). As discussed in Appendix D, Eq. (44) has the structure of an eigenvalue problem, with χ playing the role of a complex eigenvalue. For every choice of parameters, there exists an infinite set of nontrivial solutionsh n to Eq. (44) if χ is equal to the corresponding eigenvalue χ n , n = 1, 2, ... We order these values according to their amplitude: |χ n | ≤ |χ n+1 |. The functional dependence of the Eigenvalues χ n on the model parameters is completely determined by dimensionless functionsχ n (σ,ω) according to whereσ ≡ σL 2 /κ andω ≡ ωξ ⊥ L 4 /κ are a dimensionless tension and frequency, respectively. 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 details of the problem such as the structural complexity within the axoneme. The motion is given by one of the modesh n which only depend onσ andω.   Fig. 5 for boundary conditions A-C as snapshots taken at different times. The oscillation amplitudes are chosen in such a way that the maximal value of ψ ≃ π/2. The direction of wave propagation depends on the boundary conditions. For clamped head (case A), waves propagate towards the head whereas for fixed or free head (cases B and C), they propagate towards the tail. The positions of the lowest eigenvalues χ n in the complex plane for boundary conditions B are displayed in Fig. 7 forσ = 0 and varyingω.

B. 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, explicit knowledge of the linear response function χ(Ω, ω) is necessary. As a simple example, we use χ as given by Eq. (41) for a two-state model. 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, whereω c = ξ ⊥ L 4 ω c /κ. Since χ is complex, both the real and the imaginary part of Eq. (46) represent independent conditions. Therefore, Eq. (46) 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 |χ(Ω, ω)| typically 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 .

C. Axonemal vibrations for different lengths
We choose the parameters of our model to correspond to the axonemal structure. We estimate the bending rigidity κ ≃ 4 · 10 −22 Nm 2 which is the rigidity of about 20 microtubules. Furthermore, we choose a microtubule separation a ≃ 20nm, motor density ρ ≃ 5 · 10 8 m −1 , friction per unit length λ ≃ 1kg/ms, rate constant α ≃ 10 3 s −1 , cross-bridge elasticity k ≃ 10 −3 N/m and a perpendicular friction ξ ⊥ ≃ 10 −3 kg/ms 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 Fig. 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 Ldependent and the system propagates bending waves of a wave-length shorter than the filament length.
The coefficient γ depends on boundary conditions but not on any model parameters. For a clamped head (A), γ ≃ 0.184, for a fixed head (B), γ ≃ 0.008, see Appendix D. The criterion for a Hopf bifurcation for small L is Together with Eq.(41), 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µm for our choice of the parameters. The maximal frequency of oscillations occurs for L = L min . Within the range of lengths between 1µm and 50µm, we find frequencies ω c /2π which vary between tens of Hz and about 20kHz. The parameters chosen in Fig.  8 lead to frequencies above 150Hz. Lower frequencies are obtained for larger values of the friction λ or smaller chemical rate α.

D. 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 eigen-valuesχ n depend onσ, therefore the tension affects shape and frequency of an unstable mode. Inspection of Eq. (44) 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 nontrivial. However, for a filament with 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 Eq. (44) and in the boundary conditions and 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 Eq. (41) by σ/a 2 . Therefore, we can include the effects of tension on the critical frequency in Eq. (48) 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 Eq. (49), we find For our choice of parameters, ρka 2 ≃ 200pN which indicates that forces of the order of 10 2 pN should have a significant effect on the bifurcation. Furthermore, Eq.
(51) 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.

V. DISCUSSION
In the previous sections, we have shown that many aspects of axonemal eating can be described by a simple model based on the idea of local sliding of microtubule doublets driven by molecular motors inside the axoneme. This model represents a class of physical systems termed internally driven filaments (Camalet et al. 1999) which have characteristic properties that are closely linked to the geometric constraint that couples global bending and local filament sliding.
Our work shows that axonemal beating can be studied both numerically and analytically using a coarse-grained description which ignores many details of the proteins involved and is based on effective material properties such as bending rigidities, internal friction coefficient and elastic modulus per unit length as well as the frequency dependent linear response function of the active elements. Our main results are: (i) The filament pair introduces a geometric coupling of the active elements at different position along the filaments. This coupling is suited for the generation of periodic beating motion by a general selforganization mechanism of the system. (ii) The generation of oscillations and propagating bending waves occurs via a Hopf bifurcation for a critical value of a controlparameter such as the ATP concentration. The patterns of motion generated close to this bifurcation can be calculated without any knowledge of the internal force generating mechanism if the oscillation frequency is known. (iii) The frequency depends on chemical rates of the motors, coefficients of internal elasticities and frictions, the solvent viscosity and the microtubule rigidity. Our calculations using a simple two-state model for molecular motors suggest significant variations of the oscillation frequency if the axonemal length is varied. Short axonemes are predicted to be able to oscillate at frequencies of several kHz. This fact is particularly interesting in the case of kinocilia which are located in the hair bundles of many mechanosensitive cells that are involved in the detection of sounds. If cilia are able to vibrate at high frequencies, this suggests that the kinocilium could play as active role in sound detection (Camalet et al. 2000). (iv) Forces applied to the axoneme can be an important tool to study the mechanism of force generation. An external force can play the role of a second control parameter for the bifurcation which has a stabilizing effect for increasing forces.
These results have been obtained using several simplifying assumptions. Our model represents the solvent hydrodynamics by an anisotropic local friction acting on the rod-like filament. This approach ignores hydrodynamic interactions between different parts of the filament which lead to logarithmic corrections. These hydrodynamic effects do not change the basic physics but can lead to corrections to the numerical results. It is straightforward to generalize our model to incorporate the effects of hydrodynamic interaction, but this does not change our results qualitatively. A more serious simplification of our model is the restriction to a two-dimensional system and to filament configurations which are planar. This choice is motivated by the fact that observed bending waves of flagella are planar in many cases (Brokaw 1991). If filament configurations are planar then our model applies and is complete. However, this model cannot explain why motion is confined to a plane and it misses all non-planar modes of beating. In order to address such questions, a three dimensional generalizations have to be used. Such generalizations introduce additional aspects. In particular, torsional deformations become relevant. The local sliding displacements depend in the three dimensional case both on bending and torsional deformations and the full tor-sional dynamics has to be accounted for in the dynamic equations. Finally, we have restricted most calculations to the limit of small deformations which corresponds to filament shapes that are almost straight. This regime has several important features. The filament dynamics can be studied by an analytic approach by a systematic expansion in the deformation amplitude. This allows us to characterize linear and nonlinear terms both for the filament dynamics as well as for the properties of the active elements. Patterns of motion close to a Hopf bifurcation are fully characterized by the linear terms of this expansion. These linear terms are given by the structure of the problem and their form does not depend on molecular details of the system. For example, most details of the operation of the molecular motors are unimportant for the shape of the filament oscillations at the bifurcation, only the linear response function of the active material plays a role.
If filament beating with larger amplitude is of interest, nonlinearities become relevant. Nonlinearities arise due to nonlinear geometric terms in the bending energy or via non-Hookian corrections of the elasticity. Furthermore, nonlinearities in the force-generation process of molecular motors exist. All these nonlinearities determine the large amplitude motion and could give rise to new types of behavior such as additional dynamic instabilities. However, in contrast to the linear terms, the form of the dominant nonlinearities does depend on structural details of the axoneme. Therefore, an analysis of large amplitude motion is difficult and knowledge of the nature of the dominating nonlinearities is required. However, if no new instabilities after the initial Hopf-bifurcation occur, the principal effect of nonlinear terms is to fix the amplitude of propagating waves. We therefore expect that propagating waves with larger amplitude are in many cases well approximated by our calculation to linear order.
Our work shows that propagating bending waves and oscillatory motion of internally driven filaments can occur naturally by a simple physical mechanism. Complex biochemical networks to control the system are not required for wave propagation to occur. This suggests that the basic axonemal structure intrinsically has the ability to oscillate. Biochemical regulation systems are thus expected to control this activity on a higher level but are not responsible for oscillations. Our work shows that experimental studies of axonemal beating close to a Hopf bifurcation would be very valuable. Such experiments could be e.g. performed by using demembranated flagella and ATP concentrations not far from the level where beating sets in. The observed behavior close to the bifurcation would give insight in the self-organization at work. Furthermore, externally applied forces and manipulation of boundary conditions could be sensitive tools to test some of the predictions of our work. The variation δG of the enthalpy G given by Eq. (8) under variations δr of a shape r(s) is Note, thatṙ 2 = 1 but (ṙ+δṙ) 2 = 1 in general. Under such a variation, δC = δ(n·r) = n·δr sincer·δn = Cn·δn = 0. Therefore, Two subsequent partial integrations lead to The functional derivative δG/δr given by Eq. (11) can be read off from the integrand, the boundary terms provide expressions for forces and torques at the ends given by Eqs. (19) and (20).