key: cord-0303902-9fy5ska6
authors: Cicconofri, Giancarlo; Noselli, Giovanni; DeSimone, Antonio
title: The biomechanical role of extra-axonemal structures in shaping the flagellar beat of Euglena
date: 2020-03-15
journal: bioRxiv
DOI: 10.1101/2020.03.15.991331
sha: 1c94d7a308acd1d96bd486746e13c6d7d274762f
doc_id: 303902
cord_uid: 9fy5ska6

We propose and discuss a model for flagellar mechanics in Euglena gracilis. We show that the peculiar non-planar shapes of its beating flagellum, dubbed “spinning lasso”, arise from the mechanical interactions between two of its inner components, namely, the axoneme and the paraflagellar rod. The spontaneous shape of the axoneme and the resting shape of the paraflagellar rod are incompatible. The complex non-planar configurations of the coupled system emerge as the energetically optimal compromise between the two antagonistic components. The model is able to reproduce the experimentally observed flagellar beats and their characteristic spinning lasso geometric signature, namely, travelling waves of torsion with alternating sing along the length of the flagellum.

Flagella and cilia propel swimming eukaryotic cells and drive fluids on epithelial tissues of higher organisms [1] . The inner structure of the eukaryotic flagellum is an arrangement of microtubules (MTs) and accessory proteins called the axoneme (Ax). A highly conserved structure in evolution, the Ax typically consists of nine cylindrically arranged MT doublets cross-bridged by motor proteins of the dynein family. An internal central pair of MTs is connected by radial spokes to the nine peripheral doublets, determining the typical "9+2" axonemal structure. Motor proteins hydrolyze ATP to generate forces that induce doublet sliding. Due to mechanical constraints exerted by linking proteins (nexins) and the basal body, dynein-induced sliding of MTs translates into bending movements of the whole structure. Motor proteins are thought to self regulate their activity via mechanical feedback, generating the periodic beats of flagella, see e.g. [6] and [19] .

Despite a general consensus on the existence of a self-regulatory mechanism, the inner working of the Ax is still not fully understood and it is still a subject of active research [33] . While bending-through-sliding is the accepted fundamental mechanism of flagellar motility, how specific flagellar shapes are determined is not yet clear. For the most studied swimming microorganisms, such as animal sperm cells and the biflagellate alga Chlamydomonas reinhardtii, the flagellar beat is, to a good approximation, planar. For these organisms, beat planarity is thought to be induced by the inter-doublet links between one pair of MTs, typically those numbered 5 and 6 [18] . These links inhibit the relative sliding of the 5-6 MTs pair, thus selecting a beating plane that passes through the center of the Ax and the midpoint between the inhibited MTs. Figure 1 . a) A specimen of freely swimming Euglena gracilis, and b) a sketch of the cross section of its flagellum seen from the distal end. The flagellar inner structure is composed by the paraflagellar rod (PFR, textured), and the axoneme (Ax). The PFR is connected via bonding links to the axonemal doublets 1, 2, and 9 (our numbering convention). The inner structure of the flagellum is enclosed by the flagellar membrane (dotted line). By inhibiting MTs' sliding the PFR selects the spontaneous bending plane of the Ax (dashed line). The solid line that joins the Ax center r a and the PFR center r p cross at an angle φ p the spontaneous bending plane.

A remarkable deviation from the flagellar structure of the aforementioned organisms is found in euglenids and kinetoplastids. These flagellated protists have a whole extra element attached alongside the Ax [7] , a slender structure made of a lattice-like arrangement of proteins called "paraxial" or "paraflagellar" rod (PFR), see Figure 1 . The latter name is more common, but the former is possibly more accurate [24] . PFRs are attached via bonding links to up to four axonemal MTs, depending on the species [32] . PFRs are thought to be passive but, at least in the case of Euglena, some degree of activity is not completely ruled out [21] .

In this paper, we put forward and test the hypothesis that the distinctive beating style of Euglena Gracilis, sometimes dubbed "spinning lasso" [5] , arises due to PFR-Ax mechanical interaction.

In order to put our hypothesis into context, we observe that the flagellar beat of PFR-bearing kinetoplastid organisms, such as Leishmania and Crithidia, is planar [9] . An apparent exception to beat planarity in kinetoplastids is found in the pathogenic parasite Trypanosoma brucei, which shows a characteristic non-planar "drill-like" motion [17] . However, Trypanosoma's flagellum is not free, like that of Leishmania and Crithidia, but it is attached to the organism for most of its length, wrapped helically around the cell body. According to [2] the flagellum-body mechanical interaction could alone explain Trypanosomes' distinctive motion. Confirming this conclusion, [34] showed that Trypanosoma mutants with body-detached flagellum generate fairly planar beating. It is conjectured that the PFR-Ax bonds operate as the 5-6 interdoublet links in Chlamydomonas and sperm cells, inhibiting MTs sliding and selecting a plane of 2/24 beat [36] .

Euglena's spinning lasso beat does not conform to this scenario. Indeed, the beating style of Euglena is characterized by high asymmetry and non-planarity. The full 3d flagellar kinematics of freely swimming cells has recently been revealed [25] thanks to a mixed approach based on hydrodynamic theory and image analysis. As we report in the first part of this paper, the geometry of the spinning lasso is characterized by travelling waves of torsion with alternating sign along the flagellum length.

We argue that the key to the emergence of non-planarity lies in a prominent asymmetry in the structure of PFR-Ax attachment in euglenid flagella. Figure 1 shows a sketch of the cross section of the euglenid flagellum redrawn from [20] , see also [4] . The PFR is attached to three MTs, which we number 1, 2, and 9. We consider two lines. One line (dashed) passes through the center of the Ax and MT 1, in the middle of the bonding complex. The other line (solid) connects the center of the Ax and the center of the PFR. The two lines cross each other. This is the structural feature on which we build our model.

In the model we assume that the bonding links to the PFR select the local spontaneous beating plane of the Ax, from the same principle of MTs' sliding inhibition discussed above. The local spontaneous beating plane so generated passes through the dashed line in Figure 1 . We follow closely [14] and [27] in our modeling of the Ax, while we use a simple elastic spring model for the PFR. We show that, under generic actuation, the two flagellar components cannot be simultaneously in their respective states of minimal energy, and this crucially depends on the offset between the spontaneous beating plane of the Ax (dashed line in Figure 1 ) and the line joining the PFR-Ax centers (solid line in Figure 1 ). Instead, the typical outcome is an elastically frustrated configuration of the system, in which the two competing components drive each other out of plane. Under dyneins activation patterns that, in absence of extra-axonemal structures, would produce an asymmetric beat similar to those of Chlamydomonas [23] , or Volvox [26] , the model specifically predicts the torsional signature of the spinning lasso, which we discuss in the following Section.

Interestingly, the lack of symmetry of the spinning lasso beat produces swimming trajectories with rotations coupled with translations [25] . In turn, cell body rotations have a role in phototaxis (see Discussion and Outlook Section and [12] ). In the light of these observations, our analysis shows that the beat of the euglenid flagellum can be seen as an example of a biological function arising from the competition between antagonistic structural components. It is not dissimilar from the body-flagellum interaction in Trypanosoma, which generates 3d motility. But the principle is much more general in biology and many other examples can be found across kingdoms and species and at widely different scales. For instance in plants, a mechanism of seed dispersal arises from the mechanical competition between the two valves of the seed pods, see, e.g., [3] for Bauhinia variegata and [15] for Cardamine hirsuta. Contraction by antagonistic muscles is key for animal movement and, in particular, for the functioning of hydrostatic skeletons (used from wormlike invertebrates to arms and tentacles of cephalopods, to the trunk of elephants, see [16] ). A similar principle of antagonistic contraction along perpendicularly oriented families of fibers is at work at the sub-cellualr level, for example in the antagonistic motor protein dynamics in contractile ring structures important in cell division and development (see, e.g., [8] ). At the same sub-cellular scale, competing elastic forces arising from lipid-protein interactions are often crucial in determining the stability of complex shapes of the cellular membrane [31] , and in the case of the overall structure of the coronavirus envelope [29] . 

We first analyze the experimental data from the 3d reconstruction of the beating euglenid flagellum obtained in [25] for freely swimming organisms. Swimming Euglenas follow generalized helical trajectories coupled with rotation around the major axis of the cell body. It is precisely this rotation that allows for a 3d reconstruction of flagellar shapes from 2d videomicroscopy images. Euglenas take many beats to close one complete turn around their major body axis. So, while rotating, Euglenas show their flagellar beat to the observer from many different sides. Stereomatching techniques can then be employed to reconstruct the flagellar beat in full (assuming periodicity and regularity of the beat). Figure 2 shows N = 10 different curves in space describing the euglenid flagellum in different instants within a beat taken from [25] . The reconstruction fits well experimental data from multiple specimens. The figure also illustrates the calculated torsion of the flagellar curve at each instant (not previously published). Torsion, the rate of change of the binormal vector, is the geometric quantity that measures the deviation of a curve from a planar path (see the Results Section below for the formal definition). The spinning lasso shows here a distinct torsional signature characterized by torsion peaks of alternate sign, traveling from the proximal to the distal end of the flagellum.

To further investigate Euglena's flagellar beat we observed stationary cells trapped at the tip of a capillary. In this setting the flagellum is not perturbed by the hydrodynamic forces associated with Euglena's rototranslating swimming motion. The beat can then manifest itself in its most "pristine" form. We recorded trapped Euglenas during periodic beating. Then, we rotated the capillary and recorded the same beating cell from different viewpoints. Videomicroscopy images from one specimen are shown in . While with fixed specimens we cannot reconstruct reliably the 3d flagellar shapes, Figure 3 shows that there is a high stereographical consistency with the flagellar shapes obtained from swimming organisms. Flagellar non-planarity is thus not intrinsically associated with swimming, which reinforce the idea that the mechanism that generates non-planar flagellar shapes might be structural in origin. Moreover, these observations justify the choice we made in our study to focus on a model of flagellar mechanics for stationary organisms, allowing for substantial simplifications. As a final remark we observe how, from a simple geometric construction, we can show that the torsional pattern in Figure 2 is consistent with Euglena's flagellar shapes as seen from common 2d microscopy, for either swimming or trapped organisms. Typically, the 2d outline (i.e. the projection on the optic plane) of a beating euglenid flagellum is that of a looping curve, see e.g. [30] for independent observations. Consider now an idealized 3d model of the spinning lasso geometry: a curve with two singular points of concentrated torsion with opposite sign, such as the one shown in Figure 3 . If we move along the curve, from proximal end to distal end, we first remain on a fixed plane (blue). Then the plane of the curve abruptly rotates by 90 degrees (red plane) first, and then back by 90 degrees in the opposite direction (yellow plane). These abrupt changes correspond to concentrated torsional peaks of opposite sign. When seen in a two dimensional projection, this torsion dipole generates a looping curve that closely matches euglenid flagella's outlines during a spinning lasso beat.

We model Ax and PFR as cylindrical structures with deformable centerlines, see Figure 4 . The euglenid flagellum is the composite structure consisting of Ax and PFR attached together. We suppose that the Ax is the only active component of the flagellum, whereas the PFR is purely passive. Our mechanical model builds on the definition of the total internal energy of the flagellum

which is given by the sum of three terms: the passive (elastic) internal energy W a pas of the Ax, the active internal energy W a act of the Ax (generated by dynein action), and the (passive, elastic) internal energy W p of the PFR. The passive internal energy of the Ax is given by

which is formally identical to the classical expression for the energy of elastic (inextensible) rods. We denoted with U 1 and U 2 the bending strains, U 3 is the twist, while B a and C a are the bending and twist moduli, respectively. L is the total length of the Ax centerline r a . Bending strains and twist depend on the arc length s of the centerline, and they are defined as follows. We associate to the curve r a an orthonormal frame d i (s), with i = 1, 2, 3, which determines the orientation of the orthogonal sections of the Ax (enclosed by light blue circles in Figure 4 ). The unit vectors d 1 (s) and d 2 (s) define the plane of the orthogonal section at s. The unit vector d 3 (s) = ∂ s r a (s) lies perpendicular to the section. Bending strains and twist are then given by

Thus, U 1 and U 2 measure the bending of the Ax on the local planes d 1 -d 3 and d 2 -d 3 , respectively, while the twist U 3 is given by the rotation rate of the orthonormal frame around the centerline's tangent d 3 .

We remark here that we do not consider the Ax as a "filled" beam but as a hollow tubular structure; its elasticity comes from the individual MTs lying on its outer surface. The derivation of (2) from a detailed model of the Ax is given in Appendix A.

The active internal energy of the Ax is defined as minus the total mechanical work of the dyneins

where γ 1 and γ 2 are the two variables that quantifies the shear (i.e. collective sliding) of MTs, while H 1 and H 2 are the corresponding shear forces exerted by molecular motors. Following [28] we also allow for singular shear forces, H 1 and H 2 , concentrated at the distal end of the Ax. This forces arise naturally, as we remark after (22) in the Results Section. The shear variables are related to the MTs' kinematics by the following formulas. The MTs' centerlines r j , for j = 1, 2 . . . 9 , are given by r j (s) = C(s, φ j ), where φ j = 2π(j − 1)/9, and

is the parametrization of the cylindrical surface of the Ax (ρ a is the Ax radius) in terms of the centerline arc length s and the angle φ. In Appendix A we show how the shear variables are related to the individual sliding of MTs with respect to their neighbours 6/24 (35) , and we compute the relations between the dynein forces acting on each pair of adjacent MTs and the shear forces H 1 and H 2 (36) . A special axonemal deformation with γ 2 = 0 is shown in Figure 4 . The Ax is in this case bent into a circular arc, and the centerline r a lies on the plane generated by the unit vectors d 1 and d 3 . The shear variable γ 1 (s) = 0 has here a simple geometrical interpretation. For each fixed s the curve φ → C(s, φ) describes what we call the "material" section of the Ax at s (red curves in Figure 4 ). The material section is a planar ellipse centered in r a (s) which connects points of neighbouring MTs' corresponding to the same arc length. Formula (5) says that γ 1 (s) is the tangent of the angle at which the material sections at s intersect the orthogonal sections at s. Shear variables and bending strains are coupled

The above formulas, which result from the constrained kinematics of the axonemal structure (as explained further in Appendix A), underlie the essential mechanism of axonemal motility: collective sliding of MTs generates bending of the whole Ax. We point out here that there is no coupling between the shear variables γ 1 , γ 2 and the twist U 3 , a fact that will have consequences in the remainder. The special axonemal deformation in Figure 4 shows the case in which U 1 (s) = 0 and U 2 (s) = K, so the Ax is bent into a circular arc of radius 1/K. While γ 2 (s) = 0, the shear variable γ 1 (s) = Ks increases linearly with s. Material and orthogonal sections coincide at the base (the basal body impose no shear at s = 0) and the angle between them grows as we move along the centerline towards the distal end of the Ax. In order for the Ax to bend, MTs from one side of the Ax must be driven driven towards the distal end while the others must be driven toward the proximal end.

We remark here that (4) defines the most general active internal energy generated by molecular motors, and we do not assume at this stage any specific (spatial) organization of dynein forces. We will introduce specific shear forces later in the Results Section.

The PFR is modeled as an elastic cylinder with circular cross sections of radius ρ p and rest length L. We assume that the PFR can stretch and shear. The total internal energy of the PFR is given by

where V 1 and V 2 are the shear strains, V 3 is the stretch, D p and E p are the shear and stretching moduli, respectively. We are neglecting here the PFR's bending and twisting resistance. Classical estimations on homogeneous elastic rods, see e.g. [13] , show that bending and twist moduli scale with the forth power of the cross section radius, whereas shear and stretching moduli scale with the second power and hence they are dominant for small radii. We assume that dynein forces are strong enough to induce shear in the PFR, thus PFR's bending and twist contributions to the energy of the flagellum become negligible. We are also neglecting Poisson effects by treating the PFR cross sections as rigid.

The PFR shear strains and stretch are defined as follows. The cross sections centers of the PFR lie on the curve r d , and their orientations are given by the orthonormal frame g i (s), with i = 1, 2, 3. The unit vectors g 1 (s) and g 2 (s) determine the cross section plane centered at r p (s), while the unit vector g 3 (s) is orthogonal to it. The curve r p is not parametrized by arc length and g 3 is not in general aligned with the tangent to r p . Shear strains and stretch are given by the formulas Figure 4 . a) Geometry of the Ax. MTs lie on a tubular surface C(s, φ) parametrized by generalized polar coordinates s and φ, where s is the arc length of the axonemal centerline r a . The unit vectors d 1 (s) and d 2 (s) lie on the orthogonal cross sections of the Ax (light blue circles). The material sections of the Ax are given by the curves φ → C(s, φ) (red), which connect points of neighbouring axonemal MTs corresponding to the same arc length s. Bend deformations of the axoneme are generated by the shear (collective sliding) of MTs. The shear is quantified by the angle between the orthogonal sections and the material sections of the Ax. b) Geometry of the euglenid flagellum, detail of the Ax-PFR attachment. The unit vectors g 1 (s) and g 2 (s) generate the plane of the PFR's cross sections. The vector g 1 (s) is parallel to the outer unit normal to the axonemal surface N(s, φ p ), while g 2 (s) is parallel to the tangent vector to the material section ∂ φ C(s, φ p ).

The shear strains thus depend on the orientation of the cross sections with respect to the centerline (tangent), while the stretch measures the elongation of the centerline.

The PFR-Ax attachment couples the kinematics of the two substructures, see Figure 4 . In the remainder we formalize the attachment constraint and we show how the PFR's shear strains and stretch (8) , and thus the flagellar energy (1), are completely determined by the Ax kinematic variables.

For each s, the PFR cross section centered at r p (s) is in contact with the Ax surface at the point C(s, φ p ) for for a fixed angle coordinate φ p , see Figure 1 . The PFR centerline is given by

where N(s, φ p ) ≈ d 1 (s) cos φ p + d 2 (s) sin φ p is the outer unit normal to the axonemal surface at C(s, φ p ). The normal vector N(s, φ p ) lies on the plane of the PFR cross section centered at r p (s). Indeed, we have g 1 (s) = N(s, φ p ) for the first unit vector of the PFR orthonormal frame. Only one more degree of freedom remains, namely g 2 (s), which must be orthogonal to N(s, φ p ), to fully characterize the orientations of the PFR cross sections. Here is where the bonding links attachments are introduced in the model. The bonding links of the PFR cross section centered at r p (s) are attached to three adjacent MTs at the same MTs' arc length s. The individual attachments are therefore located on the material section of the Ax at s. Given this, g 2 (s) is imposed to be parallel to ∂ φ C(s, φ p ), the tangent vector to the material section of the Ax at the point of contact C(s, φ p ), see Figure 4 . This condition critically couples MTs' shear to the orientations of the PFR cross sections, as further demonstrated below.

To summarize, we have the following formulas for the PFR orthonormal frame vectors

and g 3 (s) = g 1 (s) × g 2 (s) .

By replacing the expression (9)- (10) in (8), we obtain formulas for the shear strains and 8/24

stretch of the PFR in terms of the Ax kinematic parameters. We have that the shear strain V 1 and the stretch V 3 are of order ρ p ∼ ρ a (see Appendix A for detailed calculations). Since ρ p is small compared to the length scale L of both PFR and Ax we can neglect these quantities. The only non-negligible contribution to the PFR energy is thus given by the shear strain V 2 . After linearization, we have V 2 ≈ − sin φ p γ 1 + cos φ p γ 2 . The PFR energy in terms of Ax kinematic parameters is then given by

The shear of axonemal MTs determines the orientation of the PFR cross sections. In Figure 5 (middle pictures) we show an example of this kinematic interplay. The Ax is again bent in an arc of a circle on the plane d 1 − d 3 , with U 1 (s) = 0, U 2 (s) = K, γ 1 (s) = Ks, and γ 2 (s) = 0. PFR and Ax centerlines run parallel to each other, indeed from (9) we have that ∂ s r a ≈ ∂ s r p for every deformation. The linking bonds impose a rotation of the cross sections of the PFR as we progress from the proximal to the distal end of the flagellum, generating shear strain V 2 (s) = − sin φ p γ 1 (s) = − sin φ p Ks on the PFR. This mechanical interplay leads to non-planarity of the euglenid flagellar beat. This mechanism is controlled by the offset between the PFR-Ax joining line and the local spontaneous bending plane of the Ax, as further discussed in the Results Section.

Under generic (steady) dynein actuation, i.e. given H 1 and H 2 (not time-dependent), and in the absence of external forces, the flagellum deforms to its equilibrium configuration δW = 0. Bending strains and twist at equilibrium solve the equations

where U = (U 1 , U 2 ) is the bending vector, e p = (cos φ p , sin φ p ), and H ⊥ = (−H 2 , H 1 ). We use the symbol a ⊗ b to denote the matrix with components (a ⊗ b) ij = a i b j . The field equations (12) 

where H ⊥ = (− H 2 , H 1 ). Equations (12)-(13) can be interpreted as the torque balance equations of the Ax. The derivative of the (elastic) bending moment and the internal shear stresses balance the torque per unit length exerted by the PFR on the Ax, which is given by the D p -dependent term appearing in the first equation. The torque depends on the integral of the bending vector, making the balance equations non-standard (integrodifferental instead of differential). This dependency is due to the fact that the torque arises from the shear deformations of the PFR, which are induced by the shear of axonemal MTs, which is, in turn, related to axonemal bending strains via the integral relations (6) . The torque exerted by the PFR on the Ax is sensitive to the direction given by the unit vector e p , hence it depends on the angle φ p between the Ax-PFR joining line and the unit vector d 1 .

We derive here the dynamic equations for a flagellum beating in a viscous fluid. We consider the extended functional

9/24

where Λ is the Lagrange multiplier vector enforcing the constraint ∂ s r a = d 3 . We treat the fluid-flagellum interaction in the local drag approximation of Resistive Force Theory, see e.g. [35] . In this approximation, viscous forces and torques depend locally on the translational and rotational velocity of the flagellum, represented here for simplicity by the translational and rotational velocity of the Ax. The external viscous forces F and torques G (per unit length) acting on the flagellum are given by

where µ ⊥ , µ || , and µ r are the normal, parallel, and rotational drag coefficient (respectively), and Id is the identity tensor. The principle of virtual work imposes

for every variation δr a and δθ = δθ

, δθ 2 = (δd 3 · d 1 ), and δθ 3 = (δd 1 · d 2 ). Linearizing the force balance equations derived from (16) we obtain the following equations for bending strains and twist

and

which are decoupled from the extra unknown Λ. Equations (17) and (18) are complemented by the boundary conditions

The details of the derivation of (17)- (19) are provided in Appendix B. Once we solve for U 1 , U 2 , and U 3 either the equilibrium equations (12)- (13) or the dynamic equations (17)- (19) , the shape of the flagellum can be recovered. In particular, we obtain the orthonormal frame d i with i = 1, 2, 3 by solving while the centerline of the Ax is recovered by integrating ∂ s r a = d 3 .

We analyze the geometry of the centerline r a which, due to the slenderness of the flagellar structure, is a close proxy for the shape of the flagellum.

In general, the shape of a curve is determined by its curvature κ and torsion τ . Since r a is parametrized by arc length, the two quantities are given by the formulas ∂ s t = κn and ∂ s b = −τ n, where t = ∂ s r a , n = ∂ s t/|∂ s t|, and b = t × n are the tangent, normal, and binormal vector to the curve r a , respectively. Given κ and τ , r a is uniquely determined up to rigid motions.

From the previous definitions and from (3) we obtain the relations between curvature, torsion, bending, and twist. In compact form these relations are given by

which hold for U = 0. In (20) we introduced the angle ψ that the bending vector U = (U 1 , U 2 ) forms with the line U 2 = 0, see Figure 6 . Now, at equilibrium (12) we have

10/24 under any dynein actuation. In other words, axonemal deformations are twistless. This is, fundamentally, a consequence of the fact that shear of axonemal MTs and twist are uncoupled (6) . The torsion of the centerline r a is our main focus, since we are interested in emergent non-planarity. Combining (20) and (21) we have that torsion can arise only from the rotation rate ∂ s ψ of the bending vector U along the length of the flagellum. This last observation will be important in the following.

Under the assumptions (21) and (23), the flagellar energy (1) can be rewritten as

are the target bending strains generated by the dynein forces. The use of this terminology is clear from (22) . The effect of dynein actuation at equilibrium (when the energy is minimized) is to bring the bending strains U 1 and U 2 as close as possible to U * 1 and U * 2 , respectively. The emerging bending strains and the target bending strains might not match due to the interference by the PFR component of the energy (D p = 0). From the formulas for the target bending strains in (22) we can infer the importance of the concentrated shear forces H 1 and H 2 . Without these forces, admissible spontaneous configurations of the Ax would be ruled out. If the concentrated shear forces are null, for example, the Ax cannot spontaneously bend into a circular arc. Indeed, for a circular arc of radius 1/K on the plane d 1 -d 3 we must have U * 2 = 0 and U * 2 = K. In this case, from (22) we have that −H 1 /B a = ∂ s U * 2 = 0, which implies H 1 = 0 and H 1 = −B a K, so the concentrated forces must be non zero.

As per our main hypothesis, we suppose that the sliding inhibition exerted by the bonding links let dyneins organize so that the Ax locally bends spontaneously on the plane d 1 (s)-d 3 (s), as shown in Figure 1 . This is equivalent to require that U * 1 = 0, which leads to the following condition on the shear forces

We consider here the equilibrium equations (12) under the hypothesis (23) . We look at the equilibrium configurations for every possible value of the angle φ p between the Ax-PFR joining line and the spontaneous bending plane of the Ax, even though the value of actual interest for Euglena is φ p ≈ −2π/9. We can prove analytically the following statement: if the Ax-PFR joining line is neither parallel nor orthogonal to the spontaneous bending plane of the Ax, then the emergent flagellar shapes are non-planar. Indeed, suppose H 1 = 0. From (20) and (21) it follows that the shape of the flagellum is planar τ = 0 if (and only if) the angle ψ of the bending vector U(s) = (U 1 (s), U 2 (s)) is constant. The bending vector must therefore be confined on a line for every s. In this case there must be two constants c 1 and c 2 such that U 1 (s) = c 1 U (s) and U 2 (s) = c 2 U (s) for some scalar function U . Now, if a planar U is a 11/24 Figure 5 . Flagellar non-planarity arising from structural incompatibility. The Ax-PFR mechanical interplay is explained in a three-steps argument (left to right). Consider first the two separated structures in their spontaneous configurations (left). The Ax is bent into a planar arc while the PFR is straight. Then, the PFR is forced to match to the Ax, while the latter is kept in its spontaneous configuration (middle). The attachment constraint induces shear strains in the PFR. The composite system cannot then be in mechanical equilibrium without external forcing. When the composite system is released (right), it reaches equilibrium by the relaxation of the PFR shear, which induce additional distortion of the Ax. At equilibrium, an optimal energy compromise is reached, which is characterized by an emergent non-planarity.

solution of (12) then we must have

with c 1 U (L) = 0 and c 2 U (L) = − H 1 /B a . If φ p = 0, π/2, π, 3π/2 the system of equations (24)- (25) admits no solution. Indeed, suppose first that H 1 = 0. Since H 1 = 0 we must have (c 1 , c 2 ) = (0, 0). However, in this case, (24) admits the unique solution U = 0, which is incompatible with (25) . If H 1 = 0 then the boundary conditions impose c 1 = 0, but in this case (24) has again U = 0 as a unique solution, which is incompatible with both the boundary conditions and with (25) . Our statement is thus proved. For φ p ≈ −2π/9, the characteristic value for Euglena, the non-planarity of flagellar shapes is not just possible. It is the only possible outcome under any non-trivial dynein actuation.

Alongside the previous analysis there is a less technical way to infer the emergence of non-planarity from our model. We look here more closely to the flagellum energy, and we think in terms of structural incompatibility between Ax and PFR, seen as antagonistic elements of the flagellum assembly, see Figure 5 .

Under the assumptions (21) and (23), the flagellar energy is given by with U * 2 as in (22) . The energy has two components, W a that depends on the Ax bending modulus B a , and W p that depends on the PFR shear modulus D p . We can vary these material parameters and explore what the resulting minima of W, i.e. the equilibrium configurations (12)-(13), must look like. We consider the nondimensional parameter ν = D p /(B a L −2 ). When ν 1 the Ax component W a of the energy dominates. In this case, at equilibrium, the bending vector has to be close to the target bending vector U ≈ (0, U * 2 ). In particular, then, U(s) will be confined near the line U 1 = 0 for every s. In the case ν 1 the PFR component W p dominates, and the energy is minimized when U(s) lies close to the line generated by the vector e ⊥ p = (− sin φ p , cos φ p ). Clearly, if the latter line is different from U 1 = 0, the two extreme regimes ν 1 and ν 1, each of which favours one of the two individual components, aim at two different equilibrium configurations. In other words, Ax and PFR are structurally incompatible.

When neither of the two energy components dominates, the emergence of non-planarity can be intuitively predicted with the following reasoning. In the intermediate case ν ∼ 1 we expect the equilibrium configurations to be a compromise among the two extreme cases, with the bending vector U(s) being "spread out" in the region between the two extreme equilibrium lines. The spreading of the bending vector is aided by the concentrated shear force at the tip, which imposes U(L) = (0, U * 2 (L)) irrespectively of the PFR stiffness. The bending vector is then "pinned" at s = L on the U 1 = 0 line while it gets dragged toward the line generated by e ⊥ p for large values of ν. Hence the spreading. The bending vector will then span an area and, consequently, undergo rotations. Since torsion is determined by the rotation rate of the bending vector τ = ∂ s ψ, the resulting flagellar shapes will be non-planar. Figure 6 illustrates a critical example in which the previous intuitive reasoning effectively plays out. We consider a target bending of the kind U * 2 (s) = A 0 + A 1 sin(2πs/L)), a fair idealization of the asymmetric shapes of a Chlamydomonas-like flagellar beat [10] . We take φ p = −π/4 (larger than the Euglena value, to obtain clearer graphs). For ν = 0 the bending vector lies inside the U 1 = 0 line, and its amplitude oscillates. For positive values of ν, when the PFR stiffness is "turned on", the oscillating bending vector is extruded from the U 1 = 0 line. For large values of ν it gets closer and closer to the line generated by e ⊥ p . The bending vector spans an area and, following the oscillations, it rotates clock-wise and anti-clock-wise generating an alternation in the torsion sign. This is the geometric signature of the spinning lasso.

Our model is able to predict the torsional characteristic of the euglenid flagellum in the static case, and in absence of external forces. We test here the model in the more realistic setting of time-dependent dynein actuation in the presence of hydrodynamic interactions.

We first observe that, as in the static case, the dynamic equations for U and U 3 are decoupled (17)- (18) , and that dynein forces do not affect twist. We have then twistless kinematics under any actuation also in the dynamic case, at least after a time transient. We can simply assume (21) for all times, so the torsion of r a is still completely determined by the bending vector.

We consider a dynein actuation that generates Chlamydomonas-like shapes in a flagellum with no extra-axonemal structures. The shear forces H 1 and H 1 employed in our simulation are shown in Figure 7 . The same figure also shows the emergent bending strains of a PFR-free flagellum actuated by said forces, beating in a viscous fluid. The dynamic equations for this system are simply (17)-(18) with D p = 0. The resulting bending strains, which generate a planar beat, resembles the experimentally observed Chlamydomonas flagellar curvatures reported in [23] .

Finally, Figure 7 presents the emergent bending strains of the beating euglenid (PFR-bearing) flagellum, together with the corresponding flagellar torsion. The spinning lasso torsional signature is clearly present. Indeed, the fluid-structure interaction does not disrupt the Ax-PFR structural incompatibility, which still generates non-planar shapes with travelling waves of torsional peaks with alternating sign and the typical looping-curve outlines, cfr. Figure 2 . All the details on the methods and parameters employed in the simulations are given in Appendix C and Appendix D.

We have shown how the origin of the peculiar shapes of the euglenid flagellum can be explained by the mechanical interplay of two antagonistic flagellar components, the Ax and the PFR. Our conclusions are based mainly on the hypothesis that sliding inhibition by the PFR organizes dynein activity, and localizes the spontaneous bending plane of the Ax as the one that passes from the Ax center through the MTs bonded to the PFR. This is in agreement with the current understanding of the mechanism that generates beat planarity in other PFR-bearing flagellar systems. Non-planarity in Euglena can arise because of a marked asymmetry in the Ax-bonding links-PFR complex in the Euglenid flagellum, which is not found in kinetoplastic organisms such as Leishmania [11] or Trypanosoma [22] . In the absence of a precise knowledge of the dynein actuation pattern we tested our mechanical model under shear forces that would, in the absence of extra-axonemal structures, realize a beat similar to those found in model systems like Chlamydomonas. We appreciate that the emergent distortion of the Ax, generated by the Ax-PFR interplay, could in principle lead to different actuation patterns, consistently with the hypothesis of dynein actuation via mechanical feedback. This question will require further studies.

Along with the mechanism that let the euglenid flagellar shapes emerge, it is worth considering how this characteristic flagellar beat is integrated in the overall behaviour of the organism. As shown in [25] , the spinning lasso beat produces the typical roto-translational trajectories of swimming euglenas. Cell body rotation is in turn associated with phototaxis. Indeed, rotation allows cells to "scan" the environment, and veer to the light source direction when stimulated (or escape in the opposite direction, when the signal is too strong). Here, the key biochemical mechanism could be the one often found in nature, by which periodic signals generated by lighting and shading associated with body rotations are used for navigation, in the sense that the existence of periodicity implies a lack of proper alignment [12] . It is known that the PFR is directly 15/24 connected with the light-sensing apparatus [24] , and might even be active [21] . Further study on euglenid flagellar motility and phototaxis could lead to a more comprehensive understanding on the role of the PFR.

Strain SAG 1224-5/27 of Euglena gracilis obtained from the SAG Culture Collection of Algae at the University of Göttingen was maintained axenic in liquid culture medium Eg. Cultures were transferred weekly. Cells were kept in incubator at 15 • C at a light:dark cycle of 12:12 h under a cold white LED illumination with an irradiance of about 50 µmol · m −2 · s −1 .

An Olympus IX 81 inverted microscope with motorized stage was employed in all the experiments. Experiments were performed at the Sensing and Moving Bioinspired Artifacts Laboratory of SISSA. The microscope was equipped with a LCAch 20X Phc objective (NA 0.40) for the imaging of cells trapped at the tip of a glass capillary using transmitted brightfield illumination. The intermediate magnification changer (1.6 X) of the microscope was exploited to achieve higher magnification. Micrographs were recorded at a frame rate of 1, 000 fps with a Photron FASTCAM Mini UX100 high-speed digital camera.

Tapered capillaries of circular cross section were obtained from borosilicate glass tubes by employing a micropipette puller and subsequently fire polished. At each trial observation a glass capillary was filled with a diluted solution of cells and fixed to the microscope stage by means of a custom made, 3d-printed holder. The holder allowed for keeping the capillary in place and rotating it about its axis, so as to image a cell specimen from distinct viewpoints. Cells were immobilized at the tip of the capillary by applying a gentle suction pressure via a syringe connected to the capillary by plastic tubing.

The Ax consists of a bundle of inextensible filaments of length L (MTs) lying on a cylindrical surface of radius ρ a . For simplicity, the model ignores the mechanical effects of radial spokes and the central pair. The axonemal surface is parametrized by the generalized cylindrical coordinates z and φ via the map

with r a , d 1 , and d 2 defined as in the main text. Following [14] , we suppose that the axonemal constraints confine MTs on the Ax surface at a fixed angular distance ∆φ = 2π/9 between each other. More formally, for j = 1, . . . , 9, we define the centerline r j of the j-th MT as r j (s) = C(s, φ j ), where

The function Z(s, φ) in (28) is defined (implicitly) via the equality

From the definitions above follows ∂ s r j = 1, so MTs are indeed inextensible and s is their arc length. Moreover, the Taylor expansion of C at the first order in ρ a gives the approximated formula (5), with γ 1 and γ 2 given by (6) . We associate to the j-th MT an orthonormal frame along r j given by the unit vectors

where N = cos φ d 1 (Z) + sin φ d 2 (Z) is the (outer) unit normal to the cylindrical surface. The unit vectors (30) determine MTs' cross section orientations. The unit vectors e j 1 (s) and e j 2 (s) lie on the cross-section centered at r j (s), while e j 3 (s) = ∂ s r j (s) is orthogonal to it. The (passive) elastic energy of the Ax is given by the sum of the MTs' elastic energies

are the strains associated to the j-th MT, while B m and C m are the MTs' bending and twisting moduli (respectively). At the leading order approximation in ρ a we have

From (32) and (3) follows that W a pas , as defined in (31), can indeed be approximated by the right hand side of (2), with B a = 9B m and C a = 9C m .

We consider now MT sliding. Fixed a point r j−1 (s) on the (j − 1)-th MT's centerline we look for the nearest point to r j−1 (s) on the centerline r j of the j-th MT. Such a point r j (s * ) (we can think of it as a projection) lies at some arc length s * , which depends on s. We write s * = Π j (s). We then define the sliding σ j (s) as the difference of the two arc lengths s and Π j (s). More formally, Figure A1 illustrates the geometric idea of definition (33) . Figure A1 . Sketch of two MTs' centerlines during deformation. The sliding σ j (s) is defined as the difference between the arc lengths s and Π j (s). The latter is the arc length corresponding to the projection of r j−1 (s) on the curve r j . We have positive sliding when dyneins push the j-th MT towards the distal end of the flagellum and the (j − 1)-th MT toward the proximal end.

The active internal energy of the Ax is defined as minus the total mechanical work of the dyneins

where F j (s) are the sliding forces on the j-th MT exerted by the dyneins on the (j − 1)-th MT, and F j are the singular sliding forces (on the j-th MT exerted by the dyneins on the (j − 1)-th MT) concentrated at the distal end of the Ax. Taylor expanding (33) in ρ a we have, at the leading order,

with γ 1 and γ 2 given by (6) . From (35) we have that W a act , as defined in (34), is approximated by the right hand side of (4), with

(cos φ j − cos φ j−1 )F j (s) , H 2 (s) = ρ a 9 j=1 (sin φ j − sin φ j−1 )F j (s) ,

In the remainder we give some details of the derivation of (11) from (7). Expanding (10) at the leading order in ρ p gives

with γ p⊥ (s) = − sin φ p γ 1 (s) + cos φ p γ 2 (s). From (8) and (37), leading order calculations give V 1 , V 3 ∼ ρ p , whereas V 2 ≈ γ p⊥ / 1 + γ 2 p⊥ . Linearizing in γ p⊥ we have V 2 ≈ γ p , from which follows (11).

To derive equations (17)- (19) it is convenient to introduce the quantities M 1 , M 2 , and M 3 defined via the following variational equality 

where γ p⊥ (s) = − sin φ p γ 1 (s) + cos φ p γ 2 (s), as in the previous section. We then write the variations δU i in terms of δθ i (defined in the main text) obtaining     

Combining (38) In the calculations above we took variations with δr a (0) = δθ(0) = 0, since we consider a flagellum with a clamped end at s = 0. The principle of virtual work (16) gives us then the following force and torque balance equations

with Λ(L) = 0 , and M(L) = 0 .

Equations (17)- (19) are derived from (42) and (43), after some extra formal manipulations that we explain in the reminder. We first introduce the local angular velocities

which are related to the strains via the following compatibility equations  

Equations (45) follow from the identities ∂ s ∂ t d i = ∂ t ∂ s d i , with i = 1, 2, 3. We then rewrite the external forces and torques (15) in compact form as 

We rewrite the force balance equation as ∂ t r a = V [d 3 ] −1 ∂ s Λ and, after differentiating both sides with respect to s and then scalar multiplying by d 1 , d 2 , and d 3 , we have

From (40), (45), and (48) we obtain (18) by first differentiating with respect to s both sides of (48), and then by linearizing the resulting equation. Similarly, exploiting also equations (47) this time, we differentiate with respect to s and then linearize (49) to obtain (17) . Boundary conditions (19) follow from (43) and

which follows from the fixed end condition r a (0, t) = 0. Given U 1 , U 2 , and U 3 the orthonormal frame d 1 , d 2 , and d 3 can be recovered by integrating the system of equations

which is derived from (3). We can then recover the centerline r a by integrating ∂ s r a = d 3 . The PFR centerline r p and the orthonormal frame g 1 , g 2 , and g 3 follow from (9) and (10).

We define the nondimensional variables for arc length x, time y, strains u i , shear forces h i , and concentrated shear forcesĥ i as follows x = s/L , y = t/T , u i = U i /L −1 , h i = H i /B a L −2 , andĥ i = H i /B a L −1 .

The equilibrium equations (12) , with (21) and (23), are solved by seeking for a minimizer of the energy (26) in the adimensional form

In the formula above u = (u 1 , u 2 ) is the adimensional bending vector and u * 2 = U * 2 /L −1 is the adimensional target bending. We find the minimizer using the gradient descend method u n+1 = u n − α δw δu 

and α is a conveniently chosen step size. We initiate the algorithm with u 0 = (0, u * 2 ), and then we iterate (54) until u n+1 − u n falls below a pre-set tolerance parameter.

The dynamic equations (17) are recast, and then solved, in terms of the shear vector γ = (γ 1 , γ 2 ) where γ 1 (x) = x 0 u 2 , and γ 2 (x) = − x 0 u 1 .

(55)

Simulating the complex cell design of trypanosoma brucei and its motility

Geometry and mechanics in the opening of chiral seed pods

Euglena gracilis a model for flagellar surface assembly, with reference to other cells that bear flagellar mastigonemes and scales

Movement and locomotion of euglena

Thinking about flagellar oscillation

The paraflagellar rod: a structure in search of a function

Antagonistic behaviors of nmy-1 and nmy-2 maintain ring channels in the c. elegans gonad

Flagellar and ciliary beating in trypanosome motility

Independent control of the static and dynamic components of the chlamydomonas flagellar beat

Flagellum assembly and function during the leishmania life cycle

Green algae as model organisms for biological fluid dynamics

The mathematics and mechanics of biological growth

The chirality of ciliary beats

Morphomechanical innovation drives explosive seed dispersal

The diversity of hydrostatic skeletons

Motility and more: the flagellum of trypanosoma brucei

One of the nine doublet microtubules of eukaryotic flagella exhibits unique and partially conserved structures

Flagellar and ciliary beating: the proven and the possible

Flagellar membrane specializations and their relationship to mastigonemes and microtubules in euglena gracilis

Atpase activity in flagella from euglena gracilis. localization of the enzyme and effects of detergents

The paraflagellar rod of kinetoplastid parasites: from structure to components and function

Flagellar kinematics and swimming of algal cells in viscoelastic fluids

Ultrastructure of the apical zone of euglena gracilis: photoreceptors and motor apparatus

Kinematics of flagellar swimming in euglena gracilis: Helical trajectories and flagellar shapes

Swimming like algae: biomimetic soft artificial cilia

Curvature regulation of the ciliary beat through axonemal twist

Dynamic curvature regulation accounts for the symmetric and asymmetric beats of chlamydomonas flagella

Coronavirus envelope protein: current knowledge

Polygonal motion and adaptable phototaxis via flagellar beat switching in the microswimmer euglena gracilis

Anisotropic escrt-iii architecture governs helical membrane tube formation

A comparison of paraxial rods in the flagella of euglenoids and kinetoplastids

On the unity and diversity of cilia

Use of chiral cell shape to ensure highly directional swimming in trypanosomes

Twirling and whirling: Viscous dynamics of rotating elastic filaments

Evidence for a sliding-resistance at the tip of the trypanosome flagellum

Flexural rigidity and shear stiffness of flagella estimated from induced bends and counterbends

We thank A. Beran for his assistance with E. gracilis samples. This study was supported by the European Research Council through the ERC Advanced Grant 340685-MicroMotility.

After solving for γ, we obtain the bending strains by differentiation. The equation for γ isand h = (h 1 , h 2 ). The corresponding boundary conditions are given bywhereĥ = (ĥ 1 ,ĥ 2 ). The above (57) are point-wise conditions, which do not involve integral terms as in (19) . Avoiding this non-locality allows for an easier numerical implementation. The finite difference scheme we employ to solve (56)-(57) is illustrated in the remainder. We consider the discrete time sequence y n = n∆y with n = 0, 1, 2, . . . and we define γ n (x) = γ(y n , x), h n (x) = h(y n , x), andĥ n =ĥ(y n ). Equation (56) is discretized in time with the one-step (semi-implicit) numerical scheme η ∆yand complemented by the boundary conditionsThe adimensional arc length interval [0, 1] is discretized uniformly in M + 1 points x k = k∆x , with k = 0, 1, . . . , M = 1/∆x. We also consider the extra "ghost points"The discrete values ∂ 2 x γ n k = ∂ 2 x γ n (x k ) of the second derivative in (58) are approximated by the finite difference schemewhere γ n k = γ n (x k ). Analogous formulas are employed for the second derivative of h n . The discrete values ∂ 4x γ n+1 k = ∂ 4 x γ n+1 (x k ) of the forth derivative in (58) are given by the schemewhich involves the ghost points values γ −1 , γ M +1 , and γ M +2 . The discretized approximations of the boundary conditions (59) are given by

The previous formulas give us the expressions for the ghost points' values γ n+1 −1 , γ n+1 M +1 , and γ n+1 M +2 in terms ofĥ n+1 , γ n+1 k , γ n k , and h n k with 0 ≤ k ≤ M . These expressions are then plugged in (61). In turn, the iterative scheme (58) allows to calculate γ n+1 k from γ n k with 0 ≤ k ≤ M , while incorporating the boundary conditions (59) in the numerical solution.The scheme is iterated for several time periods until a periodic solution is reached.

We obtain the history of shear forces presented in Figure 7 by solving the following inverse dynamical problem. We assign first a history of normalized bending strains (0, u 2 (x, y)), periodic in time, that imitates the experimentally observed Chlamydomonas flagellar curvatures reported in [23] . We use the following modelwith λ(x) = 2πλ 0 (x + λ 1 sin(πx)) and ω(y) = 2π(y + ω 1 sin(πy)) .Then, we calculate the shear forces that generate said history of bending strains for an Ax beating in a viscous fluid (without extra-axonemal structures attached to it). Equation (56) with ν = 0 and h 2 = h 2 = 0 defines exactly the dynamics of this system. We can solve for h 1 and h 1 explicitly, obtaining h 1 (y) = −u 2 (1, y) andIn the dynamic simulation in Figure 7 we used the following numeric values for the physical parameters of the system. The bending modulus of the Ax is B a = 840 pN · µm 2 , taken from [37] . We set L = 28 µm, T = 25 ms, and µ ⊥ = 3.1 fN · s · µm −2 , which are all values estimated in [25] . The angle between spontaneous bending plane and the Ax-PFR joining line φ p = −2π/9 is estimated from micrographs in [20] and [4] . Without direct measurements for D p , we set ν = 20 as the value giving the most balanced mechanical interplay between Ax and PFR (the value of ν with the most evenly spread-out bending vector's solutions in Figure 6 ). For the bending strains parameters in (63) we took a 0 = 7.8, a 1 = 7.5, λ 0 = 1.85, λ 1 = 0.1, and ω 1 = −0.1.