A unified model for the dynamics of ATP-independent ultrafast contraction

Despite the model’s simplicity, we find that it reproduces the contraction kinematics of Spirostomum and Vorticella semiquantitatively. We illustrate this by comparing the trajectories of normalized length λ(T)≡(X(L,T)−X(0,T))/L from experimental measurements to fitted model curves. Contraction events for Spirostomum and Vorticella are shown in SI Appendix, Movies S1 and S2, and the experimental methods are described in Materials and Methods, Experimental methods. The second-order nonautonomous PDE model does not have an exact closed-form solution, though we describe exact series solutions in two limiting cases in Results, section C. To fit the model, we numerically integrate the PDE with different choices of parameters and compare the results to the experimental traces. See Materials and Methods, Computational methods for details on the computational methods. In Fig. 2A, we show five contraction measurements in Spirostomum along with a model solution fitted to their average, illustrating the feasibility of our model in capturing the contraction dynamics. The corresponding contraction velocities |∂_Tλ| are shown in Fig. 2B. We also fit the model with a fixed boundary condition at S = L to data for five Vorticella contraction events (Fig. 2 C and D).

In SI Appendix, Figs. S6 and S7, we present fits for each of the ten individual trajectories shown in Fig. 2 A and C, demonstrating that we can obtain excellent fits for all trajectories using fairly consistent parameters for each organism. As we explain in Results, section C, we find that these contraction events occur in a regime where inertia is unimportant, meaning that the dynamics are approximately determined only by the ratio μ/α. We fit the simplified “slow-wave” model (cf. Eq. 10) and obtain a typical set of parameters for Spirostomum: L = 0.58 mm, g_min = 0.35, μ/α = 5 ms/mm², V = 0.8 mm/ms, and W = 0.18 mm. For Vorticella, a typical set of parameters is L = 0.27 mm, g_min = 0.32, μ/α = 1 ms/mm², V = 0.23 mm/ms, and W = 0.08 mm. In SI Appendix, Supplementary Results, section B, we demonstrate that the fitted parameters for Vorticella are consistent in order of magnitude with physiological parameters used in a detailed computational model in ref. 49. According to a two-tailed t test, the distributions of quantities μ/α, V, and W each have means that are statistically different for the two organisms (t test P-value ≪0.05). In the following section, we investigate the kinetic factors that underlie the rate of chemical driving captured by the parameters V and W.

The speed V and, less intuitively, the width W of the active stretch wave g(S − VT; W) determine how rapidly the length changes over time. As W decreases, the active stretch wave approaches a traveling step function in which the springs ahead of the wave have their original rest lengths, and those behind have rest lengths which are reduced by a factor g_min. As W increases, the distribution of changing rest lengths widens, and the driving and relaxation become more spatially and temporally homogeneous. We illustrate this in Fig. 2E, where we show the length trajectory for several values of W.

We are interested in understanding which features of the Ca²⁺ dynamics control the quantities V and W, since these parameters control the rate at which the mechanical system is chemically driven out of equilibrium. To do this, we analyzed a kinetic model (58) for a traveling Ca²⁺ wave powered by calcium-induced calcium release (62, 63) from an out-of-equilibrium internal store [e.g., the endoplasmic reticulum (30, 50)] into the cytosol. It has been previously shown that in this model, the wave speed V approximately follows (58)

where J represents the Ca²⁺ influx of open channels, D_C is the diffusion constant of Ca²⁺, and C₀ is the concentration threshold of the Ca²⁺ channels (see Materials and Methods, Ca²⁺ dynamics and SI Appendix, Supplementary Methods, section B for more details). We incorporated the binding of cytosolic Ca²⁺ to myoneme fibers into the model and derived a relationship for the width W of the resulting active stretch wave:

where K₊ is the rate at which Ca²⁺ binds to myonemes. We provide numerical confirmation of this relationship in SI Appendix, Supplementary Methods, section B and the Inset of Fig. 2E. A seemingly restrictive assumption made in the derivation of Eq. 6 is that the total concentration of myoneme binding sites B_tot is much smaller than the typical cytosolic Ca²⁺ concentration. However, we find numerically that the predicted scaling holds even as B_tot increases to relatively large values.

This analysis reveals that features in the dynamics of contraction, such as the time taken for the length to change, can inform us about the magnitudes of the underlying biophysical parameters like the binding rate K₊. Assuming D_C remains the same in both organisms [though the presence of internal buffers could affect this (58)], the estimated values of V and W for Vorticella and Spirostomum suggest, through Eqs. 5 and 6, that these organisms possess different threshold Ca²⁺ concentrations C₀ (30, 64, 65) and different release rates J or binding constants K₊. This minimal model serves as a first step in connecting contraction dynamics to the organisms’ biochemistry, but more in vivo experiments are necessary to test these predictions.

The model considered so far (Eq. 2 together with 22) is written in dimensional quantities and depends on the parameters α and μ, which are related to the stiffness of the material and viscosity of the environment, and V, W, and g_min, which are the speed, width, and minimum value of the active stretch wave. However, these parameters are not all independent. We now reduce the number of parameters by one through nondimensionalization. We use lowercase variables throughout to denote dimensionless quantities.

There are two length scales, W and L, and we pick the uncontracted length L as the reference length and define s = S/L, x = X/L, and w = W/L. Next, three timescales can be defined:

Here, τ_w is the wave timescale, which represents the time taken for the center of the traveling wave to traverse the length L; τ_m is the material timescale, which represents the time taken for a sound wave to traverse the length of the spring chain; τ_d is the drag timescale, characterizing the exponential deceleration of motion due to viscous damping. Unless otherwise specified, we nondimensionalize time using t = T/τ_w. In these units, the wave speed is equal to 1, and Eq. 2 is expressed as the following PDE for x(s, t):

where f^*(s)≡g^*(Ls) is the active stretch wave profile, rewritten to take a dimensionless argument s. The size of τ_w relative to the two remaining timescales defines the nondimensional parameters

The parameter η_{w, m} describes the ratio of the myonemal speed of sound (which depends on the stiffness and density) to the wave traversal speed, and it is independent of length. The parameter η_{w, d} describes the ratio of the drag timescale (depending on viscosity and density) to the wave traversal speed and is proportional to the length.

We can further simplify these nondimensional dynamics by taking the “slow wave” limit η_{w, m}, η_{w, d} ≫ 1, in which (8) reduces to

where the only remaining parameter isEq. 10 can be viewed as an overdamped limit of the model, in which inertia has dropped out of the dynamics because only a first derivative with respect to time remains, and in which mass is no longer relevant because ζ is independent of $\tilde{m}$.

Eq. 10 is similar to the diffusion (or heat) equation with a source term. One difference is the boundary conditions: typically, the heat equation has fixed Dirichlet or Neumann boundary conditions, yet x obeys a boundary condition that depends in time on f^*(s − t) (cf. Eq. 3). However, we can define

and express Eq. 10 as

The boundary conditions are now the fixed Dirichlet conditions h(0, t)=h(1, t)=0, and thus Eq. 13 is exactly the diffusion equation with a source term −∂_tf, for which exact series solutions can be found. The mixed boundary conditions of Vorticella also result in mixed boundary conditions for the diffusion equation.

We also consider a “quench” limit η_{w, d}, η_{w, m} ≪ 1, in which the traveling wave moves at infinite speed, such that the rest lengths of all the springs are set to their new values before the organism can react. This implies that f(s) is independent of time. Because the timescale τ_w is zero if the wave is infinitely fast, we instead nondimensionalize time using τ_m in the quench limit. The dynamics consist of only the relaxation of the system to its equilibrium

The deviation from equilibrium, y(s, t)≡x(s, t)−x^eq(s), obeys the dynamics

with the initial condition y(s, 0)) ≡ x(s, 0)−x^eq(s), ∂_ty(s, 0)=0, and stress-free boundary conditions ∂_sy(0, t)=∂_sy(1, t)=0. We note that the one remaining parameter

is independent of wave speed, which is infinite. The quench limit has an exact Fourier series solution, described in SI Appendix, Supplementary Methods, section D.

In Fig. 3 A and B we compare the solutions of the full nondimensional model (Eq. 8), the slow wave model (Eq. 10), and the quench model (Eq. 15) for two sets of the nondimensional parameters η_{w, d} and η_{w, m}. In these plots, we fix w = 0.1 and g_min = 0 for simplicity. For η_{w, d} = η_{w, m} = 1 (Fig. 3A), the models all behave quite differently because neither limiting case is reached. Decaying oscillations about l = 0 are observed for the full and quench model, while the slow wave model does not oscillate because it lacks inertia. The oscillations for the quench and full model have the same frequency (which is proportional to η_{m, d}, reflecting the resonant frequency of the spring chain) yet different waveforms and amplitudes. For η_{w, d} = 20, η_{w, m} = 2 (Fig. 3B), the full model and slow wave model are in close agreement.

A stark feature of the full model compared to the quench model is the asymmetry of contraction: In the quench model, all springs instantly have new rest lengths, so all springs relax simultaneously and symmetrically about the center s = 0.5. This is not the case for the full model, where the wave first arrives at s = 0, causing the springs to begin relaxing before those at s = 1. This asymmetry of contraction can be observed in the experimental images of Spirostomum and Vorticella in Fig. 1 A and B, and we further illustrate this in Fig. 4A –C. In Fig. 4D –G we show x(s, t) at several values of t for the full and quench models corresponding to the length trajectories shown in Fig. 3 A and B.

From the fits to the experimental data of Vorticella and Spirostomum contracting in water (see SI Appendix, Supplementary Results, section B), neither organism is expected to obey the quench-like dynamics, being closer to the slow-wave limit. This implies that the organisms are relaxing appreciably while the Ca²⁺ wave traverses their lengths.

An additional significant difference between the different models is the scaling of the maximum contraction rate with the viscosity. Notably, this has been studied in previous experimental and computational works on Vorticella and Spirostomum, against which we may validate our model (9, 42, 49).

To explore this scaling behavior, we computed the maximum contraction rate ${\dot{\lambda }}_{}\text{max}\equiv {\text{max}}_t|{\partial }_t\lambda (t)|$ as we varied the nondimensional parameters η_{w, m} and η_{w, d} in the full, quench, and slow wave models. The results are displayed in Fig. 3 C and D, again setting w = 0.1 and g_min = 0 for simplicity. The dynamics of the quench model only depend on the dimensionless parameter η_{m, d} (cf. Eq. 15), and the maximum rate (in units of $L/{\tau }_m=\sqrt{\alpha }$) asymptotically scales like ${\dot{\lambda }}_{}\text{max}\sim {\eta }_{m,d}^{-1}={\eta }_{w,m}/{\eta }_{w,d}=\sqrt{\alpha }/L\mu$. This can be qualitatively understood as follows: as shown in SI Appendix, Supplementary Methods, section D, an exact solution for the quench model can be constructed by decomposing the initial condition into spatial Fourier modes and solving for the decay of each mode. Each mode acts like a mass, attached to a wall by a stretched spring with an effective drag constant η_{m, d}, which is suddenly released (see SI Appendix, Eq. 66). By Stokes’ law, the maximum rate should then scale inversely with η_{m, d} as observed.

For both the full and slow wave models, ${\dot{\lambda }}_{}\text{max}\sim {\eta }_{w,m}$ up to a point at which ${\dot{\lambda }}_{}\text{max}$ plateaus at 1. This results from the contraction speed not significantly exceeding the speed of the traveling active stretch wave, which is defined as 1 in these units. Interestingly, for both models ${\dot{\lambda }}_{}\text{max}$ (in units of L/τ_w = V) asymptotically scales like ${\dot{\lambda }}_{}\text{max}\sim {\eta }_{w,d}^{-1/2}$. To understand this, recall that the slow wave model can be formally mapped to the diffusion equation with a source term (cf. Eqs. 10 and 13). The dynamics are thus characterized by the “diffusion constant” ζ = η_{w, m}²/η_{w, d} = α/VμL. For diffusion, it is known that an initial perturbation to the system spreads out over time with a length scale that grows like $\sqrt{\zeta t}$, which implies that the corresponding rate scales like ${\dot{\lambda }}_{}\text{max}\sim {\zeta }^{1/2}=\sqrt{\alpha /V\mu L}$. This fact explains both the observed scaling ${\dot{\lambda }}_{}\text{max}\sim {\eta }_{w,m}$ up to its saturation at 1, and the asymptotic scaling ${\dot{\lambda }}_{}\text{max}\sim {\eta }_{w,d}^{-1/2}$. We emphasize that in Fig. 3, we have assumed the wave width W = 0.1 L scales with the length L; as demonstrated in SI Appendix, Fig. S2, different scaling behaviors result if one varies the uncontracted length L while holding the wave width W fixed. In SI Appendix, Fig. S4 we show that this predicted scaling fits our experimental data for Spirostomum.

In refs. 42 and 49, the scaling of ${\dot{\lambda }}_{}\text{max}$ with the solution viscosity ν for Vorticella is measured to be ${\dot{\lambda }}_{}\text{max}\sim {\nu }^{-1/2}$ (Fig. 5A). This is consistent with our observed scaling ${\dot{\lambda }}_{}\text{max}\sim {\zeta }^{1/2}$ for the full and slow wave models when holding all other parameters constant and only varying μ (which is proportional to ν). Our model thus exhibits the scaling observed in existing experimental data for Vorticella. We note that an alternative description for the ν^−1/2 scaling based on power-limited contraction was provided in ref. 42. Our explanation here does not involve a physical argument about power conservation and follows solely from the dynamical equations. It may be possible to show that the energy-based argument is equivalent to the dynamical one presented here, but doing so would likely require a more detailed treatment of thermodynamics of the Ca²⁺ wave (as in ref. 67), and we leave this to future work.

Scaling data for Spirostomum are available in refs. 8 and 9 (Fig. 5B). In these experiments, to aid imaging, the cells are confined in either a narrow agar channel (blue points) or between a microscope slide and a glass coverslip (orange points). Depending on the confinement method, the scaling appears consistent with either ${\dot{\lambda }}_{}\text{max}\sim {\nu }^{-1}$, corresponding to the quench regime, or ${\dot{\lambda }}_{}\text{max}\sim {\nu }^{-1/2}$. Our fits in Fig. 2, together with observation of the asymmetric contraction dynamics (Fig. 5A), suggest that Spirostomum should not be in the quench regime. The quench-like scaling may be an artifact of using pressure on the coverslips to prevent the organisms from swimming out of view; we speculate that this pressure could hinder contraction during the traversal of the wave, effectively quenching the springs. The imaging rate in refs. 8 and 9 (800 Hz) may also be too low to resolve the maximum in the contraction rate profile precisely. We thus performed our own measurements of the contraction of free-swimming Spirostomum using a high-speed camera, recording at a minimum of 2,000 Hz. The results, displayed in Fig. 5C, are consistent with the predicted scaling ${\dot{\lambda }}_{}\text{max}\sim {\nu }^{-1/2}$, though we cannot conclusively exclude the scaling ${\dot{\lambda }}_{}\text{max}\sim {\nu }^{-1}$. We note that the low ${\dot{\lambda }}_{}\text{max}$ at 1 cP is consistent with the predicted roll-off at low viscosity shown in Fig. 3D.

The roll-off corresponds to an approach toward a third dynamic regime, which occurs when the active stretch wave is sufficiently slow, or the springs can sufficiently overcome the drag resistance to follow the active stretch wave very closely, with little deviation from equilibrium. In this “quasi-static” limit, there is no dependence of ${\dot{\lambda }}_{}\text{max}$ on either the stiffness or drag since the contraction speed is solely determined by the speed of the active stretch wave. This regime occurs in the slow wave model when ζ ≫ 1, explaining the approach toward ${\dot{\lambda }}_{}\text{max}=1$ for large η_{w, m} and small η_{w, d} in Fig. 3 C and D. The three scaling regimes—the quench regime (${\dot{\lambda }}_{}\text{max}\sim {\nu }^{-1}$), the intermediate regime (${\dot{\lambda }}_{}\text{max}\sim {\nu }^{-1/2}$), and the quasi-static regime (${\dot{\lambda }}_{}\text{max}\sim {\nu }^0$)—fully characterize the qualitative dynamical behavior of myoneme contraction.

To summarize the dynamical behavior of the system as a function of η_{w, d} and η_{w, m}, in SI Appendix, Supplementary Methods, section E we introduce the dimensionless quantity χ, which measures the peak elastic energy stored in the system as a fraction of the maximum possible energy which is achieved in the quench limit. In Fig. 3E, we plot χ over η_{w, m} ∈ [1/32, 32] and η_{w, d} ∈ [1/4, 256], and in SI Appendix, Fig. S1 we show corresponding trajectories of the energy and length for these conditions. Let us imagine fixing all parameters except the spring stiffness and the viscosity. Fig. 3E then illustrates the system’s ability to respond to the active stretch wave under different mechanical and environmental conditions: for small η_{w, m}, corresponding to soft springs, the system stores a large fraction of the possible amount of internal energy (χ ≈ 1), but this fraction decreases as the drag, proportional to η_{w, d}, decreases. By contrast, for large η_{w, m}, the springs closely follow the active stretch wave and never accumulate much internal energy (χ ≈ 0), and the amount that does accumulate decreases as the drag decreases. This energetic perspective thus provides a unified way to understand the dynamic regimes for organisms using Ca²⁺-based motion.

We present a minimal mathematical model of myoneme contraction that captures the kinematic features of the organisms Spirostomum and Vorticella, which exhibit different shapes, sizes, and myoneme arrangements. This model suggests that a fast-traveling wavefront of active stretch, underpinned by rapid Ca²⁺ kinetics, synchronizes contraction across the organism’s length. Using an auxiliary reaction–diffusion model for the Ca²⁺ dynamics, we derive scaling relationships for the speed and width of the active stretch wave that enters the mechanical model. These parameters dictate the rate of chemical driving of the myoneme network. Through nondimensionalization, we determine the minimal independent set of parameters controlling contraction and further examine two limiting cases: a quench limit and a slow wave limit. The latter approximately applies to both Spirostomum and Vorticella in water, as confirmed by fitting the model to experimentally obtained length trajectories. We formally map the dynamics in this limit to the well-studied diffusion equation, with the unitless parameter ζ playing the role of the diffusion constant. The key geometric and mechanical factors that control contractile dynamics include length, stiffness, wave speed, and drag on the organism. These factors manifest through various nondimensional combinations in the different dynamical regimes outlined above. Our analysis reveals three such regimes, each with distinct kinematic signatures such as the symmetry of contraction about the organism’s midpoint and the exponent with which the maximum contraction rate scales with viscosity. We further differentiate these regimes by a quantity representing the peak elastic energy stored in the organism relative to the maximal possible value attained in the quench limit.

It is worthwhile to consider the possible biological implications of the three dynamic regimes. In the quasi-static regime, the myoneme contractile apparatus remains at mechanical equilibrium while the Ca²⁺ wave traverses the organism’s length. In this limit, the contractile speed is fully determined by the wave speed V and does not scale with other parameters. The dissipation of chemical energy is minimized at the expense of the speed of motion. The marine protist Acantharia reportedly has myoneme-based spicules that are used to regulate buoyancy by contracting in a Ca²⁺-powered process that is two orders of magnitude slower than contraction in Spirostomum or Vorticella (19, 29). In this case, contraction may occur quasi-statically. Conversely, the quench regime attains the highest contractile velocities. These high speeds consequently dissipate the most energy because the myoneme apparatus deviates the furthest from mechanical equilibrium. This regime would be suited to processes in which large speeds are physiologically important. For example, Spirostomum excites contraction of nearby organisms through “hydrodynamic trigger waves,” requiring it to move sufficiently rapidly to disturb the fluid, dissipating energy in the process (10). The intermediate regime, where we place the contraction of Vorticella and Spirostomum, is expected to achieve a balance between thermodynamic efficiency and biological functionality. While myoneme-based systems are predominantly found in ciliated protists, similar ATP-independent contractile protein networks exist in a wide range of other organisms across taxa, facilitating a diverse set of biological functions (8–10, 15–27). If similar mechanochemical dynamics apply in those organisms, then various selection pressures may favor a particular regime for each function.

The theory of latch-mediated spring actuation (LaMSA) provides a general description of motion characterized by a fast unlatching event that releases stored elastic energy (68–70). It has been applied to analyze macroscopic biomechanical processes like the striking of a Mantis shrimp’s claw and a human finger snap (71, 72). One may expect that LaMSA also explains myoneme-based contraction, but here energy is stored in the form of a difference in Ca²⁺ chemical potential between the cytosol and the endoplasmic reticulum, and release of a “chemical latch” corresponds to initiation of a traveling wave of Ca²⁺ ions. In addition to this difference in the energy storage medium, LaMSA is designed to apply to macroscopic systems, neglecting the strong dissipative forces which are unavoidable at the cellular level. Furthermore, in myonemes, contraction occurs over an approximately continuous spatial domain of many small elastic elements connected together, compared to the discrete mechanical elements considered in LaMSA. These qualitative features of traditional LaMSA systems compared to the myoneme-based “chemical latch” are illustrated in Fig. 6. Outside the quench limit, the Ca²⁺ wave, which loads mechanical energy into the system, has a finite traversal time. This means that the system can begin relaxing before the energy has been fully loaded. This is another key distinction from traditional LaMSA systems, which are assumed to be fully loaded before relaxing. LaMSA systems also include a constraint force which is lacking in free-swimming organisms like Spirostomum. Exploring the implications of adding constraint forces to myoneme-based systems would be an interesting direction for future research. For example, fixing the ends of Spirostomum and then releasing them after the Ca²⁺ wave has passed would allow artificially pushing the contraction process into the quench regime. Although how physical latches mediate energy flow has only recently been appreciated in biology (69, 73), our model offers a guiding framework for how a chemical latch dynamically controls and tunes force output in ultrafast single-cell organisms.

Finally, our model offers valuable insights for engineering myoneme-based synthetic systems. Researchers have recently demonstrated the feasibility of reconstituting myonemes in vitro and exogenously controlling them using light-patterned Ca²⁺ release (74). Myoneme-based systems, therefore, provide a mechanism for exerting mechanical forces at subcellular scales, with orthogonal chemistry to light-controllable actomyosin and kinesin-microtubule systems (75–79). Our model uncovers several critical parameters that can modulate the dynamics of myoneme contraction. For example, increasing cross-linking would enhance the effective α, elevating the solution viscosity would raise the effective μ, and increasing the effective light-activated Ca²⁺ release rate or reducing Ca²⁺ buffering (thus increasing D_C) would adjust V and W of the active stretch wave. By controlling these experimental parameters, researchers can explore the various regimes identified in this work, paving the way for designing synthetic cells equipped with myoneme-based artificial cytoskeletons that harness the full potential of ultrafast contraction dynamics.

The dependence of the active stretch g(S, T) on the local bound Ca²⁺ concentration may be modeled using the instantaneous linear relation (49)

where B(S, T) is the local concentration of Ca²⁺ bound to myoneme fibers and B_tot is the total concentration of binding sites (which is assumed to be uniform over the spatial domain). Finding g(S, T) then requires specifying the dynamics for the bound Ca²⁺ profile B(S, T). To avoid introducing extra variables, we first analyze the complete dynamics for B(S, T) by extending a tractable model for a traveling Ca²⁺ wave developed by Kupferman et al. ((58)), and we then use results of this analysis to determine an approximate form for g(S, T) to use in the mechanical model.

The Kupferman model for a traveling Ca²⁺ wave is

Here, C(S, T) is the concentration profile of cytosolic Ca²⁺ which diffuses with a diffusion constant D_C, leaks from the cytosolic compartment with a first-order rate constant Γ, and is released into the cytosol by internal stores at a spatially dependent rate JΘ(C − C₀). This last term involves a Heaviside function Θ(C − C₀), which approximates the kinetics of calcium-induced calcium release (63). It “turns on” when the local Ca²⁺ concentration C(S, T) exceeds a spatially uniform threshold C₀, at which point Ca²⁺ is locally released into the cytosol at a rate J. The authors of ref. 58 show that the speed of the traveling wave is given by Eq. 5.

We extend the Kupferman model to allow cytosolic Ca²⁺ to bind to the myoneme assembly. This introduces the variable B(S, T) for the local concentration of myoneme-bound Ca²⁺, and the coupled dynamical equations for C and B are

Here, K₊ is the second-order rate constant for binding of cytosolic Ca²⁺ to open myoneme binding sites, which have concentration B_tot − B. K₋ is the first-order rate constant for the subsequent unbinding of myoneme-bound Ca²⁺.

In SI Appendix, Supplementary Methods, section B, we analyze Eqs. 19 and 20 under mild simplifying assumptions to show that: i) the function B(S, T) is well approximated by a traveling wave with a sigmoidal profile:

ii) its speed V is equal to the speed of diffusing Ca²⁺ wave (cf. Eq. 5), and iii) its width W can be expressed in terms of the chemical parameters according to Eq. 6.

Combing Eqs. 17 and 21 gives the following form for the nonautonomous active stretch:

The additional parameter O is a time that allows an adjustable offset between T = 0 and the time T = O when the center of the wave arrives at S = 0. This allows shifting time to achieve acceptable numerical agreement between boundary and initial conditions, as described in Materials and Methods, Computational methods.

We obtained Spirostomum organisms from Carolina Biological and cultured them at room temperature. We prepared culturing medium by mixing one part boiled hay medium solution (Ward’s Science, 470177-390), four parts boiled spring water, and two boiled wheat grains (Carolina, Item #132425) as described previously (36). We inoculated organisms into fresh media approximately twice per week. For imaging, we pipetted several organisms at a time into microscopy slides with multiple channels (Ibidi, μ-Slide VI 0.4) or with multiple wells (Ibidi, μ-Slide 8-well). The channels are 400 μm high, allowing the organisms to swim and contract without confinement.

Our Vorticella samples were collected from Morse Pond near the Marine Biological Laboratory. We dissected small pieces of plant material to which individual or small groups of Vorticella were attached. We added these pieces of plant matter, along with enough imaging medium to cover them, into 8-well coverslip-bottomed sample chambers (Ibidi, μ-Slide 8-well).

To alter the viscosity, we prepared stock solutions of Ficoll PM 400 (Sigma F4375) in water and diluted them to concentrations between 3% and 12% (weight per volume). In the case of Spirostomum, we pipetted organisms directly into the Ficoll solution. In the case of Vorticella, we rinsed the samples several times with the solution within the sample chamber to ensure that the appropriate viscosity was reached. In both cases, slow pipetting was important to prevent shearing of the organisms in the high-viscosity solutions. To calculate viscosity ν as a function of percentage p of Ficoll in water, we used the calibration function ν = e^0.162p reported in refs. 80 and 81.

We performed brightfield imaging on either a Zeiss inverted microscope equipped with a 0.55 NA condenser, a 0.63× defocusing lens, DIC optics, and a 10× air objective or a Nikon Ti2 equipped with a 2× (MRD70020) or 10× (MRD70170) air objective. In both cases, we collected images using a pco.dimax CS camera equipped with at least 2,000 and up to 10,000 Hz.

For quantifying Spirostomum images, we first used Fiji (82) to threshold the images to identify the pixels associated with the organisms, and then we used Matlab to run Morphometrics (83) to obtain the cell contours. We used the pill-mesh option of Morphometrics to obtain cell shape and length. See SI Appendix, Fig. S3 and Movie S3 for a visualization of this mesh generation procedure.

For quantifying Vorticella length over time, we manually tracked the position of the end of the stalk attached to the cell body using home-built Python software.

All numerical integration was performed with Mathematica’s NDSolve function, using the method of lines option with a spatial discretization based on tensor product grids. The accuracy and precision goals were set high enough to ensure reasonable numerical convergence. To achieve agreement between the initial conditions Eq. 4 and the boundary conditions Eq. 3, we increased the parameter O in Eq. 22 so that g(S, 0)≈1 across the computational domain. To solve the quench model Eq. 15, we evaluate the exact solution given in SI Appendix, Eq. 60 using 50 terms in the Fourier expansion.

C.F., A.T.M., X.L., J.E.H., F.C., M.W.E., S.V., A.R.D., and M.S.B. designed research; C.F., A.T.M., M.W.E., S.V., A.R.D., and M.S.B. performed research; C.F., A.T.M., M.W.E., S.V., A.R.D., and M.S.B. analyzed data; and C.F., A.T.M., M.W.E., S.V., A.R.D., and M.S.B. wrote the paper.

The authors declare no competing interest.