Ambient Fluid Rheology Modulates Oscillatory Instabilities in Filament-Motor Systems

Semi-flexible filaments interacting with molecular motors and immersed in rheologically complex and viscoelastic media constitute a common motif in biology. Synthetic mimics of filament-motor systems also feature active or field-activated filaments. A feature common to these active assemblies is the spontaneous emergence of stable oscillations as a collective dynamic response. In nature, the frequency of these emergent oscillations is seen to depend strongly on the viscoelastic characteristics of the ambient medium. Motivated by these observations, we study the instabilities and dynamics of a minimal filament-motor system immersed in model viscoelastic fluids. Using a combination of linear stability analysis and full non-linear numerical solutions, we identify steady states, test the linear stability of these states, derive analytical stability boundaries, and investigate emergent oscillatory solutions. We show that the interplay between motor activity, filament and motor elasticity, and fluid viscoelasticity allows for stable oscillations or limit cycles to bifurcate from steady states. When the ambient fluid is Newtonian, frequencies are controlled by motor kinetics at low viscosities, but decay monotonically with viscosity at high viscosities. In viscoelastic fluids that have the same viscosity as the Newtonian fluid, but additionally allow for elastic energy storage, emergent limit cycles are associated with higher frequencies. The increase in frequency depends on the competition between fluid relaxation time-scales and time-scales associated with motor binding and unbinding. Our results suggest that both the stability and oscillatory properties of active systems may be controlled by tailoring the rheological properties and relaxation times of ambient fluidic environments.

A striking and well-studied feature of these filament-motor systems is the emergence of stable oscillations under favorable conditions. Examples include oscillatory instabilities in muscles and sarcomeres [19][20][21] and in microtubule assays involving NCD motors [22,23], the periodic undulations of eukaryotic cilia [7,8], and spontaneous oscillations during asymmetric cell division that involve microtubules interacting with cortical force generators [24][25][26]. These temporally varying patterns may follow complicated periodic functions, and are furthermore very tunable. Often, the onset of these oscillations may be controlled by quenching the system by depleting the molecular motors (the active agents driving these oscillations) of their energy source (ATP or Ca 2+ ). All these systems are open systems with continuous energy input, and highly non-linear. In all, these observations suggest that the onset of these oscillations may be interpreted as instabilities to a base state triggered by system changes.
Biological filaments and motors usually inhabit fluidic or gellike media with non-Newtonian properties and viscoelastic properties. Such complex fluid environments can directly impact these active systems by exerting time-dependent stresses. For instance, a filament moving in a Newtonian fluid feels viscous drag stresses exerted by the ambient medium. A filament moving in viscoelastic environments will be subject to more complicated time dependent forces that include dissipative (viscous) and non-dissipative (elastic) components. Recent work on cilia in the green alga Chlamydomonas reinhardtii that enables whole-cell locomotion suggests a strong influence of the viscosity and complex rheology of the ambient fluidic medium [9,27]. Sketched in Figure 1A is a schematic summarizing experiments from [9] that depicts the variation in ciliary beat frequency of an algal cell in Newtonian (red) and polymeric (blue) fluids. Frequencies are seen to be significantly higher in the viscoelastic polymeric fluid compared to a Newtonian fluid with the same viscosity. Similarly strong effects of viscoelasticity on ciliary beat frequency are seen in the human mucociliary tract. In diseased states or when there is ciliary dysfunction-as seen for instance in cystic fibrosis (CF), primary ciliary dyskinesia (PCD) or chronic obstructive pulmonary disease (COPD)-the mucus layer submerges cilia in an extremely viscoelastic environment that impacts their beating [10,[12][13][14].
Thus the effect of the rheology of the ambient fluid needs to be understood in order to predict the dynamical response in these biological filament-motor systems. There is however a significant gap in our understanding of these aspects since it is experimentally difficult to independently vary fluid viscosity and elasticity. Abstracting the ingredients crucial to the examples of actively driven motor-filament systems allows us to identify three main ingredients that may control system-level mesoscale and macroscale dynamics-1) passive elasticity, from intrinsic filament-motor elasticity, 2) rheological properties of the embedding medium, and 3) associated motor mechanochemistry and kinetics. Here, we will combine these elements to build a simple phenomenological model for a filament-motor aggregate interacting with model fluidic media. To enable analytical progress, we assume the ambient fluid to correspond to one of three canonical types-a Newtonian fluid, a Kelvin-Voigt gels or a Maxwell medium. Using our minimal model we investigate the onset, stability and control of oscillations and address associated questions: 1) how do variations in the fluid induced mechanical stresses affect the manner in which the internal biomechanical state of the molecular motors in the aggregate is modulated, 2) how do these modulations impact the onset of oscillations, and when may they stabilize an otherwise unstable system, and 3) how does fluid rheology affect the frequency of emergent oscillations and their amplitude? Biological filament-motor systems are of course much more complex than the minimal model analyzed here. Nonetheless, minimal models offer an essential tool for understanding basic biophysical mechanisms. In this sense, the simple minimal model analyzed here using model fluids represents an important step towards more complete computational models.

MINIMAL MODEL
Our model system is shown in Figure 1A (bottom pane). We model the filament bundle-referred to henceforth as (composite) filament-as a rigid sheet that is at a certain fixed height from an underlying flat substrate. In the reduced setting analyzed here, the sheet has length ℓ A , thickness b A ≪ ℓ A and lateral extent w A satisfying ℓ A ≫ w A ≫ b A . The view shown in the figure is from the side, the axial dimension is the length and the lateral dimension into the paper (corresponding to the lateral extent, b) is treated as a neutral direction. The area of the segment that faces the underlying substrate is ℓ A w A . Biofilaments and aggregates may have slender cylindrical geometries; the area ℓ A w A is thus to be interpreted as the projected area of the aggregate that faces the substrate.
The filament is held a fixed distance away from an underlying flat substrate by passive, permanently attached linear springs (blue springs) with areal density ρ p , stiffness K p , and equilibrium length ℓ p . These passive springs resist shearing deformations and lateral translation of the filament. Additionally, the filament also interacts with a collection of model molecular motor proteins that are shown in purple in the figure. These motors are grafted to the lower plate with an areal density ρ m . We treat these molecular motors as active linear elastic springs with spring constant k m and equilibrium rest length ℓ m . The tail of each motor is permanently attached to the substrate and thus we do not allow for diffusion of the motors. Motor heads meanwhile can periodically adhere to the filament and when attached generate forces displacing the filament laterally. We choose the inter-motor spacing to be small enough, and densities high enough that a continuum description of filament-motor interactions suffices. Finally, for ease of analysis without loss of generality, we set ℓ m = ℓ p = 0.
To aid in the analysis, we model motor-filament interactions using two-state cross-bridge models [8,19,21]. Each motor is either attached state or detached. Detached motors can attach to the filament in a forward-leaning position with specified attachment probabilities. They then quickly undergo a conformational change, which makes them strained. As the filament moves, so does the motor and this results in the extension of the motor and thus extension of the internal spring (with spring constant k m ). Attached motors can detach any time, and the statistics of this process is embodied via microscopic strain dependent detachment probabilities. The transitions between attached and detached statesthe kinetics of the mechanochemical cycleare thus determined by motor kinetics, and specifically the attachment and detachment probabilities (rates). Detailed microscale equations modeling the mechanochemistry and motor-filament interactions in systems such as these have been derived and used by others and by us (c.f [19,28] and references therein). Briefly summarized, motor kinetics may be described by a set of population balances relating the attached and the detached probability densities to the attachment and detachment fluxes via microscopic transition rates. For simplicity we neglect motor diffusion, and consider the noiseless mean-field limit, where the number of motors (and the density) is large (high) enough that fluctuations in the time average of the density of attached motors are small compared to the mean value. When the distribution in extension of attached motors is sharply peaked about the typical (average) length, transients to this distribution occur over times very small compared to the averaged macroscopic time scale. The extension of attached motors is then peaked about the mean value with small deviations from the mean. We assume that detached motors relax to a delta function with the change occurring instantaneously. Coarse-graining the microscale dynamical equations provides a mean-field continuum description for the (mean) attached motor fraction, and the mean attached motor extension.

GOVERNING EQUATIONS FOR THE MINIMAL MODEL
Following the physical picture described in §2, we next derive equations governing the dynamics for our minimal system -an elastically constrained filament-motor assembly embedded in a viscoelastic fluid. To enable analytical progress, we assume that the rheology of the embedding fluidic medium behaves in one of three ways: as 1) a Newtonian fluid, as 2) a Kelvin-Voigt material, or as a 3) Maxwell medium. Based on these, filament-motor and filament-medium interactions may be modeled analytically.
Consider again the schematic in Figure 1A. As mentioned earlier, the (composite) filament interacts with motors via a projected area ℓ A w A where w A is the width in the neutral direction. To model the effect of external elastic constraints on the filament, we choose to anchor the filament to a boundary with a linear elastic anchoring spring with stiffness K eff and rest length L o . We ignore edge effects and choose a reference frame with origin at the left attachment point where the spring is anchored. Due to the action of the motors, the rigid filament moves horizontally; with its motion resisted by the anchoring spring and the passive substrate attached springs. Filament motion is also resisted by fluid drag forces arising from the ambient viscoelastic medium. As a result of these restoring forces, when forced by attached molecular motors, the equilibrium stable state of the filament can either be one of rest or a time dependent state.
We define the displacement of material points on the plate relative to the rest state (where no motors are attached to the filament and the spring is not strained) by Up = L − L o ; the speed of the plate is _ U p . In the absence of fluid and solid inertia, the sum of forces on the filament must vanish. If F p m is the effective active force density per attached motor, force balance on the filament provides the equation The active force F p m in Eq. 1 is obtained from a knowledge of two variables describing the state and quantity of activity-the fraction of attached motors N and a mean motor extension Yp. We follow a model that has been studied by others and by us previously [19,28] that incorporates an experimentally motivated mean-field description of motor attachment and detachment kinetics. Motor kinetics is captured by the meanfield attachment rate ω on , and the mean-field detachment rate ω off . The overall motor density of motors that interacts with the filament is ρ m , a fraction N of which is attached at any time.
In the context of this simplified picture, the evolution of the fraction of attached motors, N follows Here the mean attachment rate ω o on is constant while the mean detachment rate ω off ω o off F incorporates the state dependent function F (Y p ) that depends on motor extension (stretch). This quantity-the stretch of attached motors, Y p -follows Given N and Y p , the effective active force per motor in Eq. 1 is given by We next model the fluid interaction term, F p d , in Eq. 1 by specifying constitutive laws for the manner in which the fluid interacts with the filament.

Newtonian fluid:
In the simplest case where the medium is a Newtonian fluid, the fluid drag per unit length acting on the plate depends just on the local shear rate at the plate surface. Since our intention is to capture qualitative aspects we proceed with scaling laws. Let us assume that the shear rate induced in the fluid at the filament surface is uniform along the filament and the velocity field penetrates to a distance ℓ N so that an aerage measure of the local strain rate is (1/ℓ N )dU p /dt p . An approximate equation for the drag force is then with μ N being the Newtonian viscosity. The reduction is consistent with an effective drag coefficient approximation used for motion of bodies at low Reynolds number (with drag force proportional to the instantaneous speed). Note that the segment length ℓ A , segment width w A and length ℓ N are all assumed to be constants. This is a major simplification since this penetration depth can be time dependent.
2. Kelvin-Voigt fluid: When the embedding medium is a networked gel-like medium with both elastic and viscous features, the drag forces depend both on plate velocities as well as displacements. One may consider these displacement contributions to arise as a result of restoring forces as the plate moves relative to a fixed meshed elastic network. We choose a Kelvin-Voigt like model with a single retardation time scale and equilibrium mesh size ℓ KV represented by relating the strain to the stress. Evaluating this at the surface of the plate and assuming no slip, we obtain (measuring strain in mesh lengths, and assuming a continuum formulation holds) 3. MaxwelL gel: For a viscoelastic Maxwell-like with a single relaxation time scale and a mesh length ℓ M , we have the relationship Evaluating this at the surface of the plate and assuming no slip, we obtain (again measuring strain in units of mesh length and assuming a continuum formulation holds so that Equations 1,5,6,7 assume that the ambient medium is always quasi-statically coupled to motion of the upper plate. A similar analysis may be done for the Zener standard solid model of the type illustrated in Figures 1B-iv, but we do not pursue this in the current work.
We emphasize that Equations 5-7 ignore transient stress fields induced as the viscous/viscoelastic boundary layer forms near the moving filament and associated time-dependent features. A more accurate description would necessitate solving the full Navier-Stokes equation with the appropriate constitutive equation relating the stress in the fluid to the strain and strain rate and is beyond the aim of this paper.
To obtain scaled equations that allow the identification of important non-dimensional parameters, we scale time with 1/ω o on , and motor/filament extensions with A. Defining t p T(ω o on ) −1 , Y p ≡ AY and U p ≡ AU allows us to recast (1), in dimensionless form. In this process, we find that re-scaling anchoring spring stiffness by defining K eff ≡ K eff /ℓ A w A provides a further reduction. Since the fraction of attached motors N ≡ ρ a / ρ m where ρ a is the density of attached motors, the density of detached motors is ρ m − ρ a . Tables 1 and 2 summarize the scaled parameters in our model. Equations 8-10 constitute the final scaled nonlinear ODE's governing the evolution of the motorfilament system embedded in the viscoelastic medium.
with dimensionless parameter β quantifying motor density while dimensionless parameters B 1 to B 4 quantifying the viscoelasticity of the fluid and fluid specific relaxation times. We would like to point out here that Equation 10 is applicable to Newtonian, Maxwell, Kelvin-Voigt and the standard solid models (illustrated in Figures 1B,I-1iii.
The anti-Zener Jeffrey element depicted in Figures 1B-1iv involves a more complicated relation with a second derivative of U as follows In this article we will consider not consider models resulting in a second derivative of U. The linear stability analysis can however be extended to the Jeffrey element, albeit with slightly more involved algebra.
Following previous work, we model motor kinetics as responding to the deformation experienced by attached motors [1,2,19]. Attached motors behaving as slip bonds with a strain dependent detachment probability. This functionality is encoded in the detachment rate ω off (Y) ω 0 off F (Y) that modifies the bare detachment rate by a multiplicative factor F (Y). For the analytical investigations, we restrict F (Y) to functions that is at least twice differentiable but otherwise general. For plotting, we will use the following forms Equations 8-11 admit a stationary state given by the set (N, In Equation 11, the dimensionless parameter E s quantifies the characteristic force attached motors can withstand before detachment. Note that the ratio B 1 /B 2 compares the motor time scale with the relaxation time scale of the Kelvin-Voigt medium. Dimensionless parameter B 3 is similarly a ratio of motor time-scales to fluid relaxation times scales for the Maxwell fluid. Physical interpretations of variables (U, N, Y) and that of dimensionless parameters are summarized in Tables 1 and 2.
In our model, bending is ignored and the filament composite/ aggregate is considered rigid and inextensible. While our model is not directly applicable to ciliary oscillations, using some numerical estimates relevant to filaments and motors in cilia provides understanding of the magnitudes of the parameters in the model. Our model also may allow us to interpret how oscillations may arise within and along an initially straight passive cilium. Localized elastically weak regions may result in cilia fragments driven by dynein arrays moving against their neighbors. In this scenario, we estimate k m~1 0 −3 N/m, the effective linear density ρ m w A~O (10 8 ) m −1 , w A~4 0−60 nm, k N~1 6-100 pN μm −1 and w A ρ N~1 0 5 -10 7 m −1 [1,29]. The value of K eff depends on the material. Microtubules are relatively stiff with large persistent lengths and the Young's modulus E1 .2 GPa providing the stiffness per length K eff~G w A . For microtubules driven by dynein, previous studies report a

Symbol Definition Interpretation
A -Typical dimensionless motor displacement and strain E s k m A/k B T Strain dependent force scale for motor detachment β ρ m k m /(ρ p k p + K eff ) Scaled motor density and a measure of active elasticity Ratio of motor kinetic constants and related to motor duty ratio TABLE 2 | Dimensionless parameters quantifying the ambient medium stress on moving filament. Here, ℓ A and w A are the length and width of the filament aggregate interacting with motors and allow one to go from the stress to a force description. Displacements are scaled with A. We assume a continuum description is valid and that strains in the fluid/medium network may be appropriately defined using elastic and viscous moduli. We consider three cases:

Symbol Newtonian
Kelvin-Voigt Maxwell Frontiers in Physics | www.frontiersin.org June 2022 | Volume 10 | Article 895536 5 persistence length~5mm at room temperature, and for actin driven by myosin,~16 μm at room temperature. These length scales provide an upper bound for ℓ A . For polymeric fluids, relaxation times can vary significantly depending on the type of polymer and concentration, Carboxymethylcellulose (CMC) solutions for instance, with MW = 7 × 10 5 can have relaxation times~20−50 ms as concentration is varied from 225 to 500 ppm [30]. At very high concentrations (4000 ppm), the relaxation time is as large as 0.2 s. Shear viscosity meanwhile for these can vary from~1 mPa.s to~18.8 mPa.s for concentrations ranging from 0-500 ppm. Polymeric solutions made using polyacrylamide also allows us to access a range of relaxation times. Specifically, in [9] viscoelastic fluids prepared by adding small amounts of high molecular weight polyacrylamide at concentration in solution ranging from 5 to 80 ppm resulted in fluid relaxation times that ranged from 6 ms to 0.12 s, respectively.

Dynamics of the Unconstrained System
Before analyzing the dynamics of the constrained filament-motor system, we briefly address the dynamics of an unconstrained system. Accordingly, we set ρ p = 0, K eff = 0 and study the dynamics of the animated filaments on a rigid surface on which motors are grafted, as in motility assays. Short fragments in these assays are typically rigid enough that they move without bending. Long fragments on the other hand are susceptible to buckling instabilities as previous studied experimentally and analytically [18,[31][32][33][34][35][36] in these and similar systems involving driven active filaments.
Here, we invert the problem and study the dynamics of the motor variables (Y, N) as functions of the speed Z = dU/dT. In other words, we restrict ourselves to time scales where the dynamics of the motor kinetics is slaved to the speed of the filament. The speed Z may be envisioned as an independently controlled variable; we then enquire how steady values of N and Y depend on Z and if these solutions are linearly stable. Eq. 10 is irrelevant in this limit, and Equations 8, 9 become.
When Z = 0, Equations 13, 14 admit steady state solutions (N s , Y s ) given by For each value of Z ≠ 0, other steady solutions different from Eq. 15 exist. To study both the steady solutions and their stability, we use an in-house numerical continuation implemented in Matlab @ . Treating Z as the free variable, we plot steady solutions to Equations 13, 14 and investigate the stability of the solutions for various values of Ψ and E s and for different functional forms of F the detachment function. Figure 2 shows an example case where we keep E s and Ψ the same while using two forms of F . Solid curves (maroon and blue) represent stable solutions while dashed curves are unstable solution branches. Turning points, where these branches meet, are shown as circles.
Numerical investigation reveals that the existence and extent of the unstable region depends on the form of the detachment curve. When the detachment function is constant or a decreasing function of extension Y, no unstable regions exist. In these cases, typically each value of Z is associated with a unique value of N and Y. However when the detachment function is an increasing function of Y, some parts of the steady state solutions becomes unstable as seen in Figure 2. This is manifest clearly in bistability, evident for instance in the curves of N vs. Z with three possible values of Z (two stable, one unstable) corresponding to a single value of Z. A value of N = 0.4 in Figure 2A-(ii) and 2(b)-(ii) for instance is unstable and is susceptible to disturbances that can push the state to either the upper branch or the lower branch each with the same value of Z. The origin of this instability is the emergence of effective negative spring constants due to motor activity; such features have been studied previously [19].

Analytical Results: Oscillatory Instabilities and Emergent Frequencies
The question that arises naturally next is how K eff > 0 impacts this response. To answer this question, we turn next to the stationary base state of the full Equations 8-11 and the linear stability of the base state given by Eq. 12. We find that the base state exhibits linear instability for certain regions in parameter space. We derive analytical expressions for the locus of these critical points-the neutral stability curves. Classical linear stability analysis shows at these critical points, stable oscillatory solutions emerge via the Hopf-Andronov-Poincare bifurcation [37]. Limit cycles emanate at these points with fixed frequencies and small amplitude. We calculate the frequencies of these emergent solutions analytically and investigate their dependence on dimensionless parameters defined earlier. Predictions of our linear stability analysis are confirmed by full non-linear solutions of the equations.
To identify the stability of the stationary (base) state (Eq. 12), we analyze how small imposed perturbations to this state evolve in time. We introduce a small (amplitude) parameter ϵ ≪ 1 and write Substituting (Eq. 16) in Eqs 8-11, expanding all terms, and retaining only terms that are O(ϵ), we get the three coupled linear ODE's: in Eqs 17-19, we obtain coupled algebraic equations for (N,Ŷ) that determine the growth rate σ. Imposing solvability constraints, we get a cubic algebraic equation for σ: where The growth rates σ of Eq. 21 depend crucially on the signs of b and c. We note that a and d are always positive. Since coefficients a − d are real, Equation 21 has-from the fundamental theorem of algebra-three roots. One of these is purely real and the other two may be complex conjugate. Since we set B 4 1, and are more interested in varying the rheological parameters of the fluid, the dependence of b on B 3 is the more important dependence. Here we will restrict our analysis to the physically relevant case where b > 0 and c > 0 (see Supplementary Appendix A).
Restricting ourselves to the subspace of solutions that comprise 1 real root and a complex conjugate pair, we write (σ 1 ), σ 2 = σ + iω and σ 3 = σ − iω. When the base state is stable we have σ 1 < 0 and σ < 0. We seek incipient oscillatory instabilities such that emergent solutions may be stable and oscillating. At the onset of instability-that is on the neutral stability curve-this implies the condition σ = 0. The imaginary part ω may then be identified as the frequency of the emergent solutions.

Results From Stability Analysis
For constant Ψ and E s , Equation 22 provides the critical value of β above which one can expect oscillatory limit cycles to bifurcate from the stationary steady state. Note that an essential condition here is F 0 ′ > 0 and that the magnitude be large enough for the equation to be satisfied. Formally, we require F to be a monotonically increasing function of Y consistent with previous work (see [2] for instance) in similar systems. The periodic solutions at onset have frequency 2. Newtonian fluid: For the simplest drag-producing medium -a Newtonian fluid -we have B 1 B 3 0, B 2 > 0 and B 4 1 with specific forms listed in Table 1. The force balance reduces to with first two left hand side terms being resisting (passive) forces and the third term the driving (active) force; the right hand is zero here due to the absence of inertia as explained earlier.The frequency of the bifurcating limit cycles at onset is given by Combining (Eqs 23) and (24) yields Thus for very high viscosity μ N , the frequency scales as 1/ μ N √ . The equation for the neutral stability curve simplifies to From the analytical expression for the neutral stability curve, we deduce that for oscillatory instabilities to exist we require that F 0 ′ > 0 and be a suitably large number. For large B 2 such that B 2 ≫ (βN 0 + 1), the critical value of β for onset of instability is larger than for the case of no drag and is consistent with expectation, since extra viscous drag implies higher dissipation rates. At constant motor properties-that is, at constant values of Ψ and E s and with the form of F (Y) held fixed-the only way to enable this is by having a larger number (density) of active motors.
3. Kelvin-Voigt medium: Here, the drag force depends on the displacement and the rate of displacement. Thus we have B 4 1, B 1 > 0, B 2 > 0 while B 3 0. In this case again, oscillatory instabilities are possible at critical points provided F 0 ′ > 0. The frequency of the bifurcating limit cycles at onset is obtained as For fixed motor parameters and fixed value of B 2 , the frequency in a KV medium is thus higher than the purely viscous medium (with the same viscosity/dissipative components, B 2 ) consistent with experiment. The increment comes from the elastic component of rheology of the fluid but is augmented by the motor activity parameter Ψ and the strength of the detachment function but not by its slope. To compare with the two limiting cases analyzed earlier, we note that.
where we have defined the shift in frequency relative to a purely Newtonian fluid, Δω KV ≡ (ω KV − ω N ). We can also derive the expression for the neutral stability curve 4. MaxwelL medium: For the Maxwell medium, the displacement rate is linearly coupled to the drag force and the rate of change of the drag force. We have for this case, B 1 0, B 2 > 0, B 3 > 0 and B 4 1.
From an estimation of constants a − d, we rule out the possibility of a fold-Hopf bifurcation, that is a fold bifurcation coincident with the Andronov-Hopf-Poincare bifurcation. However we can still satisfy conditions for the existence of an oscillatory instability. Further analysis (see Supplementary Appendix A) indicates a mathematical possibility for nonoscillatory instability for very high values of B 3 for which b turns negative. We take this as evidence that the model is ill-posed for these high values of B 3 and possibly inapplicable. We do not consider these unphysical solutions, and instead focus on the physically relevant oscillatory instabilities.
We therefore look for a subset of critical points for which exactly 1 real negative root (σ 1 ) and a complex conjugate pair (±iω) with zero real part (σ = 0) exists for the cubic. The neutral stability curve is defined by the locus of points that simultaneously satisfy ad − bc = 0. The frequency at onset is given by Since B 3 > 0, the overall sign of the term multiply B 3 depends on the magnitude of the term ΨF 0 ′N 0 β that has a negative sign ahead of it. We have deduced earlier that for oscillations to exist we require that the detachment function be an increasing function of Y evaluated at the steady point Y 0 . Combining these observations, we deduce that the frequency in a Maxwell medium can be higher than the Newtonian medium provided the term in red is negative. In fact the magnitude of this term (when negative in sign) can be readily modified two ways. First, for fixed B 3 , the magnitude can be adjusted by solely changing the form of the detachment function (through ω 0 off and/or E s and independent of the rheology parameters and the motor attachment rate. Second, for fixed motor mechanochemistry as long as the term in the parenthesis is negative, increasing B 3 will increase the frequency of oscillations. Before concluding this section, we address the biophysical origin of the oscillations. The simplest scenarios-the case of no fluid, and that of a Newtonian fluid-provide hints as to the destabilizing mechanism. Linearizing (Eqs 8)-(11) about the base state (Eq. 12) provides (in frequency space) a relationship between the displacement U and an effective compliance that includes contributions from the active motor-based interactions with the filament as well as the fluid response. When there is no drag/resistive force at all from the medium, we find that the instability arises when the active compliance has a negative effective elastic coefficient i. e, active motor induced elasticity overcomes the passive elastic resistances. This provides a mechanism by which system-level oscillations are initiated. Meanwhile the dissipative components of the system-i. e, viscous effects-limit the amplitude of the oscillations. Thus stable limit cycles can emerge. Interestingly no fluid or external medium is needed to initiate and sustain oscillations in this model system.
We have three independent parameters here that determine if oscillations exists-β, Ψ and E s . In general, we need the motors to detach so that the filament can slide back and forth. Very high motor attachment rates with no detachment also impedes oscillations by increasing the effective shear stiffness of the filament. Thus there is a lower limit for oscillations in terms of Ψ. For very high Ψ, the linear stable phase-space occupies larger extent. Similarly, there is a minimum value of E s below which oscillations cannot be sustained-this is because motor detachment is not sufficient to sustain reversal of motion of the filament. Very high values of E s on the other hand imply that motors detach almost as soon as they attach. Finally, there is also a minimum value of β, or equivalently a minimum density of active motors, needed to initiate oscillations.

Amplitudes and Frequencies of Limit Cycles Vary Away From Critical Point
The linear stability analyses in Section 4.1, Section 4.2 provide information on the stability boundaries for the stationary steady state (Eq. 12) and also on the emergent frequencies of the bifurcating limit cycles at critical points on these boundaries, the neutral stability curves. The frequency when parameters correspond to positions in phase space far away from the neutral stability curves will of course differ from the values at criticality. When all parameters are kept fixed, and just one parameter-say for example, β varies, the amplitude of the limit cycles is expected to scale~ β − β c where β c is the critical value at which the instability manifests. While this scaling certainly holds close to criticality, actual amplitudes may deviate from this far from the critical point. Physically, oscillations and amplitudes are limited by a power balance. The active energy the motors generate and input into the system has to be dissipated viscously by fluid drag. Elasticity-from the anchoring spring, the underlying passive springs and from the fluid mediates oscillation amplitudes as well as provides a mechanism by which energy can be stored temporarily. Thus both the amplitude and the frequency of oscillations may vary significantly from linear stability predictions. Nonetheless, the emergent frequencies provide a good estimate and qualitative understanding of how the different biophysical parameters impact the frequency. Specifically, we are able to analytically calculate how motor density and elasticity (β), motor kinetics (Ψ and E s ), and fluid rheology (parameters B 1 to B 3 ) impact the neutral stability curve and frequencies. The predicted trends are line with experimental observations of more complicated systems with both sliding and bending deformations such as that shown in Figure 1A. While a direct comparison is not possible, it is heartening to note that significantly higher frequencies (increase of up to a factor of 10) can be achieved by varying the rheological parameters of the model fluids, such as by keeping viscosity via B 2 fixed while changing B 1 in the Kelvin-Voigt model and B 3 in the Maxwell model.
Our analysis also suggests that smooth modulation of the frequencies may be implemented once past the critical point but linear stability alone is insufficient to validate this. Furthermore the linear stability analysis does not show how the amplitude of the oscillations changes with parameters. The amplitude is Frontiers in Physics | www.frontiersin.org June 2022 | Volume 10 | Article 895536 9 dominated by non-linear terms. Thus to understand how nonlinearity affects frequencies and amplitudes far from critical points, we solved Equations 8-11 numerically using Matlab @ . Equations 8-11 constitute a stiff set of equations especially in the limit of small B 2 or large B 3 and hence the suite ODE 15s that solves stiff equations using a variable order method was employed. Sample time-dependent solutions are shown as follows: no drag case - Figure 3B, for Newtonian fluid solutions are shown in Figures 4A,B. Features of the nonlinear solution for the Kelvin-Voigt material is highlighted in Figures 5E,F that demonstrate important features of the full limit cycles. We also conducted systematic studies of how the amplitude and frequency varies for the fully non-linear solutions far from the critical point. Results are summarized in  Figures 8-10; here the focus was to study the characteristics and shape of these nonlinear solutions.

Supercritical Hopf Bifurcations and Fluid Limited Oscillation Amplitudes
Let us start with the two reference cases first. For the case of no fluid drag, the active energy input into the system by continuous motor activity is eventually dissipated away by internal motor friction with the anchoring spring as well as the passive links serving as intermediate storage agents. The motor activity derived internal friction is associated with the energy loss when motors detach and eventually go back to their rest length. We note that FIGURE 4 | (A) Phase-plots over complete cycles obtained from fully non-linear solutions for the Newtonian case are shown. The detachment function is cosh (2Y) and the motor activity parameter (density) β =100. We note that as B 2 increases from 10 −3 to 10 −1 , the amplitude of oscillations decreases in sync with smaller variations in N and in Y over a cycle. (B) When motor density and kinetic parameters (β and Ψ) are kept constant, the value of E s still controls if oscillations can ensure. Shown here are non-linear solutions-at B 2 10 1 -demonstrating the quenching of oscillations as E s is reduced from 1.7 to 1.6 (C,D) Neutral stability curves shown here for the Newtonian case when F cosh(E s Y ). Typically there are critical values of Ψ (for fixed β), or of β (for fixed Ψ) below which the base state is always stable.
Frontiers in Physics | www.frontiersin.org June 2022 | Volume 10 | Article 895536  That the bifurcating solutions arise via a supercritical Hopf-Poincare-Andronov bifurcation is deduced by the amplitude going to zero at criticality. For large deviations from the critical value of β. the amplitude decreases slowly. The inset shows sample steady periodic solutions at β =100≫18.9 (far from the critical point) for B 2 10 −3 (blue) and 1.0 (red). We note that as one goes far from the critical point, oscillations become highly asymmetric and the maximum value of Y attained deviates significantly from the minimum value. (Kelvin-Voigt, right) Here we plot Y max −1 as a function of β. Note that the amplitude of the oscillation tends to zero as one approaches the critical point. Comparison with the Newtonian case indicates a slower increase in the amplitude, as well as a lower magnitude of the amplitudes for large β. Examination of the time-series indicates that the frequencies (inset) are substantially higher than for the Newtonian case with the same value of B 2 and β.
Frontiers in Physics | www.frontiersin.org June 2022 | Volume 10 | Article 895536 oscillations can feature very sharp gradients in Y. Over a cycle, the attached motor fraction N can become very small reaching relative values of O (10 −4 ) as seen in Figure 4A (albeit for non-zero B 2 10 −3 ). Also rather large values of Y are attained relative to the stationary state (Y 0 = 1). Increasing E s enhances these features. When external viscosity starts to play a role and . In (C), we plot the (dimensionless) oscillation period. Note that the period for the KV case shown in (i) is typically lower than that for the Newtonian case (ii) indicating higher frequencies. becomes comparable to motor friction, the parameter range over which instability is exhibited becomes smaller. This is seen in Figure 4C where increasing B 2 is seen to push the neutral stability envelope to the right. When at a point where oscillations are the stable state, increasing B 2 is seen to reduce variations in (N, Y) (evident when plotted as a phase-plot, Figure 4A). For very large viscosity, oscillations become damped and then die out. Examination of the time-dependent behavior of N and Y just before and just after criticality confirms that the bifurcation is supercritical [37], with oscillations beginning with very small amplitudes. The farther the parameters are from their critical values, the larger the amplitude as is expected for a Hopf-Androvov-Poincare bifurcation. Thus external fluid viscosity, even in the purely Newtonian case, modulates the internal (mean-field) motor state captured by the two variables Y and N.
Introducing viscoelastic effects changes the dynamics in a variety of ways. When the ambient medium follows Kelvin-Voigt like response (a reasonable simple approximation to a gel), the ability to store energy increases for B 1 > 0. When 0 < B 2 ≪ 1, and with increasing B 1 surprisingly increases the minimum value of Ψ (for fixed β) or the minimum value of β (for fixed Ψ) for the onset of oscillations. There is also a significant effect on the frequency on emergent solutions that we have analytically derived. Let us consider what happens when the motor parameters Ψ, E s and β are held fixed, while B 1 and B 2 are allowed to vary ( Figure 5C). We note modest increases in the frequency with lower values of β yielding larger increments with increasing B 1 with B 1 0 reducing to the case of a pure Newtonian fluid medium. More dramatic increases in the frequency are seen when the kinetic parameter Ψ is allowed to vary and attain values larger than 1. Figure 5D illustrates a case with fixed B 2 where changing Ψ from 0.5 to 2.0 and simultaneously changing B 1 from 0 to 1 results in the frequency increasing by nearly a factor of 7. The amplitude of the oscillations however decreases with increasing B 1 with oscillations vanishing at sufficiently large values. Taken together, these results suggest that tuning the rheology of the fluid by adjusting B 1 allows for modest increases in frequencies with reduced solution amplitudes.
We further probe the variation in amplitudes away from the critical point and present results in Figures 6, 7. Figure 6 focusses on the motor variables (or the motor state) for both the Newtonian model and the Kelvin Voigt model. In both cases, we confirm that limit cycles emerge via supercritical bifurcations [37]. Figure 6A illustrates this feature; plotted here is the quantity |Y max (T) − Y 0 | and we find that this tends to 0 as β → β c . Emergent oscillations also onset, as required for this class of bifurcations, with fixed non-zero frequency. When β ≫ β c , the amplitude saturates. This saturation is likely due to a combination of two features. First, even though the number of active motors increases with increasing β, attached motors generate passive resistance when pulled backwards and limit the amplitude of the oscillations. Secondly, increased amplitudes result in a zipper effect-almost all motors rapidly disengage and detach together.   Figure 6B illustrates that the emergent oscillations when the ambient medium is a Kelvin-Voigt medium also arise from supercritical bifurcations. We note increased frequencies as well as slightly larger amplitudes. Figures 6A,B suggest that in the absence of any energy constraints (an open system where energy limitations are not present), limit cycles can in viscoelastic media can exhibit larger amplitudes and higher frequencies than in the Newtonian case. Until now we have focussed on how the motor variables Y and N are impacted by the fluid rheology dependent parameters. In experimental realizations of the filament motor system, one typically finds it easier to track filament related positions. For instance we may envision marking one end of the filament via florescent markers and tracking the position as a function of time. In Figure 7, we therefore illustrate how the filament displacement U varies over a period of a stable oscillation. For simplicity we treat all variables except β as constants and plot U max ≡max [U(T)], U min ≡min [U(T)], the characteristic amplitude |U max − U min | and the period T p defined implicitly via U (T + T p ) = U(T) for stable oscillations. For the Newtonian case, we fix B 2 0.1. For the Kelvin-Voigt case, we keep the viscous contribution the same but add an additional elastic component to the fluid response by setting B 1 1. The following features can deduced from Figures 7A,B. Fluid elasticity as expected, shifts the critical point at which oscillations emerge via supercritical bifurcations. Very close to onset the amplitudes compare well with the expected square root scaling with significant deviations seen as one goes to larger values of β as the oscillations become highly non-sinusoidal. Interestingly, we note that U min is nonmonotonic in β demonstrating significant non-linear effects. As far as the oscillation periods are concerned, our numerical solutions reveal that, consistent with the linear stability results, oscillations in the Newtonian fluid have larger time periods than in the KV fluid and thus lower frequencies. At β = 60, the value for the VE case is ≈ 10 while for the Newtonian case is ≈ 17. For β = 100, the difference is more enhanced with a period of ≈ 17 for the KV case and ≈ 29 for the Newtonian case. We also see some evidence of a qualitative difference in the way the time period changes with β between the two cases. For the Kelvin-Voigt case; the slope of the curve is smaller initially than at higher values of β, a trend missing in the case of motion in a Newtonian fluid.

Differences Between the Maxwell and Kelvin-Voigt Models
Numerical solutions for oscillations in a Maxwell fluid indicate that the onset of oscillations and frequencies may be more readily tuned than for the Kelvin-Voigt case. Oscillations are enabled at even O (1) values of β provided B 3 O(1); with these oscillations persisting for even large values of the viscosity B 2 . For fixed motor kinetics, in general, the higher the value of β, the lower B 3 needs to be to sustain oscillations. With Ψ, β and B 1 fixed, oscillations decrease in amplitude but increase in. frequency as B 2 increases. We hypothesize that the tuning of oscillations is made possible by changing the value of B 3 which is the ratio of the relaxation time for the Maxwell fluid and the mean time a motor is attached to the filament and is able to animate the filament. The viscosity on the other hand controls the energy dissipate rate and limits the amplitude of the oscillations. Figure 8 shows sample limit cycles seen for various parameter values. For B 3 1 and keeping parameters Ψ = 0.8, E s 2, B 1 0 fixed and restricting 0 < B 2 < 1, oscillations are seen for β = 5 but not for β = 50. Decreasing B 3 to 0.1 however moves the steady state to the linear stable phase space and thus oscillations are seen.
We have discussed the frequencies as well as the amplitudes (measured via suitable norms) of emergent oscillations until now. We conclude this section by focussing on the shape of the limit cycle and the asymmetries that may manifest. A motivation for looking at this in closer detail comes from the possibility of using synthetic versions of motor-filament systems as switching or triggering devices. One may imagine connecting this active viscoelastic unit in series to a dynamical system that uses the input dU/dT and/or U to provide downstream responses/signals. In Figures 9, 10 we examine the variation in the shape of the limit cycles for the specific case of the ambient fluid being a Maxwell fluid. In Figure 9 close to onset, the variations in Y(T) are smooth, and in fact almost sinusoidal as one moves towards the critical point. Increasing β makes the profile asymmetric; the mean motor extension attains large extensions (Y > 1) and compressions but the profile develops three extrema with two maxima and 1 minimum as seen by comparing the curves for β = 80 and β = 20. Figure 10 illustrates how the asymmetric double peaked shaped may be smoothened by increased viscosity (here by increasing the value of B 2 from 0.08 to 5.00.

SUMMARY AND CONCLUSION
Here, we studied a minimal system comprised of an active biofilament segment working against an elastic spring and immersed in model complex materials. A mean-field macroscale continuum description of filament-motor interactions was employed based on a cross-bridge model-however, other models can be readily used instead (see [28] for a general model that allows motors to undergo initial conformations as well as move relative to the filament). We used a combination of analytical and numerical methods to extract the conditions under which rheology of the fluid may be exploited to both shift the onset of instabilities, and modulate ensuing frequencies and amplitudes. We found that steady base states yield to stable oscillatory states or limit cycles. These solutions arise via supercritical Hopf-Andronov-Poincare' bifurcations with exponentially damped oscillation change to small limit cycle oscillation about the steady state. Linear stability results were validated and complemented with full time-dependent solutions that provide an indication of the manner in which the amplitude of these oscillations changes.
A shortcoming in our theoretical approach is that in the process of formulating analytically tractable equations, complex biochemical and biomechanical aspects are highly simplified. We have not considered noise and stochastic aspects of the problem that may be important for small to moderate values of β and/or low attachment rates [28]. While our model does not include motor driven filament bending which dominates in oscillatory dynamics of biological and synthetically devised ciliary structures, reduced dimensional and spatially dependent versions have been used to investigate oscillations of cilia in Newtonian fluids [8,9,29,31,36]. However, ciliary or flagellar motion in viscoelastic fluids remains an open challenge. Elastic and viscous responses of for a viscoelastic fluid made by adding polymers to a Newtonian solvent fluid may depend on concentration differently. Thus in a polymeric fluid the ratio B 1 /B 2 may scale with concentration with different exponents. Furthermore, the rheological behavior of a real polymeric fluid is more complex than that we have modeled. Nonetheless, the results here provide a foundation for more comprehensive analysis using detailed computational techniques and more realistic models for polymeric viscoelastic fluids.
The governing equations we utilized employed simplifications that allow for analytical and independent exploration of fluid viscoelastic properties (both its effective viscosity and elasticity, measured appropriately and encapsulated in parameters B 1 to B 4 ), motor kinetics and specifically the ratio of detachment and attachment constants (via parameter Ψ), effective importance of elastic (restoring) and viscous (dissipative) forces that includes possible coupling between the filament and its neighbors or anchoring substrates (through ratios B 1 and B 2 ), the ratio of viscoelastic relaxation times and the motor kinetic time-scales, and the strain dependent motor detachment rate (parameter S s . A further advantage of the simple model we have proposed is highlighted by examining the steady state solutions given by Eq. 5. We note that the steady state motor extension Y 0 and the steady state fraction of attached motors N 0 are independent of the rheological characteristics of the bulk medium. Thus coupling between the fluid and the motor aggregate is purely dynamical in nature. In other words, the aggregate is pre-strained internally by motor activity; but this pre-strained state is independent of the ambient fluid. It is thus possible to identify clearly changes in the dynamical state of the aggregate induced purely due to viscoelastic contributions, and to compare these with the results from the Newtonian case. Our results are also relevant to understanding the rheology and dynamical response of biological and synthetic muscles, and in devising actuators for soft-robots [38][39][40][41]. In general, active responses to imposed external stimuli or perturbations are analyzed using an effective frequency-dependent dynamic moduli in these systems. The minimal approach we have used here-modeling the ambient fluid in terms of minimal models with spring-like and dashpot-like elements-provide a clear means to evaluate these dynamic response functions (or transfer functions).
While direct comparison to ciliary experiments requires a spatially varying version of our model, there are experimental setups that can be used to check the theoretical predictions we have made. Motility assays offer an easy way to study the spatial dynamics of gliding microtubules or actin filaments. These filaments are animated by molecular motors that are grafted on the substrate and periodically attach to the filaments propelling them via tangential forces. Oscillations have been observed in filaments that are temporarily pinned in these systems with frequencies comparable to the values calculated in this paper for Newtonian fluids. In [18] we interpreted and analyzed bending oscillations for long filaments in Newtonian settings via a combination of linear stability analysis and nonlinear computations. By using 1) very short rigid (nonbendable) filaments, 2) model viscoelastic fluids as the nutrient media in which the gliding takes place, and 2) additional polymeric binders that can attach to the filament and play the role of passive adhesive springs [42,43], one can adapt current protocols and establish a setup that mimics closely the model analyzed in this paper. A natural next step would be the test of the model predictions from in these adapted gliding motility assays.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.