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.

1 Introduction

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].

FIGURE 1

FIGURE 1. (A) (Top) A schematic summarizing experimental observations of ciliary frequencies on bi-ciliated alga Chlamydomonas reinhardtii moving in model viscoelastic fluids [9] made by adding polymer to a Newtonian solvent. The observed ciliary beat frequency in polymeric fluids is much higher than in Newtonian fluids with the same viscosity. (Bottom) Schematic of the animated motor-filament segment that forms the basis of the minimal motor-filament system. Active motors (purple) may be in one of two states–attached to the segment (in green) or detached. Passive cross-linkers are shown in blue and are treated as linearly elastic springs. The resistance to motion of the test segment exerted by its neighbors are combined into an effective elastic resistance (spring constant K _eff). The ambient medium in which the assembly is embedded is viscoelastic. In the minimal model this medium may be treated as a Maxwell medium or as a Kelvin-Voigt medium. (B) Representations of the viscoelastic ambient media considered in this article—here the spring stiffness K (if present) for each material is related to the elastic modulus of the material comprising the medium. (i) A Kelvin-Voigt material, (ii) a Maxwell fluid, and (iii) a generalized standard solid element. The anti-Zener Jeffrey element in (iv) results in a more complicated relationship involving second derivatives of the strain rather than first derivatives as in (i)-(iii).

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.

2 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 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.

3 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 U∗ = L − L _o ; the speed of the plate is $\dot{U^{\ast }}$ . In the absence of fluid and solid inertia, the sum of forces on the filament must vanish. If $F_m^{\ast }$ is the effective active force density per attached motor, force balance on the filament provides the equation

$F_{\mathrm{d}}^{\ast }+{\ell }_A{w}_A{\rho }_{\mathrm{p}}{k}_{\mathrm{p}}{U}^{\ast }+K_{\text{eff}}{U}^{\ast }-{\ell }_A{w}_A{\rho }_m{F}_m^{\ast }=0(1)$

The active force $F_m^{\ast }$ 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

$\dot{N}={\omega }_{\text{on}}^o\left(1-N\right)-{\omega }_{\text{off}}\left(Y^{\ast }\right)N.(2)$

Here the mean attachment rate ${\omega }_{\text{on}}^o$ is constant while the mean detachment rate ${\omega }_{\text{off}}={\omega }_{\text{off}}^o\mathcal{F}$ incorporates the state dependent function $\mathcal{F}(Y^{\ast })$ that depends on motor extension (stretch). This quantity—the stretch of attached motors, Y ^∗—follows

$\dot{Y^{\ast }}=-\dot{U^{\ast }}+{\omega }_{\text{on}}^o\left(\mathcal{A}-Y^{\ast }\right)\left(\frac{1-N}{N}\right)(3)$

Given N and Y ^∗, the effective active force per motor in Eq. 1 is given by

$F_m^{\ast }=k_{m}N{Y}^{\ast }.(4)$

We next model the fluid interaction term, $F_d^{\ast }$ , 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

$F_{\mathrm{d},\mathrm{N}}^{\ast }-\left({\ell }_A{w}_A/{\ell }_N\right){\mu }_{\mathrm{N}}\dot{U^{\ast }}=0(5)$

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

${\sigma }_{\text{KV}}^{\ast }=E^{\text{KV}}{ϵ}_{\text{KV}}^{\ast }+{\mu }^{\text{KV}}{\dot{ϵ}}_{\text{KV}}^{\ast }$

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)

$F_{d,KV}^{\ast }={\ell }_A{w}_A{\sigma }_{\text{KV}}^{\ast }={\ell }_A{w}_A\left(E^{\text{KV}}/{\ell }_{\text{KV}}\right)U^{\ast }+{\ell }_A{w}_A\left({\mu }^{\text{KV}}/{\ell }_{\text{KV}}\right)\dot{U^{\ast }}(6)$

3. MaxwelL gel: For a viscoelastic Maxwell-like with a single relaxation time scale and a mesh length ℓ ^M, we have the relationship

$E_M{\mu }_M{\dot{ϵ}}_M^{\ast }=E_M{\sigma }_M^{\ast }+{\mu }_M{\dot{\sigma }}_M^{\ast }$

or

${\dot{ϵ}}_{\mathrm{M}}^{\ast }=\left(1/{\mu }^{\mathrm{M}}\right){\sigma }_{\mathrm{M}}^{\ast }+\left(1/E^{\mathrm{M}}\right){\dot{\sigma }}_{\mathrm{M}}^{\ast }.$

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

${\ell }_{\mathrm{A}}{w}_{\mathrm{A}}\dot{U^{\ast }}=\left({\ell }^{\mathrm{M}}/{\mu }^{\mathrm{M}}\right)F_{\mathrm{d},\mathrm{M}}^{\ast }+\left({\ell }^{\mathrm{M}}/E^{\mathrm{M}}\right){\dot{F}}_{\mathrm{d},\mathrm{M}}^{\ast }(7)$

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/{\omega }_{\text{on}}^o$ , and motor/filament extensions with $\mathcal{A}$ . Defining $t^{\ast }=T{({\omega }_{\text{on}}^o)}^{-1}$ , $Y^{\ast }\equiv \mathcal{A}Y$ and $U^{\ast }\equiv \mathcal{A}U$ allows us to recast (1), in dimensionless form. In this process, we find that re-scaling anchoring spring stiffness by defining ${\bar{K}}_{\text{eff}}\equiv K_{\text{eff}}/{\ell }_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.

$\frac{dN}{dT}=\left(1-N\right)-\mathrm{\Psi }\mathcal{F}\left(Y\right)N,(8)$

$\frac{dY}{dT}=-\frac{dU}{dT}+\left(1-Y\right)\left(\frac{1-N}{N}\right)(9)$

$-F_{\mathrm{d}}=U-\beta NY,\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\mathcal{B}}_{1}U+{\mathcal{B}}_2\frac{dU}{dT}={\mathcal{B}}_3\frac{d{F}_d}{dT}+{\mathcal{B}}_4{F}_d(10)$

with dimensionless parameter β quantifying motor density while dimensionless parameters ${\mathcal{B}}_1$ to ${\mathcal{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

${\mathcal{B}}_{1}U+{\mathcal{B}}_2\frac{dU}{dT}+{\mathcal{B}}_3\frac{d^{2}U}{d{T}^2}={\mathcal{B}}_3\frac{d{F}_d}{dT}+{\mathcal{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.

TABLE 1

TABLE 1. Motor-related parameters and scaled variables in the minimal model.

TABLE 2

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 $\mathcal{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: (a) a Newtonian fluid, (b) Kelvin-Voigt liquid-like viscoelastic medium, and (c) Maxwell solid-like medium. Models (b) and (c) serve as limiting cases of the more general linear Jeffreys’ model that belongs to the class of anti-Zener models and the standard solid model (SSM) of classical viscoelasticity. The ratio ${\mathcal{B}}_2/{\mathcal{B}}_1$ for the Kelvin-Voigt model is a ratio of rheological to motor kinetics time-scales. Parameter ${\mathcal{B}}_3$ has a similar interpretation for the Maxwell model. The anti-zener Jeffrey element results in a relationship involving the second derivative of the strain/displacement and results in the additional term with coefficient ${\mathcal{B}}_5$ , and is not analyzed in this paper.

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 ${\omega }_{\text{off}}(Y)={\omega }_{\text{off}}^0\mathcal{F}(Y)$ that modifies the bare detachment rate by a multiplicative factor $\mathcal{F}(Y)$ . For the analytical investigations, we restrict $\mathcal{F}(Y)$ to functions that is at least twice differentiable but otherwise general. For plotting, we will use the following forms

$\mathcal{F}\left(Y\right)=\cosh \left({\mathcal{E}}_s{Y}^2\right),\mathrm{o}\mathrm{r}\mathcal{F}\left(Y\right)=\cosh \left({\mathcal{E}}_{s}Y\right).(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=\left(\frac{1}{1+\mathrm{\Psi }{\mathcal{F}}_0}\right),U_0=\left(\frac{\beta N_0}{1+{\mathcal{B}}_1/{\mathcal{B}}_4}\right).(12)$

In Equation 11, the dimensionless parameter ${\mathcal{E}}_s$ quantifies the characteristic force attached motors can withstand before detachment. Note that the ratio ${\mathcal{B}}_1/{\mathcal{B}}_2$ compares the motor time scale with the relaxation time scale of the Kelvin-Voigt medium. Dimensionless parameter ${\mathcal{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 $\sim 5$ mm at room temperature, and for actin driven by myosin, $\sim 16\mu$ 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 $\sim 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 $\sim 1$ mPa.s to $\sim 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.

4 Results and Discussion

4.1 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–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.

$\frac{dN}{dT}=\left(1-N\right)-\mathrm{\Psi }\mathcal{F}\left({\mathcal{E}}_s,Y\right)N,(13)$

$\frac{dY}{dT}=-Z+\left(1-Y\right)\left(\frac{1-N}{N}\right)(14)$

When Z = 0, Equations 13, 14 admit steady state solutions (N _s , Y _s ) given by

$N_s\left(Z=0,\mathrm{\Psi },{\mathcal{E}}_s\right)={\left(1+\mathrm{\Psi }\mathcal{F}\left(1\right)\right)}^{-1},Y_s\left(Z=0,\mathrm{\Psi },{\mathcal{E}}_s\right)=1(15)$

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 ${\mathcal{E}}_s$ and for different functional forms of $\mathcal{F}$ the detachment function. Figure 2 shows an example case where we keep ${\mathcal{E}}_s$ and Ψ the same while using two forms of $\mathcal{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.

FIGURE 2

FIGURE 2. The mean motor extension Y and fraction of motors that are attached N as a function of Z steady solutions to Eqs 13, 14 for two forms of the detachment function (A) $\mathcal{F}=\cosh ({\mathcal{E}}_{s}Y)$ and (B) $\mathcal{F}=\cosh ({\mathcal{E}}_s{Y}^2)$ .

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].

4.2 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

$\left(N,Y,U,F_d\right)=\left(N_0,1,U_0,F_d^0\right)+ϵ\left(\bar{N},\bar{Y},\bar{U},{\bar{F}}_d\right).(16)$

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:

$\frac{d\bar{N}}{dT}=-\left(1+\mathrm{\Psi }{\mathcal{F}}_0\right)\bar{N}-\left(\mathrm{\Psi }{\mathcal{F}}_0^{\prime }{N}_0\right)\bar{Y}(17)$

$\frac{d\bar{Y}}{dT}=-\frac{d\bar{U}}{dT}-\left(\mathrm{\Psi }{\mathcal{F}}_0\right)\bar{Y}(18)$

$0=-\left({\mathcal{B}}_1+{\mathcal{B}}_4\right)\bar{U}-\left({\mathcal{B}}_2+{\mathcal{B}}_3\right)\frac{d\bar{U}}{dT}+\beta {\mathcal{B}}_3\left(N_0\frac{d\bar{Y}}{dT}+\frac{d\bar{N}}{dT}\right)+\beta {\mathcal{B}}_4\left(N_0\bar{Y}+\bar{N}\right)(19)$

Substituting next

$\left(\bar{N},\bar{Y},\bar{U}\right)=\exp \left(\sigma T\right)\left(\widehat{N},\widehat{Y},\widehat{U}\right)(20)$

in Eqs 17-19, we obtain coupled algebraic equations for $(\widehat{N},\widehat{Y})$ that determine the growth rate σ. Imposing solvability constraints, we get a cubic algebraic equation for σ:

$a{\sigma }^3+b{\sigma }^2+c\sigma +d=0(21)$

where

$\begin{array}{ccc}\hfill a& =\hfill & \left(\beta N_0+1\right){\mathcal{B}}_3+{\mathcal{B}}_2,d=\left(1+\mathrm{\Psi }{\mathcal{F}}_0\right)\mathrm{\Psi }{\mathcal{F}}_0\left({\mathcal{B}}_1+{\mathcal{B}}_4\right)\hfill \\ \hfill b& =\hfill & \left(\beta N_0+1\right){\mathcal{B}}_4+{\mathcal{B}}_1+\beta N_0{\mathcal{B}}_3\left(1+\mathrm{\Psi }{\mathcal{F}}_0\right)+\left({\mathcal{B}}_2+{\mathcal{B}}_3\right)\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)-\mathrm{\Psi }{\mathcal{F}}_0^{\prime }{N}_0\beta {\mathcal{B}}_3\hfill \\ \hfill c& =\hfill & \beta N_0{\mathcal{B}}_4\left(1+\mathrm{\Psi }{\mathcal{F}}_0\right)+\left({\mathcal{B}}_1+{\mathcal{B}}_4\right)\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)+\left({\mathcal{B}}_2+{\mathcal{B}}_3\right)\mathrm{\Psi }{\mathcal{F}}_0\left(1+\mathrm{\Psi }{\mathcal{F}}_0\right)-\mathrm{\Psi }{\mathcal{F}}_0^{\prime }{N}_0\beta {\mathcal{B}}_4\hfill \end{array}$

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 ${\mathcal{B}}_4=1$ , and are more interested in varying the rheological parameters of the fluid, the dependence of b on ${\mathcal{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 (σ ₁), σ ₂ = σ + iω and σ ₃ = σ − iω. When the base state is stable we have σ ₁ < 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.

4.2.1 Results From Stability Analysis

1. No drag from medium: Consider the simplest degenerate case where the medium does not exert any forces on the motor-filament system. In this case, we have ${\mathcal{B}}_1={\mathcal{B}}_2={\mathcal{B}}_3=0$ and ${\mathcal{B}}_4=1$ . The equation for the growth rates is a quadratic equation for this degenerate case of the form bσ ² + cσ + d = 0 where σ is the linear growth rate of infinitesimally small disturbances. For oscillatory instabilities, we require that c < 0 (so the growth rate, σ is non-negative); with this holding, c = 0 defines the locus of critical points. The neutral stability curve follows

$\beta \left(\mathrm{\Psi },{\mathcal{E}}_s\right)=\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right){\left[\left(\frac{\mathrm{\Psi }{\mathcal{F}}_0^{\prime }}{1+\mathrm{\Psi }{\mathcal{F}}_0}\right)-1\right]}^{-1}.(22)$

For constant Ψ and ${\mathcal{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 ${\mathcal{F}}_0^{\prime }>0$ and that the magnitude be large enough for the equation to be satisfied. Formally, we require $\mathcal{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

${\omega }^{\ast }=\sqrt{\frac{\left(1+\mathrm{\Psi }{\mathcal{F}}_0\right)\mathrm{\Psi }{\mathcal{F}}_0}{\beta N_0+1}}(23)$

2. Newtonian fluid: For the simplest drag-producing medium - a Newtonian fluid - we have ${\mathcal{B}}_1={\mathcal{B}}_3=0$ , ${\mathcal{B}}_2>0$ and ${\mathcal{B}}_4=1$ with specific forms listed in Table 1. The force balance reduces to

$-{\mathcal{B}}_2\frac{dU}{dT}-U+\beta NY=0$

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

${\omega }_N=\sqrt{\frac{\left(1+\mathrm{\Psi }{\mathcal{F}}_0\right)\mathrm{\Psi }{\mathcal{F}}_0}{\left(\beta N_0+1\right)+{\mathcal{B}}_2\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)}}=\sqrt{\frac{\beta +\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)+{\mathcal{B}}_2\mathrm{\Psi }{\mathcal{F}}_0\left(1+\mathrm{\Psi }{\mathcal{F}}_0\right)-\mathrm{\Psi }{\mathcal{F}}_0^{\prime }{N}_0\beta }{{\mathcal{B}}_2}}(24)$

Combining (Eqs 23) and (24) yields

$\left({\omega }_N/{\omega }^{\ast }\right)=\sqrt{\frac{\beta N_0+1}{\left(\beta N_0+1\right)+{\mathcal{B}}_2\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)}}$

and so,

$\left({\omega }_N/{\omega }^{\ast }\right)\approx 1-\frac{1}{2}{\mathcal{B}}_2\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right){\left(\beta N_0+1\right)}^{-1}\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}{\mathcal{B}}_2\ll 1$

while

$\left({\omega }_N/{\omega }^{\ast }\right)=\sqrt{\frac{1}{{\mathcal{B}}_2}}\sqrt{\frac{\beta N_0+1}{1+2\mathrm{\Psi }{\mathcal{F}}_0}}\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}{\mathcal{B}}_2\gg \left(1+\beta N_0\right).$

Thus for very high viscosity μ _N, the frequency scales as $1/\sqrt{{\mu }_{\mathrm{N}}}$ . The equation for the neutral stability curve simplifies to

$\beta \left[\left(\frac{\mathrm{\Psi }{\mathcal{F}}_0^{\prime }}{1+\mathrm{\Psi }{\mathcal{F}}_0}-1\right)\right]=\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)+{\mathcal{B}}_2\mathrm{\Psi }{\mathcal{F}}_0\left(1+\mathrm{\Psi }{\mathcal{F}}_0\right)\left(\frac{\beta N_0+{\mathcal{B}}_2\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)}{\beta N_0+1+{\mathcal{B}}_2\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)}\right).(25)$

From the analytical expression for the neutral stability curve, we deduce that for oscillatory instabilities to exist we require that ${\mathcal{F}}_0^{\prime }>0$ and be a suitably large number. For large ${\mathcal{B}}_2$ such that ${\mathcal{B}}_2\gg (\beta 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 ${\mathcal{E}}_s$ and with the form of $\mathcal{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 ${\mathcal{B}}_4=1$ , ${\mathcal{B}}_1>0$ , ${\mathcal{B}}_2>0$ while ${\mathcal{B}}_3=0$ . In this case again, oscillatory instabilities are possible at critical points provided ${\mathcal{F}}_0^{\prime }>0$ . The frequency of the bifurcating limit cycles at onset is obtained as

${\omega }_{KV}=\sqrt{\frac{\left(1+\mathrm{\Psi }{\mathcal{F}}_0\right)\mathrm{\Psi }{\mathcal{F}}_0\left({\mathcal{B}}_1+1\right)}{\left(\beta N_0+1\right)+{\mathcal{B}}_1+{\mathcal{B}}_2\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)}}.(26)$

For fixed motor parameters and fixed value of ${\mathcal{B}}_2$ , the frequency in a KV medium is thus higher than the purely viscous medium (with the same viscosity/dissipative components, ${\mathcal{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.

$\frac{{\omega }_{KV}}{{\omega }^{\ast }}=\sqrt{\frac{\left(\beta N_0+1\right)\left({\mathcal{B}}_1+1\right)}{\left(\beta N_0+1\right)+{\mathcal{B}}_1+{\mathcal{B}}_2\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)}}(27)$

$\mathrm{\Delta }{\omega }_{KV}=\left(\frac{{\mathcal{B}}_1}{{\mathcal{B}}_2}\right)\frac{\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)}{\left({\omega }_{KV}+{\omega }_N\right)}(28)$

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

$\beta \left[\left(\frac{\mathrm{\Psi }{\mathcal{F}}_0^{\prime }}{1+\mathrm{\Psi }{\mathcal{F}}_0}-1\right)\right]=\left(1+{\mathcal{B}}_1\right)\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)+{\mathcal{B}}_2\mathrm{\Psi }{\mathcal{F}}_0\left(1+\mathrm{\Psi }{\mathcal{F}}_0\right)\left(\frac{\beta N_0+{\mathcal{B}}_2\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)}{\beta N_0+1+{\mathcal{B}}_1+{\mathcal{B}}_2\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)}\right).(29)$

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, ${\mathcal{B}}_1=0$ , ${\mathcal{B}}_2>0$ , ${\mathcal{B}}_3>0$ and ${\mathcal{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 non-oscillatory instability for very high values of ${\mathcal{B}}_3$ for which b turns negative. We take this as evidence that the model is ill-posed for these high values of $\mathcal{B}{\left.\right]}_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 (σ ₁) 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

${\omega }_M=\sqrt{\frac{\left(1+\mathrm{\Psi }{\mathcal{F}}_0\right)\mathrm{\Psi }{\mathcal{F}}_0}{\left(\beta N_0+1\right)+{\mathcal{B}}_2\left(1+2\mathrm{\Psi }{\mathcal{F}}_0\right)+{\mathcal{B}}_3\left(1+\beta +2\mathrm{\Psi }{\mathcal{F}}_0-\mathrm{\Psi }{\mathcal{F}}_0^{\prime }{N}_0\beta \right)}}(30)$

Since ${\mathcal{B}}_3>0$ , the overall sign of the term multiply B₃ depends on the magnitude of the term $\mathrm{\Psi }{\mathcal{F}}_0^{\prime }{N}_0\beta$ 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 ₀. 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 ${\mathcal{B}}_3$ , the magnitude can be adjusted by solely changing the form of the detachment function (through ${\omega }_{\text{off}}^0$ and/or ${\mathcal{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 ${\mathcal{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 ${\mathcal{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 ${\mathcal{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 ${\mathcal{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.

4.3 Comparison to Non-linear Solutions

4.3.1 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 $\sim \sqrt{\beta -{\beta }_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 ${\mathcal{E}}_s$ ), and fluid rheology (parameters ${\mathcal{B}}_1$ to ${\mathcal{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 ${\mathcal{B}}_2$ fixed while changing ${\mathcal{B}}_1$ in the Kelvin-Voigt model and ${\mathcal{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 dominated by non-linear terms. Thus to understand how non-linearity 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 ${\mathcal{B}}_2$ or large ${\mathcal{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 non-linear 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 6, 7 for the Newtonian and Kelvin-Voigt model. Finally results for the Maxwell fluid are shown in Figures 8–10; here the focus was to study the characteristics and shape of these nonlinear solutions.

FIGURE 3

FIGURE 3. (A) Neutral stability curves separating linearly stable from unstable regions are shown for the no drag case. Plotted as solid curves are results for different values of ${\mathcal{E}}_s$ (as indicated, from 1.4–2). The dash-dot lines are the limiting values of Ψ below which the base state is always stable. The shaded yellow region indicates the stable region for Ψ =2. The dotted green region indicates linearly unstable space susceptible to oscillatory instabilities for Ψ =1.5. The yellow rhombus point indicated in (A) is unstable for some parameter values, as demonstrated by full non-linear solutions when Ψ =2 (shown alongside). (B) For even larger values of ${\mathcal{E}}_s>2$ , the curves shift rightward as Ψ increases. In all these cases, $\mathcal{F}({\mathcal{E}}_s,Y)=\cosh ({\mathcal{E}}_{s}Y)$ .

FIGURE 4

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 ${\mathcal{B}}_2$ increases from 10⁻³ to 10⁻¹, 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 ${\mathcal{E}}_s$ still controls if oscillations can ensure. Shown here are non-linear solutions—at ${\mathcal{B}}_2=10^1$ —demonstrating the quenching of oscillations as ${\mathcal{E}}_s$ is reduced from 1.7 to 1.6 (C,D) Neutral stability curves shown here for the Newtonian case when $\mathcal{F}=\cosh ({\mathcal{E}}_{s}Y)$ . Typically there are critical values of Ψ (for fixed β), or of β (for fixed Ψ) below which the base state is always stable.

FIGURE 5

FIGURE 5. Results for the case of a Kelvin-Voigt medium (A,B) Neutral stability curves from the linear stability analysis shown here for the Newtonian case when $\mathcal{F}=\cosh ({\mathcal{E}}_{s}Y)$ . Typically there are critical values of Ψ (for fixed β), or of β (for fixed Ψ) below which the base state is always stable. As plotted, the region on the concave side of the surface is linearly unstable. (C) and (D) Plots of the emergent frequencies (at onset) corresponding to Equation 26, also obtained from the linear stability results. (E) and (F) Time-dependent long-time solutions obtained from a full non-linear solution to the ode system. In (E) we show phase-plots in Y − N space that demonstrates how see oscillations exist only for a range of Ψ values (with other parameters fixed). Subplot (F) demonstrates how increasing ${\mathcal{B}}_1$ increases the frequency, and reduces the amplitude of the oscillations. Eventually at large values there is no instability.

FIGURE 6

FIGURE 6. Representative metrics for the fully non-linear periodic solutions for the (left) Newtonian, and (right) Kelvin-Voigt cases (Newtonian, left) We calculate the maximum value of Y _max over a cycle and plot Y _max −1 as a function of β − β _c where the critical value of beta for these parameters $({\mathcal{E}}_s=2,{\mathcal{B}}_1=0.1,\mathrm{\Psi }=0.8)$ is β _c ≈18.9. 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 ${\mathcal{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 ${\mathcal{B}}_2$ and β.

FIGURE 7

FIGURE 7. Non-linear solutions for Newtonian and Kelvin-Voigt media compared with the focus on the filament displacements U and period of oscillations. (A) Here we plot U _max and U _min, the maximum and minimum values attached by U(t) over a cycle when the oscillations have stabilized, as a function of β. Parameters are ${\mathcal{E}}_s=2$ and Ψ =0.8. The Newtonian case (N) corresponds to ${\mathcal{B}}_2=0.1$ while the Kelvin-Voigt (KV) case corresponds to ${\mathcal{B}}_2=0.1$ and ${\mathcal{B}}_1=1$ . (B) Amplitude of the filament oscillations |U_max − U_min| here plotted as a function of the deviation from the critical point β − β _c . Here the critical values β _c are different for the N and KV cases and may be read off from (A). 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.

FIGURE 8

FIGURE 8. Fully non-linear solutions for the Maxwell medium. Here to focus on the rheological parameters ${\mathcal{B}}_2$ and ${\mathcal{B}}_3$ , we study solutions for constant Ψ (=0.8) with $\mathcal{F}=\cosh ({\mathcal{E}}_{s}Y)$ (here ${\mathcal{E}}_s=2.0$ for all cases). Results are shown for increasing values of β (from left to right) (A–C) and increasing values of ${\mathcal{B}}_2$ (top to bottom).

FIGURE 9

FIGURE 9. Time series of the mean motor extension Y for the Maxwell case shown for β =20,40,60,80,100. Here we focus on the shape of the function close to the point where Y = Y _max. When β to close to its critical value, the oscillations are harmonic as seen from the form of Y(T) for β =20. As the value of β increases, the periodic profiles becomes highly asymmetric and also develop double maxima and double minima, especially evident for β =100. The high values of Y are associated with very low values of N when very few motors are attached, with these motors extended significantly.

FIGURE 10

FIGURE 10. The effect of increased dissipation (viscosity) on the form of Y(T) for the Maxwell model is illustrated here. The value of β is much larger than the critical value for the onset of oscillations. Fixing ${\mathcal{B}}_3$ at 0.1, and Ψ at 0.8, we investigate how the profiles change as ${\mathcal{B}}_2$ increases from 0.08 to 5.0. Note that each time-series has a different frequency and we have shown just a part of the period for each.

4.3.2 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 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⁻⁴) as seen in Figure 4A (albeit for non-zero ${\mathcal{B}}_2=10^{-3}$ ). Also rather large values of Y are attained relative to the stationary state (Y ₀ = 1). Increasing ${\mathcal{E}}_s$ enhances these features. When external viscosity starts to play a role and becomes comparable to motor friction, the parameter range over which instability is exhibited becomes smaller. This is seen in Figure 4C where increasing ${\mathcal{B}}_2$ is seen to push the neutral stability envelope to the right. When at a point where oscillations are the stable state, increasing ${\mathcal{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 ${\mathcal{B}}_1>0$ . When $0<{\mathcal{B}}_2\ll 1$ , and with increasing ${\mathcal{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 Ψ, ${\mathcal{E}}_s$ and β are held fixed, while ${\mathcal{B}}_1$ and ${\mathcal{B}}_2$ are allowed to vary (Figure 5C). We note modest increases in the frequency with lower values of β yielding larger increments with increasing ${\mathcal{B}}_1$ with ${\mathcal{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 ${\mathcal{B}}_2$ where changing Ψ from 0.5 to 2.0 and simultaneously changing ${\mathcal{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 ${\mathcal{B}}_1$ with oscillations vanishing at sufficiently large values. Taken together, these results suggest that tuning the rheology of the fluid by adjusting ${\mathcal{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 ₀| 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 ${\mathcal{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 ${\mathcal{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 non-monotonic 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 $\approx 10$ while for the Newtonian case is $\approx 17$ . For β = 100, the difference is more enhanced with a period of $\approx 17$ for the KV case and $\approx 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.

4.3.3 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 ${\mathcal{B}}_3=O(1)$ ; with these oscillations persisting for even large values of the viscosity ${\mathcal{B}}_2$ . For fixed motor kinetics, in general, the higher the value of β, the lower ${\mathcal{B}}_3$ needs to be to sustain oscillations. With Ψ, β and ${\mathcal{B}}_1$ fixed, oscillations decrease in amplitude but increase in. frequency as ${\mathcal{B}}_2$ increases. We hypothesize that the tuning of oscillations is made possible by changing the value of ${\mathcal{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 ${\mathcal{B}}_3=1$ and keeping parameters Ψ = 0.8, ${\mathcal{E}}_s=2$ , ${\mathcal{B}}_1=0$ fixed and restricting $0<{\mathcal{B}}_2<1$ , oscillations are seen for β = 5 but not for β = 50. Decreasing ${\mathcal{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 ${\mathcal{B}}_2$ from 0.08 to 5.00.

5 Summary and Conclusion

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 ${\mathcal{B}}_1/{\mathcal{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 ${\mathcal{B}}_1$ to ${\mathcal{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 ${\mathcal{B}}_1$ and ${\mathcal{B}}_2$ ), the ratio of viscoelastic relaxation times and the motor kinetic time-scales, and the strain dependent motor detachment rate (parameter ${\mathcal{S}}_s\left.\right)$ . 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.

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.

Author Contributions

The project was conceived by AG. Theoretical calculations were done by AM and AG. Numerical analysis was conducted by AM and JT. All authors contributed to the final form of the manuscript.

Funding

AG acknowledges funding from the National Science Foundation via awards NSF-CBET-2047210 and NSF-MCB-2026782. JT was supported by funding from NSF-MCB-2026782.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2022.895536/full#supplementary-material

References

1. Howard J. Mechanics of Motor Proteins and the Cytoskeleton. Sunderland: Sinauer Associates (2001).

Google Scholar

2. Grill SW, Kruse K, Jülicher F. Theory of Mitotic Spindle Oscillations. Phys Rev Lett (2005) 94(10):108104. doi:10.1103/physrevlett.94.108104

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Bizzi E, Chapple W, Hogan N. Mechanical Properties of Muscles: Implications for Motor Control. Trends Neurosciences (1982) 5:395–8. doi:10.1016/0166-2236(82)90221-1

CrossRef Full Text | Google Scholar

4. Hogan N. Adaptive Control of Mechanical Impedance by Coactivation of Antagonist Muscles. IEEE Trans Automat Contr (1984) 29:681–90. doi:10.1109/TAC.1984.1103644

CrossRef Full Text | Google Scholar

5. Walcott S. Muscle Activation Described with a Differential Equation Model for Large Ensembles of Locally Coupled Molecular Motors. Phys Rev E (2014) 90:042717. doi:10.1103/PhysRevE.90.042717

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Witman GB. Introduction to Cilia and Flagella. In: RA Bloodgood, editor. Ciliary and Flagellar Membranes. New York: Plenum (1990). p. 1–30. doi:10.1007/978-1-4613-0515-6_1

CrossRef Full Text | Google Scholar

7. Machin KE. The Control and Synchronization of Flagellar Movement. Proc Roy Soc B (1963) 158(970):88–104.

Google Scholar

8. Brokaw CJ. Molecular Mechanism for Oscillation in Flagella and Muscle. Proc Natl Acad Sci U.S.A (1975) 72(8):3102–6. doi:10.1073/pnas.72.8.3102

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Qin B, Gopinath A, Yang J, Gollub JP, Arratia PE. Flagellar Kinematics and Swimming of Algal Cells in Viscoelastic Fluids. Sci Rep (2015) 5:9190. doi:10.1038/srep09190

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Sears PR, Thompson K, Knowles MR, Davis CW. Human Airway Ciliary Dynamics. Am J Physiology-Lung Cell Mol Physiol (2013) 304:L170–L183. doi:10.1152/ajplung.00105.2012

CrossRef Full Text | Google Scholar

11. Kempeneers C, Seaton C, Chilvers MA. Variation of Ciliary Beat Pattern in Three Different Beating Planes in Healthy Subjects. Chest (2017) 151:993–1001. doi:10.1016/j.chest.2016.09.015

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Knowles MR, Leigh MW, Carson JL, Davis SD, Dell SD, Ferkol TW, et al. Mutations ofDNAH11in Patients with Primary Ciliary Dyskinesia with normal Ciliary Ultrastructure. Thorax (2012) 67:433–41. doi:10.1136/thoraxjnl-2011-200301

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Papon J-F, Bassinet L, Cariou-Patron G, Zerah-Lancner F, Vojtek A-M, Blanchon S, et al. Quantitative Analysis of Ciliary Beating in Primary Ciliary Dyskinesia: a Pilot Study. Orphanet J Rare Dis (2012) 7:78. doi:10.1186/1750-1172-7-78

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Thomas B, Rutman A, O'Callaghan C. Disrupted Ciliated Epithelium Shows Slower Ciliary Beat Frequency and Increased Dyskinesia. Eur Respir J (2009) 34:401–4. doi:10.1183/09031936.00153308

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Ringers C, Olstad EW, Jurisch-Yaksi N. The Role of Motile Cilia in the Development and Physiology of the Nervous System. Phil Trans R Soc B (2019) 375:20190156. doi:10.1098/rstb.2019.0156

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Sanchez T, Welch D, Nicastro D, Dogic Z. Cilia-like Beating of Active Microtubule Bundles. Science (2001) 333(6041):456–9. doi:10.1126/science.1203963

CrossRef Full Text | Google Scholar

17. Sanchez T, Chen DTN, DeCamp SJ, Heymann M, Dogic Z. Spontaneous Motion in Hierarchically Assembled Active Matter. Nature (2012) 491(7424):431–4. doi:10.1038/nature11591

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Fily Y, Subramanian P, Schneider TM, Chelakkot R, Gopinath A. Buckling Instabilities and Spatio-Temporal Dynamics of Active Elastic Filaments. J R Soc Interf (2020) 17(165):20190794. doi:10.1098/rsif.2019.0794

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Vilfan A, Frey E. Oscillations in Molecular Motor Assemblies. J Phys Condens Matter (2005) 17(47):S3901–S3911. doi:10.1088/0953-8984/17/47/018

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Nguyen KD, Sharma N, Venkadesan M. Active Viscoelasticity of Sarcomeres. Front Robot AI (2018) 5:69. doi:10.3389/frobt.2018.00069

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Vilfan A, Duke T. Synchronization of Active Mechanical Oscillators by an Inertial Load. Phys Rev Lett (2003) 91:114101. doi:10.1103/physrevlett.91.114101

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Endow SA, Higuchi H. A Mutant of the Motor Protein Kinesin that Moves in Both Directions on Microtubules. Nature (2000) 406:913–6. doi:10.1038/35022617

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Badoual M, Jülicher F, Prost J. Bidirectional Cooperative Motion of Molecular Motors. Proc Natl Acad Sci U.S.A (2002) 99:6696–701. doi:10.1073/pnas.102692399

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Grill SW, Gönczy P, Stelzer EHK, Hyman AA. Polarity Controls Forces Governing Asymmetric Spindle Positioning in the Caenorhabditis elegans Embryo. Nature (2001) 409:630–3. doi:10.1038/35054572

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Colombo K, Grill SW, Kimple RJ, Willard FS, Siderovski DP, Gönczy P. Translation of Polarity Cues into Asymmetric Spindle Positioning in Caenorhabditis elegans Embryos. Science (2003) 300:1957–61. doi:10.1126/science.1084146

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Grill SW, Howard J, Schäffer E, Stelzer EHK, Hyman AA. The Distribution of Active Force Generators Controls Mitotic Spindle Position. Science (2003) 301:518–21. doi:10.1126/science.1086560

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Li C, Qin B, Gopinath A, Arratia PE, Thomases B, Guy RD. Flagellar Swimming in Viscoelastic Fluids: Role of Fluid Elastic Stress Revealed by Simulations Based on Experimental Data. J R Soc Interf (2017) 14(135):20170289. doi:10.1098/rsif.2017.0289

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Gopinath A, Chelakkot R, Mahadevan L. Filament Extensibility and Shear Stiffening Control Persistence of Strain and Loss of Coherence in Cross-Linked Motor-Filament Assemblies. bioRxiv (2020). doi:10.1101/423582

CrossRef Full Text | Google Scholar

29. Camalet S, Jülicher F. Generic Aspects of Axonemal Beating. New J Phys (2000) 2:24.1–24.23. doi:10.1088/1367-2630/2/1/324

CrossRef Full Text | Google Scholar

30. Patteson AE, Gopinath A, Goulian M, Arratia PE. Running and Tumbling with E. coli in Polymeric Solutions. Sci Rep (2015) 5(1):15761–11. doi:10.1038/srep15761

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Hamilton E, Pellicciotta N, Feriani L, Cicuta P. Motile Cilia Hydrodynamics: Entrainment versus Synchronization when Coupling through Flow. Phil Trans R Soc B (2019) 375:20190152. doi:10.1098/rstb.2019.0152

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Sangani A, Gopinath A. Elastohydrodynamical Instabilities of Active Filaments, Arrays and Carpets Analyzed Using Slender Body Theory. Phys. Rev. Fluids (2020). 5(8):083101.

CrossRef Full Text | Google Scholar

33. Fatehiboroujeni S, Gopinath A, Goyal S. Nonlinear Oscillations Induced by Follower Forces in Prestressed Clamped Rods Subjected to Drag. J Comp Non Dyn (2018) 13(12):121005. doi:10.1115/1.4041681

CrossRef Full Text | Google Scholar

34. Fatehiboroujeni S, Gopinath A, Goyal S. Follower Forces in Pre-stressed Fixed-Fixed Rods to Mimic Oscillatory Beating of Active Filaments. In: Proceedings of the ASME 2018 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. Volume 6: 14th International Conference on Multibody Systems, Nonlinear Dynamics, and Control; August 26–29, 2018; Quebec City, Quebec, Canada (2018). ASME DETC2018-85449, V006T09A033. doi:10.1115/detc2018-85449

CrossRef Full Text | Google Scholar

35. Chelakkot R, Gopinath A, Mahadevan L, Hagan MF. Flagellar Dynamics of a Connected Chain of Active, Polar, Brownian Particles. J R Soc Interf (2014) 11(92):20130884. doi:10.1098/rsif.2013.0884

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Chakrabarti B, Saintillan D. Spontaneous Oscillations, Beating Patterns, and Hydrodynamics of Active Microfilaments. Phys Rev Fluids (2019) 4:043102. doi:10.1103/physrevfluids.4.043102

CrossRef Full Text | Google Scholar

37. Strogatz SH. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. Reading, Mass: Addison-Wesley (1994).

Google Scholar

38. Anderson IA, Gisby TA, McKay TG, O’Brien BM, Calius EP. Multi-functional Dielectric Elastomer Artificial Muscles for Soft and Smart Machines. J Appl Phys (2012) 112:041101. doi:10.1063/1.4740023

CrossRef Full Text | Google Scholar

39. Buerger SP, Hogan N. Complementary Stability and Loop Shaping for Improved Human-Robot Interaction. IEEE Trans Robot (2007) 23:232–44. doi:10.1109/TRO.2007.892229

CrossRef Full Text | Google Scholar

40. George NT, Irving TC, Williams CD, Daniel TL. The Cross-Bridge Spring: Can Cool Muscles Store Elastic Energy? Science (2013) 340:1217–20. doi:10.1126/science.1229573

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Hines L, Petersen K, Lum GZ, Sitti M. Soft Actuators for Small-Scale Robotics. Adv Mater (2017) 29:1603483–43. doi:10.1002/adma.201603483

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Sharma P, Ghosh S, Bhattacharya S. Microrheology of a Sticking Transition. Nat Phys (2008) 4:960–6. doi:10.1038/nphys1105

CrossRef Full Text | Google Scholar

43. Mani M, Gopinath A, Mahadevan L. How Things Get Stuck: Kinetics, Elastohydrodynamics, and Soft Adhesion. Phys Rev Lett (2012) 108(22):226104. doi:10.1103/physrevlett.108.226104

PubMed Abstract | CrossRef Full Text | Google Scholar