Effect of fluid elasticity on the emergence of oscillations in an active elastic filament

Many microorganisms propel themselves through complex media by deforming their flagella. The beat is thought to emerge from interactions between forces of the surrounding fluid, the passive elastic response from deformations of the flagellum and active forces from internal molecular motors. The beat varies in response to changes in the fluid rheology, including elasticity, but there are limited data on how systematic changes in elasticity alter the beat. This work analyses a related problem with fixed-strength driving force: the emergence of beating of an elastic planar filament driven by a follower force at the tip of a viscoelastic fluid. This analysis examines how the onset of oscillations depends on the strength of the force and viscoelastic parameters. Compared to a Newtonian fluid, it takes more force to induce the instability in viscoelastic fluids, and the frequency of the oscillation is higher. The linear analysis predicts that the frequency increases with the fluid relaxation time. Using numerical simulations, the model predictions are compared with experimental data on frequency changes in the bi-flagellated alga Chlamydomonas reinhardtii. The model shows the same trends in response to changes in both fluid viscosity and Deborah number and thus provides a possible mechanistic understanding of the experimental observations.


I. INTRODUCTION
Many microorganisms, such as sperm, propel themselves through complex media by deformations of their flagella.It has long been observed that the rheology of the surrounding fluid alters the shape and frequency of the flagellum beat of mammalian sperm [1][2][3][4][5], of sea urchin and related marine animal sperm [6,7], and of the bi-flagellated alga Chlamydomonas reinhardtii [8,9].In addition to beat changes, fluid viscosity has been show to affect the coordination in arrays of cilia that drive cell locomotion [10] and transport mucus [11].
Human sperm and sperm from marine invertebrates exhibit different gait changes in response to high viscosity environments [5].It has been hypothesized that the gait changes in mammalian sperm are important for fertilization because they must swim through viscoelastic mucus [2,5].There have been multiple observations of sperm gaits in mucus and other viscoelatic fluids [1][2][3][4], but there has been no systematic documentation of how gradual changes in fluid elasticity affect the gait.Changes in the beat frequency, shape, and swimming speed of Chlamydomonas reinhardtii in response to both viscous and elastic properties of the surrounding fluid were recently documented [8]; it was observed that the beat frequency was enhanced by fluid elasticity, and the frequency changed nonmonotonically with fluid viscosity in viscoelastic fluids.
Fluid elasticity clearly influences the flagellum beat, but the physical mechanism for how fluid rheology shapes the beat is not known.There have been theoretical studies of how fluid elasticity affects the gait for prescribed active motor forces [12,13].These studies considered the active forces as a traveling wave with a given frequency, and they examined how the fluid rhelology affected the resulting shapes.This approach was able to explain the qualitative shape change observed in some sperm species in viscoelastic fluids [3], and the increased amplitude of the beat, and thus increased swimming speed, in artificial swimmers with flexible tails [14].However because the frequency of the active forces was prescribed, this approach cannot be used to understand how the beat frequency changes with fluid rheology as observed in [8].
The flagellum beat is powered by dynein motors that form crossbridges between the microtubule doublets that make up the axoneme.These motors generate active shear forces between adjacent doublets that through interactions with other passive forces and constraints at the base leads to bending [15].It is not understood how the motors along the flagellum are coordinated spatially and temporally to produce the observed waves of bending.There are different hypotheses about how mechanical feedback on motor activity from deformations of the flagellum lead to emergent coordination of the whole system.Some of the leading feedback mechanisms that have been explored are that the motors respond to changes in curvature [16][17][18], tangential deformations (i.e.sliding control) [19,20], or normal forces (i.e."geometric clutch") [21,22].All of these mechanisms have been explored thoroughly in models, and they are all capable of producing emergent waves in the flagellum.Analysis of the bifurcation structure of the three models cast doubt on the sliding-control mechanism [23], though sliding control was capable in matching experimental data on bull sperm [20].A comparison of all three models on data on Chlamydomonas reinhardtii data favors curvature control [24].Despite years of theoretical effort, it is not clear if any of these feedback mechanisms are involved in producing the flagellum beat.
In [25] an alternative mechanism was proposed and analyzed for producing the flagellum beat that does not require dynein regulation nor spatiotemporal organization of dynein activity.The mechanism suggested in [25] is related to a dynamic instability known as flutter that results when an elastic structure in fluid is subject to axial loading.Dynein generates a tension that buckles the filament.Unlike static buckling, the direction of the motor forces remains tangent to the filament as it deforms which results in an oscillatory motion.
The instability analyzed in [25] is similar to the instability of filaments under external load in the tangent direction known as a "follower force" (because it follows the direction of the filament) [26].There have been several recent analyses of filaments subject to a follower force at low Reynolds number inspired by the motion of biological filaments driven by molecular motors [27][28][29][30][31][32][33].These models have been used to understand observations in in vitro motility assays [27], cytoplasmic streaming [31], beating flagella and cilia [29], and coordination between pairs [32] and arrays of cilia [33].These works have thoroughly analyzed the Hopf bifurcation from rest to a beating pattern including how different boundary conditions and restrictions (i.e.planar vs. 3D) result in different dynamics in viscous fluids [28][29][30].
In this paper we analyze the emergence of oscillations of a planar elastic filament pinned at the base subject to a follower force at the tip in a viscoelastic fluid.We examine how the elasticity of the surrounding fluid affects the strength of the applied force needed to produce the oscillation and the emergent frequency at the bifurcation.Our results show that the critical value of the force at which oscillations occur is greater in a viscoelastic fluid than in a Newtonian fluid, and the frequency of the beat is always increased by the elasticity of the fluid.We compare the model predictions for how the frequency changes with relaxation time and total viscosity with the experimental measurements from [8].Our analysis captures the observed frequency increases with relaxation time and the nonmonotonic frequency response to changes in viscosity, and thus offers a possible mechanistic explanation for how the beat frequency is affected by fluid elasticity.

II. MODEL AND EQUATIONS
We consider the motion of a slender, inextensible, planar, elastic filament, clamped at one end and subject to a compressive follower force of strength Γ at the tip of the filament in a viscoelastic fluid at zero Reynolds number.The mathematical model is analogous to that presented in [28] with the addition of fluid viscoelasticity.Let 0 ≤ s ≤ L be the arclength where s = 0 corresponds to the clamped base.The position of the filament is X(s, t) = (x(s, t), y(s, t)); see Figure 1.
The instantaneous force balance for the filament is where the first term is the force per unit length from bending, the second term represents the tension that enforces the inextensibility constraint |X s | = 1, and F fluid is the drag force from the surrounding fluid.We assume that the fluid drag can be expressed as the drag in a viscous fluid plus a drag accounting for the viscoelastic effects.One can think of the fluid as composed of a Newtonian solvent with the addition of polymers which are responsible for viscoelastic stresses.Thus, the drag force is where F sol and F pol are the drag forces due to the solvent and polymers, respectively.
< l a t e x i t s h a 1 _ b a s e 6 4 = " / + j Y J u c / + y t Q T I w 7 b g 8 t < l a t e x i t s h a 1 _ b a s e 6 4 = " e Y w C c h 9 u T g 5 q X q 0 x 5 G 4 7 1. Schematic of a horizontal flexible filament clamped at one end with a follower force of strength Γ applied at its tip.The filament position is defined as X(s, t), where 0 ≤ s ≤ L is the arc length coordinate.The local tangent vector is Xs = t(s, t).
At the clamped end, we have the boundary conditions X(0, t) = 0, and X s (0, t) = êx , while at the free end which capture the fact that the filament is torque-free and that the force at the tail and the external force must balance.The compressive follower force, ΓX s , is applied tangentially to the tip of the filament.This non-conservative force drives the motion of the filament.
The viscous drag force due to the solvent F sol acting on the filament from the surrounding flow is given by resistive force theory [34] which provides a local relation between the local filament velocity, X t , and the hydrodynamic force per unit length.The viscous drag force per unit length is defined as where t and n are the local tangent and normal unit vectors.The drag coefficients in the perpendicular and parallel direction are ξ s ⊥ and ξ s || are proportional to the solvent viscosity µ s .For example, for a cylinder of radius b and length L, ξ s ⊥ = µ s α, where α = 4π/[ln(L/b) + 1/2] and ξ s ⊥ /ξ s || → 2 as L/b → ∞ [34].
We are interested in the effects of viscoelasticity at and near the bifurcation at which oscillations first emerge, and are thus small in amplitude.Given this, we utilize a linear viscoelastic model to describe the polymeric force F pol : where ξ p ⊥ = µ p α, µ p is the polymeric viscosity, and τ is the fluid relaxation time.In the limit of zero fluid relaxation time (τ = 0) or zero polymer viscosity, (µ p = 0), the polymeric force F pol = 0 and we recover the viscous Newtonian fluid.The assumption of a linear viscoelastic model for the drag force on a deforming filament at small amplitude was utilized by previous authors to analyze the effect of viscoelasticity on shape changes of flagellum shapes and bending filaments [12,13].Further, this assumption was numerically validated in [13] by comparing this linear model with numerical simulations that involve the nonlinear viscoelastic stress.
Equations ( 1), (6), and ( 7) are nondimensionalized by rescaling lengths by L, time by the viscous time scale L 4 (ξ s ⊥ + ξ p ⊥ )/k b , tension (the Lagrangian multiplier) by k b /L 2 , and polymeric force by Here is the ratio of tangential and normal drag coefficients, and is the dimensionless relaxation time.Note that the ratio of tangential and normal drag coefficients depends only on the aspect ratio of the filament, not on the viscosity, and thus has the same value for the solvent and polymerirc fluid.
Similarly rescaling equations ( 4) and ( 5) yields the following dimensionless free boundary conditions where the dimensionless ratio between the strength of the force at the tip and the elastic force is defined as Since the force is compressive (Γ > 0), σ is always positive.

III. BIFURCATION ANALYSIS IN A VISCOELASTIC FLUID
In a viscous fluid, the strength of the follower force, σ, is the only nondimensional parameter.As analyzed in [28], there is a critical strength of the follower force, σ 0 , below which the straight filament at rest is stable.At σ = σ 0 there is a supercritical Hopf bifurcation so that for σ > σ 0 the filament exhibits sustained oscillations.In a viscoelastic fluid there are three dimensionless parameters to consider: the follower force strength, σ, the relaxation time, λ, and the viscosity ratio, β.In this section we analyze how the critical follower strength and the frequency of the emergent oscillation in a viscoelastic fluid depend on λ and β.
We consider small amplitude deviations from the rest state of a straight filament.In the small deformation regime, the tension to leading order is constant and therefore equal to the external force applied to the tip [19,28], so that T (s, t) = σ.The leading order equations for the vertical displacement, y(s, t), and vertical component of the polymer force, f pol (s, t), are Deformations in the horizontal direction and changes to the tension occur at higher order in deformation.The resulting boundary conditions are A. Relationship to Stability in a Viscous Fluid We assume solutions of the form y(s, t) = ŷ(s)e ηvet and f pol (s, t) = f pol (s)e ηvet in equations ( 13)- (15).The real part of η ve quantifies the growth (or decay if negative) rate of perturbations in viscoelastic fluid, and its imaginary part gives the frequency of oscillations.Eliminating f pol we obtain We let L = −∂ ssss − σ∂ ss denote the operator acting on the space of functions that satisfy boundary conditions (15).We express (16) as the eigenvalue problem where denotes the eigenvalues of L, the real part of which quantifies the growth rate of perturbations in a viscous fluid.Consistent with this idea, note that when either β = 1 or λ = 0, η ve = η v because the viscoelastic fluid reduces to a viscous fluid.

Instability in a Viscous Fluid is Necessary for Instability in a Viscoelastic Fluid
We first show that for a given strength of the follower force, the system is unstable in a viscoelastic fluid only if the system is unstable in a viscous fluid.Equation ( 18) relating the eigenvalues of the follower problem in the viscoelastic fluid to those in a viscous fluid can be expressed as The real parts of η v and η ve are related by where α 0 and α 1 are real valued, nonnegative functions of η ve .Therefore, if Re(η ve ) > 0, then Re(η v ) > 0. This establishes that instability in a viscous fluid is a necessary condition for instability in a viscoelastic fluid.

More Force Required for Instability in Viscoelastic
Let σ denote the value of the follower force at which the rest state become unstable in the viscoelastic fluid.Because Re(η ve ) = 0 when σ = σ, Re(η v ) = c 0 (η ve ) > 0 from (20).Because Re(η v ) < 0 for all σ < σ 0 , it follows that σ > σ 0 .Therefore it follows that the follower force required for instability in the viscoelastic case is always larger than the force required in the viscous case.

Higher Frequency in Viscoelastic Fluid
We show that for the same follower force, the frequency of the oscillation in the viscoelastic fluid is always larger the frequency in the viscous fluid.From equation (19), the imaginary parts of η v and η ve are related by Assume that Re(η ve ) > 0 so that the rest state is unstable.If follows that |1 + λη ve | 2 > 1, and thus from ( 21) Because the emergent frequency of the oscillation near the bifurcation is approximately the imaginary part of the eigenvalue, we conclude that viscoelasticity increases the frequency of the oscillation.

B. Numerical Calculation of Eigenvalues
We obtain the viscous eigenvalues using a second-order, centered finite-difference discretization of the operator L appropriately modified near the ends to account for the boundary conditions.The real and imaginary parts of the eigenvalue with largest real part are plotted in Figure 2(a,b).Consistent with [28] we find the critical strength of the follower force at which oscillations emerge in the viscous fluid is σ 0 ≈ 37.7.At this critical value, the eigenvalues corresponding to the bifurcation are η v ≈ ±191i.We define ω 0 ≈ 191 as the angular frequency of the emergent oscillation in a viscous fluid.Given the viscous eigenvalues as a function of σ, for a given relaxation time, λ, and viscosity ratio, β, we identify the corresponding viscoelastic eigenvalues by solving equation ( 18) for η ve .To identify instability in the viscoelastic case, we only need to compute the viscoelastic eigenvalues when the corresponding viscous eigenvalues have positive real part.We find that for σ 0 < σ < 174.6 there is only the single pair of eigenvalues with positive real part for the viscous fluid as shown in Figure 2(c).In the remainder of this work we restrict σ < 174.6 which simplifies the stability analysis in the viscoelastic case because we only need to consider the viscoelastic eigenvalues corresponding to a single viscous eigenvalue.

Limit of Large Relaxation Time
As λ → ∞ equation ( 18) at leading order is Thus in this limit of large relaxation time the viscoelastic eigenvalue is proportional to the viscous eigenvalue, and two conclusions follow immediately from this relation.First, the critical follower force strength at which oscillations emerge in the viscoelastic fluid (i.e. the Hopf bifurcation point) is identical to the critical follower force for a viscous fluid.Second, the emergent frequencies in the two fluids are related by Because 0 < β ≤ 1, the frequency in the viscoelatic case is always larger than the frequency in the viscous case.In summary, in the limit λ → ∞, the bifurcation occurs at the same follower force strength and the emergent frequency is greater in the viscoelastic case by a factor inversely proportional to the viscosity ratio.We note that many biologically relevant media such as respiratory and cervical mucus have relatively large relaxation times [35].

Limit of vanishing polymer viscosity: β → 1
There are two limits in which the viscoelastic fluid reduces to a viscous fluid: vanishing relaxation time (λ → 0) and vanishing polymer viscosity (β → 1).Both limits are useful in considering perturbation from a viscous fluid.We consider the limit β → 1, which physically corresponds to the limit of small polymer viscosity and is typical of dilute polymeric solutions.In this limit we are able to examine how the critical follower force strength and emergent frequency depend on both the relaxation time and follower force strength for all relaxation times.As we explain later, this analysis also captures the limit λ → 0.
We next consider the limit β → 1, which physically corresponds to the limit of small polymer viscosity and is typical of dilute polymeric solutions.In this limit we are able to examine how the critical follower force strength and emergent frequency depend on both the relaxation time and follower force strength.
At the bifurcation point the viscoelastic eigenvalue is pure imaginary, i.e. η ve = iω, where ω is the angular frequency at the bifurcation point.The relationship between the viscous and viscoelastic eigenvalues, equation (18), at the bifurcation point is where the notation η v (σ) is used to emphasize that the viscous eigenvalue depends on the unknown follower force strength, σ.This is a single complex valued equation involving the two real-valued unknowns σ and ω.
For β close to 1, σ will be close to σ 0 , the critical follower force strength in a viscous fluid.We linearize η v about this point so that where σ = σ − σ 0 and (dη v /dσ)| σ0 = a + bi.Eliminating η v from (S6) using ( 26) and then equating the real and imaginary parts results in where ϵ = 1 − β.We seek a solution in the limit ϵ → 0 by using the expansions σ = ϵσ 1 + ϵ 2 σ 2 . . .and ω = ω 0 + ϵω 1 + . . .and matching the O(ϵ) terms.The resulting expansions for critical follower force and corresponding frequency are In these expressions, the relaxation time appears paired with the frequency in the product λω 0 .This quantity is similar to the Deborah number, but the frequency is fixed at the emergent frequency in the viscous limit.In this section we consider how the bifurcation location and emergent frequency depend on the scaled relaxation time λω 0 .Before continuing, we remark that one can recover the λ → 0 expansions from these expressions.Equations ( 27) and ( 28) are valid near the viscous bifurcation point, and thus when λ → 0 for all β.Equation (29) follows directly from (27).The expressions for the frequency in the limit λ → 0 that one obtains from expanding (28) or (30) are equivalent at O(λ).
In Figure 3 we plot the asymptotic and numerical solutions for critical follower force and emergent frequency in a viscoelastic fluid relative to their respective values in a viscous fluid for β = 0.95 as a function of relaxation time.Plots for other values of β are shown in the Supplementary Information.As the relaxation time increases, the critical force increases and then decreases.The peak critical force occurs at λω 0 = 1, and as expected from the large relaxation time analysis σ → σ 0 as λ → ∞.
For small polymer viscosity, there remains a single critical value of the follower force above which oscillations emerge.However, for the range of follower forces below the peak in Fig. 3(a) (which is approximately 1 < σ/σ 0 < 1 + (1 − β)/(2aσ 0 )), there are two Hopf bifurcation points as the relaxation time changes.For a follower force in this range, the filament oscillates at low and high relaxation times while the rest state is stable for an intermediate range of relaxation times.The frequency at the higher relaxation time is generally greater than that of the lower relaxation time.For example, for β = 0.95 and σ = 1.01σ 0 the two bifurcation points and critical frequencies are marked with red dots in Figure 3.In later sections, we examine how the frequency changes with relaxation time for fixed follower force, and we will see that generally the frequency increases with increasing relaxation time.

D. Bifurcation Location and Emergent Frequency for General β
We examine the bifurcation for general values of the viscosity ratio β by solving equation ( 18) for η ve using the numerically computed viscous eigenvalues.Note that there are two values of η ve for each viscous eigenvalue.We found that for all parameter regimes we explored, one of the two values of η ve always had negative real part.As discussed in the Supplementary Information, this additional eigenvalue scales with 1/λ and is likely related to the fluid relaxation time scale.In Figure 4 we show the location of the bifurcation in the σ-λω 0 plane for β = 0.1, 0.2, . . .0.9.For many values of β the shape of the curve denoting the location of the bifurcation is qualitatively similar to that of the asymptotic result.In the limit β → 1, the peak in the critical follower force occurs at λω 0 = 1, but as β decreases the relaxation time corresponding to this peak decreases.
For example for β = 0.5, the peak occurs near λω 0 ≈ 0.72.For small values of β the critical force strength is no longer a single valued function of the relaxation time.The inset in Figure 4 shows the bifurcation location for β = 0.1 for a smaller range of λ to highlight this feature.At β = 0.1 for 0.1749 < λω 0 < 0.2955 (end points marked with gray lines in the figure), there are three σ values at which bifurcations occur.
The quantity λω 0 that appears in the asymptotic expressions ( 29)- (30) represents the Deborah number in the limit of vanishing polymer viscosity (β → 1).More generally, we take as the Deborah number the product of the relaxation time and the emergent frequency; i.e.De = λω.In Figure 5 we examine how the critical force and the emergent frequency depend on λω 0 and De = λω for β = 0.5 and β = 0.1.For β = 0.5 the shapes of the critical force and emergent frequency are qualitatively similar when viewed as functions of either λω 0 or De.However for β = 0.1, there is a substantial difference in these curves.The critical force and emergent frequency are single-valued functions of De.Also, the shapes of these curves as functions of De are qualitatively similar to the shapes predicted by the asymptotic analysis as β → 1 and corresponding curves for β = 0.5.Both the critical force and the emergent frequency are multivalued in the same range of relaxation times ( 0.1749 < λω 0 < 0.2955 for β = 0.1).For each value of λ in this range, there are three different bifurcation points each with their own frequency, and hence their own distinct Deborah number.This explains why the critical force and emergent frequency are single-valued functions of De.
In later results we use both the scaled relaxation time λω 0 and the Deborah number λω.The former is particularly useful when examining how quantities change with relaxation time.Because the emergent frequency depends on the relaxation time the Deborah number is not simply proportional to the relaxation time.

IV. FREQUENCY ANALYSIS FOR A FIXED FOLLOWER FORCE IN A VISCOELASTIC FLUID
A. Frequency changes for varying relaxation times In the previous section, we examined how the emergent frequency at the bifurcation depends on the fluid parameters.However, the bifurcation location depends on the fluid parameters.Hence as the fluid parameters vary the follower force varies as well.Here we explore how the emergent frequency depends on fluid elasticity and viscosity at a fixed follower force.Close to the bifurcation the angular frequency is approximately the imaginary part of the eigenvalue with positive real part.In order for the results from linear stability analysis to be relevant, we consider solutions that are close to the bifurcation.For a fixed β, the critical force where oscillations emerge is non-monotonic in λ, and there is a maximum critical force as a function of λ; e.g.see Fig. 4. We choose the force to be approximately 1% higher than the maximum for a particular value of β and examine how the frequency changes as a function of λ.This choice of force is large enough to avoid bifurcations in  λ where the oscillations cease and small enough to avoid large amplitude motion where the linear analysis is less accurate.For example, for β ≥ 0.5 the local maximum in force σ ≲ 1.3σ 0 which is still relatively close to the bifurcation for all relaxation times.In Fig. 6 (a) we show the emergent frequency scaled by ω 0 at a fixed follower force as a function of the scaled relaxation time λω 0 for a range of β ≥ 0.5.The emergent frequency is monotonically increasing for a fixed follower force for each β and the frequency increase levels off for large λω 0 .The emergent frequency also increases with decreasing β.For β = 0.9 the frequency at λω 0 = 4 is about 10% higher than the viscous frequency at the same force, whereas for β = 0.5 the frequency at λω 0 = 4 is nearly double the viscous frequency at the same force.
In Fig. 6 (b) and (c) we show color-fields of the emergent frequency as a function of both λω 0 as well as σ/σ 0 for β = 0.9, 0.5.The fixed follower force strength corresponding to the figure above are highlighted with grey dashed lines.For follower forces fixed at higher values (above the grey line), the qualitative behavior of the frequency is the same, namely the frequency increases rapidly for λω 0 ≲ 1 and levels off for higher relaxation times.Quantitatively, higher forces lead to higher frequencies overall.

B. Comparing analysis and numerical simulations
To explore how well the linear stability analysis predicts the emergent frequencies away from the bifurcation we solve Eqs. ( 8) -( 9) numerically with boundary conditions given by Eqs. ( 10) - (11), which accounts for both normal and tangential deformations.Details of the numerical method are described in the Supplementary Information.With these simulations we are also able to examine how the amplitude and shape change with varying fluid rheology.In Fig. 7 (a) we compare the results for emergent frequency in the linear stability analysis and the simulations for β = 0.5, σ/σ 0 = 1.3.The corresponding amplitude of the tip of the filament is shown in Fig. 7 (b).Despite the fact that the amplitude is not particularly small, the simulations qualitatively match the frequency predicted from the linear stability analysis, with the simulations exhibiting slightly higher frequencies for moderate λω 0 .The amplitude of the oscillation is not available from the linear analysis, but the amplitude is related to the distance from the bifurcation.Because the oscillation emerges at a Hopf bifurcation, the amplitude should grow like (σ − σ c ) 1/2 , where σ c represents the follower force strength at the bifurcation.The fact that the amplitude initially decreases and then increases to a constant value as λ increases is consistent with how the distance to the bifurcation changes.
In Fig. 7 (c-h) shapes (left) and normalized curvatures (right) for the simulations at λω 0 = 0, 0.75, 3.5 (corresponding to the red dots on the amplitude figure above) are shown.The plots in Fig. 7 (c,e,g) show the filaments at the same phase in the period (phase is labeled by color) and the actual frequency of motion is listed in the figure.Other than changes to amplitude and frequency the overall shapes are similar.This can be seen more clearly in the kymographs of curvature normalized by its maximum value in Fig. 7 (d,f,g).The peak values of curvature are 3.01, 1.15, and 3.56 for λω 0 = 0, 0.75, and 3.50, respectively.The amplitude of the curvature follows the same trend as the tip amplitude, which as noted above is related to the distance from the bifurcation.The curvatures exhibit subtle differences but are similar.
Although one may expect more significant shape changes as a function of rheology, the shape of the oscillating filament near the bifurcation is determined by the eigenfuctions of the operator L in (17).These eigenfunctions do not depend on the fluid elasticity (λ nor β).It is only the eigenvalues of the system that change with fluid elasticity.We only considered parameters near the bifurcation point to remain in the low amplitude regime so that linear viscoelasticity was a reasonable approximation.
Other analyses of filaments subject to follower forces in viscous fluids demonstrated significant shape variation and different kinds of motion at large amplitudes that depend on the boundary conditions and distribution of follower forces [29,30].There may be shape changes due to fluid elasticity at higher amplitudes, but at large amplitude one must consider viscoelastic nonlinearity [13,36].

C. Comparing frequency changes with experiments
In [8] the flagellar beat pattern, beat frequency, and swimming speed of the biflagellated alga C. reinhardtii were measured in response to systematic variation of the fluid viscosity and relaxation time (or fluid elasticity).Surprisingly, it was observed that the beat frequency increased with increasing relaxation time.In a Newtonian fluid the beat frequency decreased monotonically with the fluid viscosity, but in viscoelastic fluid the frequency changed nonmonotoically with the fluid viscosity -initially decreasing, then increasing and appearing to plateau.The physical origins of the frequency response to fluid elasticity are not currently understood.Here we examine the predictions of the follower model to compare with these experimental observations.
It is reasonable to ask whether one expects the predictions from the analysis of the follower model to be relevant to flagella which are driven by dynamic motor forces along the filament.If the motor activity is constant or does not change with the rheology of the fluid, then the predictions of our analysis that follow from the relationship between the eigenvlaues in eq. ( 18) will hold.Namely, the analysis in §III predicts a higher beat frequency in a viscoelatic fluid which approaches a constant as the relaxation time increases.Several studies have examined different mechanisms of motor feedback and control in the linearized equations, and these models are capable of matching experimental data [19,20,24].While there are different proposed control mechanisms, all of them are assumed to depend on the frequency of the emergent beat.In the Supplementary Information we show that even when the motor activity changes with emergent frequency, in the limit of vanishing polymer viscosity (β → 1), the expression for the frequency at the bifurcation is of the same form as (30).Thus it is expected that the frequency generally increases with relaxation time and approaches a constant value in the limit of large relaxation time.As we show below, these two features of the frequency response to fluid elasticity are consistent with the data from [8], and thus our analysis may provide insight into the mechanisms underlying the observations.The viscoelastic fluids in [8] were prepared by adding small amounts of the high molecular weight, flexible polymer polyacrylamide to water.The addition of polymer increases the fluid relaxation time, but it also changes the total viscosity.In terms of the parameters used in this work, both β and λ change simultaneously.
In order to compare the predictions of the follower model with these experiments, we fit the rheological data from [8] to find functional forms for the polymer viscosity, µ p (c), and the (dimensional) relaxation time, τ (c), as functions of the polymer concentration c in ppm.Details of our fitting procedure are given in the Supplementary Information.Using these fits and the non-dimensionalization in Sec.II, we obtain the dimensionless relaxation time and solvent fraction λ(c) = τ (c)k b /(L 4 (µ s + µ p (c))) and β(c) = µ s /(µ s + µ p (c)), respectively, as functions of the polymer concentration.
We use these models for λ(c) and β(c) in simulations for a fixed follower force σ/σ 0 = 1.5, and k b /L 4 = 0.25.Other choices for σ/σ 0 and k b /L 4 are considered in the Supplementary Information.We vary the concentration over the range c = 0 − 80 ppm and compute the dimensional frequency of the oscillation, ω V E (c).In Fig. 8(a) we plot this frequency normalized by the frequency at c = 0 as a function of the total viscosity, (µ s + µ p )(c).We compare these results with the experimental data from [8] plotted in the inset.It is remarkable that frequency changes in a viscoelastic fluid predicted by the follower model agree qualitatively with the experimental data.Specifically, the non-monotonicity of the frequency dependence on viscosity as well as the plateau for high viscosity are all captured by the follower model.
In a Newtonian fluid, the frequency of oscillations in the follower model at fixed follower force is inversely proportional to the viscosity.This is because in the dimensionless equations the strength of the follower force is the only parameter, and the time scale in the nondimensionalization is proportional to the viscosity.On Fig. 8(a) we include a plot of the normalized oscillation frequency in a Newtonian fluid: ω N (c)/ω N (0) = (µ s +µ p ) −1 , and the corresponding experimental data is shown in the inset.Both the follower model and the experimental data show that the frequency decreases with increasing viscosity.However, the follower frequency is inversely proportional to viscosity, but as discussed in [8], the measured frequency in Newtonian fluid scaled like (µ s + µ p ) −1/2 for large viscosity, which is consistent with the frequency scaling predicted by a model that includes force-sensitive dynein motor activity [19].
The non-monotonic response of the frequency to viscosity can be explained using the asymptotic analysis in Sec.III C(a).For low polymer concentrations, β is close to 1, and the asymptotic expression for the frequency in Eq. ( 30) holds.Redimensionalizing this expression, the frequency decrease from increasing viscosity occurs at first order in concentration, but the frequency increase from elasticity occurs at second order in concentration.Therefore for small concentrations, the frequency should drop at the same rate as the frequnecy in a viscous fluid, as is observed in Fig. 8(a).Further, at high polymer concentration the relaxation time is large, and the frequency is approximately given by Eq. (24).Because β grows at the same rate that the time scale decays, when Eq. ( 24) is redimensionalized it predicts that ω V E (c) ∼ ω V E (0); i.e. the frequency should approach the frequency at zero polymer concentration.In Fig. 8(a) we see the frequency appears to approach a value 10% higher than predicted, but recall that Eq. ( 24) holds at the bifurcation point, and the results in Fig. 8 are away from the bifurcation.In summary, the analysis predicts that the frequency should initially drop at low polymer concentrations when viscous effects dominate, but it must eventually increase and approach a constant in a regime in which viscous and elastic effects counter balance each other.These same trends were observed in experiments [8], and thus this analysis provides a possible mechanistic understanding of the observed frequency response to changes in viscosity in viscoelastic fluids.
Viscosity can be varied for both Newtonian and viscoelastic fluids.In order to isolate the effects of elasticity on frequency in [8] the frequency was measured for both Newtonian and viscoelastic fluids at the same total viscosity.The frequency in a viscoelastic fluid relative to the frequency of a viscous fluid of the same viscosity was reported based on Deborah number De = λω.Similarly, in Fig. 8(b) we plot the frequency of the oscillation in a viscoelastic fluid normalized by the corresponding frequency in a Newtonian fluid of the same viscosity as a function of De.We contrast these results with the results in Fig. 6 for the frequency as a function of the relaxation time for fixed β.In both cases, the frequency increases as the relaxation time increases, but when the viscosity is fixed, as in Fig. 6, the frequency levels off for high relaxation time.When the relaxation time and viscosity change together (via the polymer concentration), as in Fig. 8(b), the frequency does not approach a constant.These results are, again, qualitatively similar to the experiments from [8] (see inset).The agreement between the follower model and the data on C. reinhardtii is remarkable given that flagella are powered by dynamic internal molecular motors and the follower model we analyzed is driven by a fixed-strength external force.This suggests that the frequency response to fluid elasticity that arises from changes in the fluid drag may apply more generally to systems with different driving forces.

V. DISCUSSION
Despite the fact that fluid rheology is known to affect the shape and frequency of beating flagella in many biological systems, there is limited data in which the fluid elasticity is systematically varied, and a mechanistic understanding of how fluid elasticity affects emergent motion does not exist.Theoretical explorations have shown how fluid elasticity can change the shape of the beat with prescribed active forces [12,13], but these models cannot be used to understand the frequency changes observed in experiments [8].In this paper we extend the model of an elastic filament driven by a follower force at the tip from [28] to examine how fluid elasticity affects the emergence of oscillations and their resulting frequency.
As in a viscous fluid [28], there is a Hopf bifurcation at a critical force at which beating emerges.Our analysis identified how the bifurcation location depends on the relaxation time and viscosity ratio.In a viscoelastic fluid the force required to induce oscillations is always higher than in a viscous fluid.Moreover, there are parameter regions where increasing fluid relaxation time will stabilize an oscillating filament, but upon further increase in the relaxation time the filament will again oscillate at a higher frequency.Our analysis predicts that fluid elasticity generically increases the frequency of beating over the same filament in a viscous fluid at the same force in agreement with the experimental observations in [8].When the relaxation time and total viscosity increase in tandem through the polymer concentration, competing effects of elasticity and viscosity lead to a non-monotonic response of frequency on viscosity that again agrees well with experiments [8].
In [8] it was observed that although the frequency of the beat increased in viscoelastic fluids the swimming speed decreased.The shape of the beat changed significantly in viscoelastic fluids, namely the maximum of the flagellum curvature increased, but the bending at the basal end was reduced.In [37] we performed numerical simulations of swimmers based on the gaits from [8] to separate the effects of changes in gait and changes in fluid rheology on the swimming speed.This work showed that the reduction in speed resulted from both the change of the shape of the beat and the nonlinear growth of elastic stress around the flagella.The model analyzed in this work does not predict shape changes in response to viscoelasticity.At low amplitude the shape is determined by the eignefunctions of the linearized operator, which do not depend on the fluid elasticity.Although the frequency response predicted from our analysis is consistent with [8], the reported shape changes in the flagella beat in response to fluid elasticity cannot be captured with the model analyzed here.Capturing shape changes due to fluid elasticity requires a more sophisticated model of the active forces from molecular motors.
The model we analyzed did not include mechanical feedback on the driving force.There are many theories about how mechanical feedback on molecular motors leads to spatiotemporal coordination of motor activity to produce the flagellum beat, but it has also been shown that the motor coordination is not necessary for producing the beat [25].Even if the mechanical feedback on motor activity is not responsible for coordination, motor regulation could play a role in modulating the flagellum beat.Our analysis captures the qualitative changes in the frequency observed in [8] in response to fluid elasticity, but quantitative agreement between the model and data would likely require a more sophisticated model that includes dynamic motor activity.For example, our analysis predicts that in a Newtonian fluid, the frequency is inversely proportional to the viscosity (see also [25]), but as discussed in [8], the motor model from [19] predicts that the frequency is inversely proportional to the square root of the viscosity, which was a better fit to the data (see Fig. 8).
One approach to analyzing how the shape and frequency of the flagellum beat is shaped by motor models that incorporate mechanical feedback is to examine time-periodic solutions of linearized equations [19,20,23,24].This approach is equivalent to examining the solutions at the bifurcation point as we have done here.The form of linearized equations, and thus the eigenvalues of the corresponding operator, depends on the motor model, but the eigenvalues in a viscous fluid are related to those in a viscoelastic fluid by (18).This suggests that the effect of fluid elasticity on frequency discussed here may be a generic effect in models of flagella which incorporate feedback or regulation from molecular motors.

Supplementary Info C: Varying Stiffness
Comparing the model predictions as functions of the dimensional relaxation time and viscosity requires choosing values for the parameter combinations k b /L 4 and σ/σ 0 .We used k b /L 4 = 0.25 and σ/σ 0 = 1.5 in the simulations presented in Fig. 8.In order to explore how these choices affect our results, we consider a range of k b /L 4 from 0.1 − 2 and σ/σ 0 from 1.25 − 2.5.To explore parameters efficiently, we use the frequencies predicted from linear stability analysis rather than from the numerical simulations.These two different methods of computing the frequency agree well near the bifurcation but start to diverge away from the bifurcation.
Before we present the results for the wider range of parameters we show a comparison of the results from the simulations and the linear stability analysis for the values presented in the manuscript.Figure S2 compares the linear stability analysis and simulation computation of the frequency for the parameters k b /L 4 = 0.25 and σ/σ 0 = 1.5 (presented in Fig. 8).In Fig. S2 (a) the simulations predict larger frequency at high polymer viscosity than the analysis.Nevertheless the qualitative feature of the non-monotonic dependence of frequency on viscosity in viscoelastic fluids is still predicted by the analysis.Fig. S2 (b) shows that the simulation and analysis show close agreement for the frequency boost seen in viscoelastic fluids over viscous fluids as the Deborah number is increased.In Figs.S3 and S4 results of the exploration of the parameters for the ranges 0.1 ≤ k b /L 4 ≤ 2 and 1.25 ≤ σ/σ 0 ≤ 2.5 are presented using the linear stability analysis calculation for frequency.In Fig. S3 we plot the normalized frequency ω(c)/ω(0) as a function of the total viscosity.
Each panel in the figure shows the results for a single value of k b /L 4 for increasing forcing strength from 1.25 to 2.5 corresponding to the colors going from light to dark.All but the smallest value of k b /L 4 , σ/σ 0 predict the non-monotonic response of frequency as a function of total viscosity.The exact parameter values will change the amount that the frequency is predicted to decrease as well as where that predicted decrease is maximized, but a wide range of parameter values show the qualitative behavior reported in Fig. 8 (a).
In     where σ = σ − σ 0 and ω = ω − ω 0 .We now look for the viscoealstic bifurcation point when 1 − β = ϵ is small.As in the follower force problem, the viscous and viscoelastic eigenvalues are related by Notice that z contains all the dependence of the motor activity on the frequency.When z = 0, we recover the result in the paper for the follower force.This expression for the frequency is of the same form as that obtained for the filament driven by a follower force given in equation (30).Thus this analysis predicts that if z < 1, then in the limit of small polymer viscosity a filament driven with internal motors consistent with (S2) will exhibit a frequency that initially increases with relaxation time and approaches and a constant in the limit of large relaxation time.These two features of the frequency response are consistent with follower model and the experimental data from [8].

FIG. 2 .
FIG. 2. Real part (a) and Imaginary part (b) of the eigenvalue with largest real part for a viscous fluid as a function of the follower force, σ.(c) Real part of the four eigenvalues with largest real parts.Blue dashed lines denote real eigenvalues, solid black lines denote pairs of complex eigenvalues, and red dots mark where the real part changes sign.

FIG. 3 .
FIG.3.Plots of the asymptotic and numerical solutions as β → 1 for the (a) critical follower force and (b) emergent frequency at the bifurcation as functions of the relaxation time for β = 0.95.The two red dots mark the bifurcation points and emergent frequency, respectively, for σ = 1.01σ0.

1 FIG. 4 .
FIG.4.Location of the bifurcation point in the σ-λω0 plane for β = 0.1, 0.2, . . .0.9.Inset: Location of the bifurcation for β = 0.1 zoomed in to show that between the two dashed vertical lines there are three bifurcation points as the follower force changes.

FIG. 5 .
FIG. 5. Bifurcation location (a,c) and frequency (b,d) at the bifurcation for β = 0.5 (a,b) and (c,d) β = 0.1.The two curves in each panel represent two different scalings for the relaxation time, λ.The red curves (bottom axes) show how the data depend on λω0, where ω0 is the angular frequency at the bifurcation point in a viscous fluid.The blue curves (top axes) show the data depend on the Deborah number De = λω, where ω is the emergent frequency in a viscoelastic fluid.Because the emergent frequency depends on the relaxation time, De represents a nonuniform scaling of the relaxation time.

FIG. 6 .
FIG.6.Emergent frequency (a) as a function of scaled relaxation time for β ≥ 0.5 at a fixed follower force that is 1% higher than the maximum force at the bifurcation for the corresponding β.Color fields of frequency for β = 0.9 (b) and β = 0.5 (c).Solid lines show location of bifurcation, and dashed line shows fixed force value for corresponding frequency values shown above.

FIG. 7 .
FIG. 7. Emergent frequency (a) and tip amplitude (b) for β = 0.5, σ/σ0 = 1.3.Frequency figure shows a comparison with simulation and linear stability analysis.Emergent shapes with 20 snapshots per period (c,e,g) and curvature normalized by its maximum (κ/κmax) (d,f,h) over a period for simulations with λω0 = 0, 0.75, 3.5, these values are highlighted with red dots above.The maximum values of curvature for these three values of λω0 are 3.01, 1.15, and 3.56, respectively.The amplitude of the curvature follows the same trend as the tip amplitude.

FIG. 8 .
FIG.8.Frequency scaled by its value at c = 0 as a function of the total viscosity (a).Newtonian viscosity follows 1/viscosity scaling for comparison.Frequency relative to viscous frequency as a function of the Deborah number (b).Simulations run using fixed value of force σ/σ0 = 1.5, and k b /L 4 = 0.25.Inset graphs using data from[8] FIG. S2.(a) Simulations and LSA for frequency vs. total viscosity at parameters in main paper.(b) Simulations and LSA for normalized frequency vs. De at parameters in main paper.
Fig. S4 the frequency ω V E /ω N is plotted as a function of De.As before each panel shows the reults different value of k b /L 4 for a range of forcing strengths from 1.25 to 2.5 corresponding to the colors going from light to dark.The difference between these graphs is what values of De are sampled for the given mechanical parameters.They all show the same qualitative behavior reported in Fig. 8 (b).
FIG. S5.The product of the real part of the viscoelastic eigenvalue with smallest real part and the relaxation time for (a) β = 0.9 and (c) β = 0.1.The black line denotes the location of the bifurcation.This same quantity as a function of λω0 for the fixed value of σ = 2σ0 for (b) β = 0.9 and (d) β = 0.1.These plots help illustrate that ηveλ → −1/β as λ → 0 and ηveλ → −1 as λ → ∞.The red dotted line on (d) represents values in the stable region where both eigenvalues have negative real parts.
FIG. S6.Plots of the asymptotic and numerical solutions for the critical follower force (left column) and emergent frequency (right column) at the bifurcation as functions of the relaxation time for different values of β.The asymptotic expressions for the critical force and frequency are given in equations (29) and(30), respectively.