Biology is replete with examples of rigid and semi-flexible filaments, that singly or as aggregate bundles, interact with molecular motors in intracellular or extracellular settings. For instance, actin-myosin systems are crucial for cell function, development and growth [1–5] and in the transport of DNA [2]. Structured filament-motor systems comprised of microtubules and axonemal dynein motor arrays drive sperm locomotion [6–8] and algal motility [9]. These also constitute crucial components of healthy respiratory and reproductive tracts [10–15]. In the synthetic realm, motor assays and reconstituted biofilament networks frequently feature molecular motors translating, deforming and bending filaments [16–18].

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–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–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²⁺). 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 gel-like 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–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.

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 states – the kinetics of the mechanochemical cycle – are 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.

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 U∗ = L − L _o ; the speed of the plate is U ∗ ̇ . In the absence of fluid and solid inertia, the sum of forces on the filament must vanish. If F m ∗ is the effective active force density per attached motor, force balance on the filament provides the equation F d ∗ + ℓ A w A ρ p k p U ∗ + K eff U ∗ − ℓ A w A ρ m F m ∗ = 0 (1) The active force F 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 Y∗. 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 mean-field 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 N ̇ = ω on o 1 − N − ω off Y ∗ N . (2) Here the mean attachment rate ω on o is constant while the mean detachment rate ω off = ω off o F incorporates the state dependent function F ( Y ∗ ) that depends on motor extension (stretch). This quantity—the stretch of attached motors, Y ^∗—follows Y ∗ ̇ = − U ∗ ̇ + ω on o A − Y ∗ 1 − N N (3) Given N and Y ^∗, the effective active force per motor in Eq. 1 is given by F m ∗ = k m N Y ∗ . (4)

We next model the fluid interaction term, F d ∗ , in Eq. 1 by specifying constitutive laws for the manner in which the fluid interacts with the filament.

1. 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 ^∗/dt ^∗. An approximate equation for the drag force is then

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 / ω on o , and motor/filament extensions with A . Defining t ∗ = T ( ω on o ) − 1 , Y ∗ ≡ A Y and U ∗ ≡ A U 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 motor-filament system embedded in the viscoelastic medium. d N d T = 1 − N − Ψ F Y N , (8) d Y d T = − d U d T + 1 − Y 1 − N N (9) − F d = U − β N Y , w h e r e B 1 U + B 2 d U d T = B 3 d F d d T + B 4 F d (10) 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 B 1 U + B 2 d U d T + B 3 d 2 U d T 2 = B 3 d F d d T + B 4 F d 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 ) = ω off 0 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 F Y = cosh E s Y 2 , o r F Y = cosh E s Y . (11) Equations 8–11 admit a stationary state given by the set (N, Y, U) = (N ₀, Y ₀, U ₀) where these satisfy Y 0 = 1 , N 0 = 1 1 + Ψ F 0 , U 0 = β N 0 1 + B 1 / B 4 . (12) 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 ∼ 10⁻³ N/m, the effective linear density ρ _m w _A ∼ O (10⁸) m⁻¹, w _A ∼ 40−60 nm, k _N ∼ 16–100 pN μm⁻¹ and w _A ρ _N ∼ 10⁵–10⁷ m⁻¹ [1, 29]. The value of K _eff depends on the material. Microtubules are relatively stiff with large persistent lengths and the Young’s modulus E ∼ 1.2 GPa providing the stiffness per length K _eff ∼ Gw _A . For microtubules driven by dynein, previous studies report a persistence length ∼ 5 mm 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⁵ 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.

Here, we studied a minimal system comprised of an active bio-filament 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 ₀ and the steady state fraction of attached motors N ₀ 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–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 non-linear computations. By using 1) very short rigid (non-bendable) 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.