Nonequilibrium dynamics of probe filaments in actin-myosin networks

TL;DR

This paper investigates how active forces from molecular motors affect the non-equilibrium bending dynamics of probe filaments in actin-myosin networks, revealing broken detailed balance and correlations among bending modes.

Contribution

It provides an analytical framework linking motor-induced forces to non-equilibrium mode correlations and probability currents in probe filaments within active cytoskeletal networks.

Findings

Motor forces induce measurable non-equilibrium bending dynamics.

Active forces create correlations among bending eigenmodes.

Probability currents in mode space break time-reversal symmetry.

Abstract

Active dynamic processes of cells are largely driven by the cytoskeleton, a complex and adaptable semiflexible polymer network, motorized by mechanoenzymes. Small dimensions, confined geome- tries and hierarchical structures make it challenging to probe dynamics and mechanical response of such networks. Embedded semiflexible probe polymers can serve as non-perturbing multi-scale probes to detect force distributions in active polymer networks. We show here that motor-induced forces transmitted to the probe polymers are reflected in non-equilibrium bending dynamics, which we analyze in terms of spatial eigenmodes of an elastic beam. We demonstrate how these active forces induce correlations among these mode amplitudes, which furthermore break time-reversal symmetry. This leads to a breaking of detailed balance in this mode space. We derive analytical predictions for the magnitude of resulting probability currents in mode space in the white-noise limit of motor activity. We relate the structure of these currents to the spatial profile of motor- induced forces along the probe polymers and provide a general relation for observable currents on two-dimensional hyperplanes.

Equations

$$\underset{-\infty }{\overset{t}{\int }}\mathrm{d}{t}^{\prime }\alpha (t-t^{\prime })r_{\perp }(t^{\prime })$$

$$\hat{\alpha }(\omega )$$

$$⟨\hat{\xi }(s,\omega )\overline{\hat{\xi }}(s^{\prime },\omega )⟩=\frac{2k_{B}T}{\omega }\delta (s-s^{\prime })\mathrm{Im}\left[\hat{\alpha }(\omega )\right].$$

$$f_{\text{M}}(s,t)=\sum _{n=1}^{N_{\text{M}}}{f}_{n}g\left(s-s_n\right){\mathcal{T}}_n(t)$$

$$⟨{\mathcal{T}}_n(t){\mathcal{T}}_m(t^{\prime })⟩=C_1+C_2{\delta }_{n,m}{e}^{-\frac{|t-t^{\prime }|}{{\tau }_{\text{M}}}}$$

$$\underset{-\infty }{\overset{t}{\int }}\mathrm{d}{t}^{\prime }\alpha (t-t^{\prime })a_q(t^{\prime })$$

$$f_{\text{M},q}(t)$$

$$=\sum _{n=1}^{N_{\text{M}}}{f}_n{y}_q(s_n){\mathcal{T}}_n(t).$$

$${\hat{a}}_q\left(\omega \right)$$

$$⟨|{\hat{a}}_q(\omega ){|}^2{⟩}^{\text{Th}}$$

$$\hat{G}=G_0+i\eta \omega$$

$$⟨a_q(t)a_w(t^{\prime })⟩=⟨a_q(t)a_w(t^{\prime }){⟩}^{\text{Th}}+⟨a_q(t)a_w(t^{\prime }){⟩}^{\text{M}}$$

$$⟨a_q(t)a_w(t^{\prime }){⟩}^{\text{Th}}$$

$$⟨a_q(t)a_w(t^{\prime }){⟩}^{\text{M}}$$

$${\mathbf{F}}_{q,w}$$

$$⟨a_q(t)⟩⟨a_w(t)⟩$$

$${\mathcal{C}}_{q,w}(t-t^{\prime })$$

$$-2\frac{{\tau }_q}{{\tau }_{\text{M}}}\frac{e^{-\frac{|t^{\prime }-t|}{{\tau }_q}}}{\left(1-{\left(\frac{{\tau }_q}{{\tau }_{\mathrm{M}}}\right)}^2\right)\left(1+\frac{{\tau }_w}{{\tau }_q}\right)})$$

$${\mathcal{C}}_{q,q}(t-t^{\prime })$$

$$⟨f_{\text{M},q}(t)f_{\text{M},w}(t^{\prime })⟩$$

$$=\frac{1}{L^2}\sum _{i,j}{f}_i{f}_j{y}_q(s_i)y_w(s_j){\delta }_{i,j}{C}_2{e}^{-\frac{|t-t^{\prime }|}{{\tau }_{\text{M}}}}$$

$$=\frac{C_2{e}^{-\frac{|t-t^{\prime }|}{{\tau }_{\text{M}}}}}{L^2}\sum _i{f}_i^2{y}_q(s_i)y_w(s_i)$$

$$=\frac{C_2{e}^{-\frac{|t-t^{\prime }|}{{\tau }_{\text{M}}}}}{L^2}{\mathbf{F}}_{q,w}.$$

$${\mathbf{F}}_{q,w}\stackrel{{\ell }_{\mathrm{M}}\to 0}{\to }\frac{1}{{\ell }_{\mathrm{M}}}\underset{0}{\overset{L}{\int }}\mathrm{d}s{f}^2(s)y_q(s)y_w(s).$$

$$⟨{\psi }_q(t){\psi }_w(t^{\prime })⟩=$$

$$\gamma \frac{\mathrm{d}{a}_q}{\mathrm{d}t}(t)$$

$${\mathbf{D}}_{q,w}$$

$$\frac{\partial \rho }{\partial t}\left(\vec{a},t\right)$$

$${\mathbf{C}}_{q,w}=⟨a_q(t)a_w(t){⟩}^{\psi },$$

$${\mathbf{C}}_{q,w}$$

$$\vec{j}(\vec{a})=(\mathbf{K}+\mathbf{D}{\mathbf{C}}^{-1})\vec{a}\rho (\vec{a})\equiv \mathbf{\Omega }\vec{a}\rho (\vec{a}),$$

$${\tilde{j}}_q^{(2)}({\tilde{a}}_q,{\tilde{a}}_w)$$

$$=\underset{-\infty }{\overset{\infty }{\int }}\mathrm{d}{\tilde{a}}_{q_1}\cdots \underset{-\infty }{\overset{\infty }{\int }}\mathrm{d}{\tilde{a}}_{q_N}\sum _v{\tilde{\mathbf{\Omega }}}_{q,v}{\tilde{a}}_v\tilde{\rho }(\vec{\tilde{a}})$$

$$={\tilde{\mathbf{\Omega }}}_{q,w}{\tilde{a}}_w(2\pi {)}^{-1}{e}^{-\frac{{\tilde{a}}_q^2+{\tilde{a}}_w^2}{2}}$$

Full text

The emergent field of active matter aims to develop systematic descriptions of stochastic out-of-equilibrium phenomena in energy-dissipating soft matter systems MacKintosh and Schmidt (2010); Marchetti et al. (2013); Prost et al. (2015). A prominent motivation for such studies are the dynamics observed in living cells and tissues. In vitro model systems based on one of the main ingredients of the cellular cytoskeleton, filamentous actin, have been playing a prominent role in pioneering experimental studies of active matter. Building on earlier microrheology experiments on equilibrium actin networks Ziemann et al. (1994); Schnurr et al. (1997); Gittes et al. (1997); Hinner et al. (1998); Gardel et al. (2003); Koenderink et al. (2006); MacKintosh and Schmidt (2010), similar experiments in reconstituted actin networks including myosin motor proteins have revealed that motor activity can drastically alter the mechanical response of actin networks Lau et al. (2003); Liu et al. (2006); Koenderink et al. (2009); Broedersz and MacKintosh (2011); Ronceray et al. (2016) and significantly enhance fluctuations Brangwynne et al. (2008); Mizuno et al. (2007, 2009). Actin-myosin model systems have also been used to study structural self-organization and pattern formation Soares e Silva et al. (2011); Murrell and Gardel (2012); Alvarado et al. (2013); Lenz (2014); Weber et al. (2010); Schaller et al. (2011).

The standard approaches to characterizing cellular dynamics involve direct microscopic imaging of the motion of whole cells or identifiable features within cells, or the monitoring of displacements of probes inserted into cells or attached to the outside of cells. Small probes are often tracked via fluorescence microscopy. Micron-sized colloidal beads attached to cells or injected into cells have been used to track fluctuations Fabry et al. (2001); Hoffman et al. (2006); Deng et al. (2006); Mizuno et al. (2009); Guo et al. (2014). However, these probes clearly lack spatial resolution and often cannot enter confined geometries such as the actin cortex of non-adherent cells. An interesting alternative is the direct fluorescent labeling of parts of the cytoskeleton itself, for example microtubules Brangwynne et al. (2008) or the insertion of high-aspect-ratio filamentous probes, such as fluorescent single-walled carbon nanotubes Fakhri et al. (2010, 2014); Tan et al. (2016) that can penetrate into tight spaces in the cell.

Once the stochastic dynamics are tracked in an active material, it can be challenging to distinguish equilibrium from non-equilibrium stochastic motions and to quantify the extent of activity. Colloids in active networks typically exhibit non-Gaussian displacement distributions Toyota et al. (2011); Ben-Isaac et al. (2015); Turlier et al. (2016), but that feature by itself does not prove non-equilibrium. Scaling regimes of non-equilibrium filament fluctuations have been theoretically investigated and are expected to deviate from equilibrium predictions Everaers et al. (1999). If fluctuations can be compared to material response properties, the fluctuation-dissipation theorem can be applied to quantify activity  Hoffman et al. (2006); Mizuno et al. (2007, 2008). Recently, a non-invasive method was introduced to discern active fluctuations based on the violation of detailed balance  Battle et al. (2016); Gladrow et al. (2016). This method can be applied on standard microscopic imaging data, and is based on analyzing probability flux patterns in phase spaces constructed from two or more degrees of freedom of the system, which are either directly coupled or share their driving forces.

In contrast to point-like probe particles or spherical beads, extended filaments offer an easily accessible spectrum of simultaneously observable variables: their bending modes Aragon and Pecora (1985); Gittes (1993); Brangwynne et al. (2007). We can decompose the instantaneous filament conformation into a sum of dynamic bending eigenmodes with the following properties: (i) Particular spatial modes act as reporters of dynamics at their respective characteristic length and time scale; since viscous relaxation times are length dependent, modes also have characteristic relaxational timescales. (ii) In equilibrium, dynamic eigenmodes are statistically independent and each mode amplitude individually fulfils a fluctuation-dissipation theorem. In a motor-activated network, by contrast, a filament will receive random kicks from the network generated by myosin motors and transmitted by the network. Such motor activity results in enhanced mode amplitude fluctuations and may, as we will show, abrogate the independence of the eigenmodes. In a preceding letter Gladrow et al. (2016), we have shown that resultant mode cross-correlations are accompanied by circular patterns of probability fluxes in mode space, which indicate a breakdown of detailed balance. In this way, the analysis of bending eigenmodes of embedded probe filaments can be used to indicate non-equilibrium dynamics in an active polymer network.

In this paper, we develop the analysis of filament motion in an active viscoelastic medium in more detail and specify how the spatial structure of actively induced network-probe interactions, as sketched in Fig. 1, translates into mode cross-correlations. We discuss the breaking of Onsager’s time reversal symmetry in this system and show how cross-correlations can provide guidance to identify motor-induced dynamics. Furthermore, we also discuss an additional non-equilibrium marker in the white-noise limit, the frequency associated with the circulatory probability currents that arise in non-equilibrium steady-states.

I The filament model

We consider the non-equilibrium dynamics of a semiflexible filament, which is embedded in a polymer network with a mesh size smaller than the probe filament’s contour length. The network is actuated by molecular motors that act as homogeneously dispersed contractile force dipoles. While we approximate the viscoelastic network as a continuum, we assume motor forces to act on the filament at discrete points, where the probe filament is assumed to be coupled mechanically to the meshwork. Without loss of generality, both the meshwork and the filament are described in two dimensions. The filament itself is modeled as an inextensible worm-like chain Kratky, O. and Porod (1949); Aragon and Pecora (1985); Broedersz and MacKintosh (2014). We therefore treat its shape as a continuously differentiable curve in space $\vec{r}(s,t)$, parametrized over time $t$ and its arclength $s$ as sketched in Fig. 2.

For every given point in time, we decompose the shape of the probe filament into a transverse and a parallel part with respect to the end-to-end vector $\vec{R}(t)$. This separation yields two internal variables, the transverse $r_{\perp}(s,t)$ and parallel $r_{\parallel}(s,t)$ relative position, which are connected by local inextensibility $\left\lVert\frac{\partial\vec{r}}{\partial s}\right\rVert=1$. The equation of motion governing both variables can be retrieved from a variation of the worm-like chain Hamiltonian $\mathcal{H}=\kappa/2\int_{0}^{L}\mathrm{d}s\,\partial^{2}\vec{r}/\partial s^{2}$, with bending rigidity $\kappa$ and $L$ denoting the total arc length of the filament. We here focus on the linear dynamics of transverse deviations at a given arc length $s$ along the polymer, described by

[TABLE]

where the usual stochastic term $\xi(s,t)$ models the net thermal force exerted by the surroundings of the filament consisting of polymer network and solvent. The viscoelastic kernel $\alpha$ in Eq. (1) is related to the bulk shear modulus $G$ via the generalized Stokes theorem Gittes et al. (1997)

[TABLE]

where $k_{0}\approx 4\pi/\ln(L/d)$ is a geometric factor, which also appears in the transverse drag coefficient $\gamma$ of a infinitesimal rod segment $\gamma\approx k_{0}\eta$ Howard (2001) in a medium with viscosity $\eta$.

The thermal force has zero mean $\langle\xi(s,t)\rangle=0$ and a temporal power spectrum satisfying Lau et al. (2003)

[TABLE]

In order to account for the impact of motor forces produced by myosin motors in the vicinity of the probe filament, we include an additional force term $f_{\text{M}}(s,t)$ in Eq. (1) Everaers et al. (1999). In a steady-state scenario, motor protein action gives rise to a fluctuating, but stationary profile of non-thermal forces along the backbone of the probe. We assume the probe filament to be stationary and to not reptate through the network. This can be achieved, for instance, by linking the probe filament at some points to the network. Reptational movement of the probe filament would cause sampling from a changing spatial motor force profile (see Fig. 1), and thus lead to a blurring of the non-equilibrium dynamics we seek to describe. The structure of a crosslinked network is characterized by its mesh-size, which determines the density $\ell_{\text{M}}$ of motor-probe interaction points $s_{n}$. We model the motor-induced force as

[TABLE]

with $N_{\text{M}}$ denoting the number of entanglement points where motors affect the probe filament, $f_{n}$ denoting the impact of motor $n$ and $g(s)$ a general spatial kernel describing how a motor impacts on the filament. For simplicity, we choose here a point-like spatial kernel, i.e. $g(s)=\delta(s)$, consistent with experiments that suggest a rather narrow force profile Brangwynne et al. (2008). The temporal profile of force generated by an individual motor $\mathcal{T}(t)$ in Eq. (4) is described as a telegraph process.

Non-muscle myosin motors act as oligomeric complexes and produce forces between actin filaments with correlation times of order 10 s Fakhri et al. (2014); Guo et al. (2014). This time scale is well separated from temporal correlations of thermal noise, but is on the order of the relaxation times of bending modes that play the main role here. To simplify our theoretical description, we assume motors to instantaneously develop maximal force with characteristic switching rates between the off and the on state $1/\tau_{\text{on}}$ and vice-versa $1/\tau_{\rm off}$, as proposed in MacKintosh and Levine (2008). A telegraph process $\mathcal{T}(t)$ switches between zero and one without memory and can therefore be considered a simple model of a molecular motor.

The activity of different motors is assumed to be uncorrelated, while the autocorrelation of a given motor telegraph process is exponential in lag time, such that we obtain Gardiner (2009)

[TABLE]

with dimensionless constants $C_{1}=\tau_{\rm off}^{2}/\left(\tau_{\rm on}+\tau_{\rm off}\right)^{2}$ and $C_{2}=\tau_{\rm on}\tau_{\rm off}/\left(\tau_{\rm on}+\tau_{\rm off}\right)^{2}$. In addition, we defined a motor timescale $\tau_{\text{M}}^{-1}=\tau_{\text{on}}^{-1}+\tau_{\text{off}}^{-1}$. A Fourier transformation of Eq. (5) yields a Lorentzian power spectrum in accord with literature MacKintosh and Levine (2008). In other words, for high frequencies the motor force spectrum follows a power law $S_{\mathcal{T}}(\omega)\sim\omega^{2}$, whereas for lower frequencies, the power spectrum becomes essentially white-noise-like $S_{\mathcal{T}}\sim\text{const}$.

II Characterization of non-equilibrium mode dynamics

We start to analyze filament mode dynamics by expanding $r_{\perp}(s,t)$ in orthogonal eigenmodes $y_{q}(s)$ of the beam operator in Eq. (1) as $r_{\perp}(s,t)=L\sum_{q}a_{q}(t)y_{q}(s)$. This expression implies a choice of units: spatial modes $y_{q}(s)$ and mode amplitudes $a_{q}(t)$ have dimensions of $1/{\rm length}^{1/2}$ and ${\rm length}^{1/2}$ respectively. Projecting Eq. (1) onto a particular spatial mode $y_{q}(s)$ leads to the equation of motion in mode space

[TABLE]

where we made the implicit assumption, that motor and thermal forces do not affect each other. Indexed quantities denote projected variables, such as the projected thermal noise $\xi_{q}(t)=L^{-1}\int_{0}^{L}\,\mathrm{d}s^{\prime}\xi(s^{\prime},t)y_{q}(s^{\prime})$. We note that thermal forces of different modes do not correlate $\langle\xi_{q}\xi_{w}\rangle\propto\delta_{q,w}$ due to the orthogonality of the spatial modes $y_{q}(s)$. Indices $q$ and $w$ refer to the corresponding wave vectors, which are usually discretized. For instance, for the relevant case of free-end boundary conditions, $q$ can assume values $q(n)\approx(n+\frac{1}{2})\pi/L$ for any integer $n$. Furthermore, the projected motor-induced force can be written as

[TABLE]

Experimentally accessible quantities, such as mode correlations may now be computed with the projected quantities. To this end, we Fourier transform equation (6) in time and obtain

[TABLE]

with the mode response function $\hat{\chi}_{q}\left(\omega\right)~{}=~{}\left(\hat{\alpha}\left(\omega\right)+\kappa q^{4}\right)^{-1}$. Although motor-induced forces may surpass thermal forces by orders of magnitude, we choose to keep both terms in Eq. (8). This allows us to smoothly transition between purely thermal and purely active dynamics. In the thermal case, the above equation together with the assumed thermal noise spectrum Eq. (3) yields the Fluctuation-Dissipation theorem of mode $q$

[TABLE]

As a first application of our model, we characterize the physiologically relevant case of motor-induced fluctuations of relatively stiff filaments in a viscoelastic medium, such as microtubules in actin-myosin networks Brangwynne et al. (2008); Fakhri et al. (2014). In particular, we calculate the deviation from the equilibrium mode variance shown in Eq. (9). The viscoelastic response of crosslinked networks of semiflexible filaments typically shows two regimes of the complex shear modulus, a high-frequency regime where $\hat{G}(\omega)\sim\omega^{3/4}$ Gittes and MacKintosh (1998) and a low-frequency plateau regime Maggs (1997); Koenderink et al. (2006); Isambert and Maggs (1996); Morse (1998). The cytoplasm of living cells exhibits an elastic shear modulus with a weak frequency dependence Fabry et al. (2001); Fakhri et al. (2014) up to  $100$ Hz. The frequency regime we need to consider here is not only determined by the relaxation times of our probe filament for the modes we can resolve, but also by the timescale of motor activity. Experiments suggest, that motor-induced fluctuations become negligible for frequencies higher than $100$ Hz Mizuno et al. (2007, 2009). The processivity time of myosin in cells was measured to be about $5$ s Fakhri et al. (2014). In in vitro model systems, values for $\tau_{\text{off}}$ and $\tau_{\text{on}}$ depend on salt and ATP concentrations. Given these facts, we focus on the low-frequency regime and assume, for simplicity, a plateau-like shear modulus $\hat{G}(\omega)$ given by

[TABLE]

with an elastic modulus $G_{0}$ and viscosity $\eta$. We can now explicitly calculate correlation functions by applying the Wiener-Khinchin theorem to $\langle\hat{a}_{q}\overline{\hat{a}}_{w}\rangle$. This yields the mode correlation function, which decomposes into the usual thermal and an additional motor-induced part

[TABLE]

respectively given by

[TABLE]

We have introduced here the mode relaxation time $\tau_{q}=1/(\kappa/\gamma q^{4}+G_{0}/\eta)$ and a coupling matrix $\mathbf{F}_{q,w}$. The elements in this coupling matrix encode the spatial structure of motor-induced forces in terms of bending modes,

[TABLE]

In contrast to the purely thermal case, Eq. (13) shows that in the active case, modes are no longer independent, but correlate with a magnitude determined by the coupling matrix $\mathbf{F}_{q,w}$ and the function $\mathcal{C}_{q,w}(\Delta t)$ that describes the temporal evolution of the motor-induced mode correlation function.

[TABLE]

For simplicity, we assume the average active forces and torques exerted on the filament to vanish, requiring $\langle f_{\rm M}(s,t)\rangle=0$ for all $s$. This leads to $\sum_{n}f_{n}y_{q}(s_{n})=0$ and thus $\langle a_{q}(t)\rangle\langle a_{w}(t)\rangle=0$.

The function $\mathcal{C}_{q,w}(t-t^{\prime})$ in Eq. (13) is limited at short times by the motor decorrelation time, is independent of the precise spatial arrangement of the motor-induced forces, and is given by

[TABLE]

where $\tau_{q}$ must be the relaxation time associated with time coordinate $t$, where $t>t^{\prime}$, and $\tau_{w}$ is associated with the mode evaluated at $t^{\prime}$. While for zero lag time, i.e. $t=t^{\prime}$, the correlator is symmetric, indices can not be exchanged for $t\neq t^{\prime}$. Thus, Onsager’s time reversal symmetry is broken in this system, as demonstrated in Fig. 3 Onsager (1931a, b). The simplest behavior of $\mathcal{C}_{q,w}(t-t^{\prime})$ can be observed for modes with relaxational time scales $\tau_{q}$, $\tau_{w}$ that are much longer than the motor time scale $\tau_{\text{M}}$. In this limit, $\mathcal{C}_{q,w}(t-t^{\prime})$ evolves as $\text{Exp}(-\Delta t/\tau_{q})$. In Fig. 3 this transition of the cross-correlation to simple exponential behavior is exemplified for the mode pair $(q,w)=(2,4)$ for three different motor time scales. By contrast, in equilibrium, these cross-correlations vanish for all lag times $\Delta t$. It is interesting to note that one could therefore directly infer non-equilibrium dynamics from asymmetries in the mode correlation functions. The breaking of the Onsager time reversal symmetry of the mode correlation function is equivalent with the breaking of detailed balance due to motor activity in our model.

As discussed previously Gladrow et al. (2016), in the diagonal case, $\mathcal{C}_{q,q}(t-t^{\prime})$ can be written in a more compact form,

[TABLE]

We now turn to the coupling matrix $\mathbf{F}_{q,w}$ defined in Eq. (14). Our motivation for defining these coefficients in this particular way can best be illustrated by regarding the active force cross-correlations $\langle f_{q}(t)f_{w}(t^{\prime})\rangle$ using Eqs. (5) and (4). Under the assumption that all motor-force processes have the same time scale and interact with the probe in only one location each, we find

[TABLE]

Thus the elements of the coupling matrix $\mathbf{F}_{q,w}$ appear in the above equation as the only index-dependent quantities.

The coupling matrix for a given fluctuating filament will depend on the spatial distribution of points through which forces are transmitted to the filament. Fig. 4 shows three examples of sets of interaction-points $\{s_{n}\}$ with the associated coupling matrices. The first two are disordered distributions, while the third and fourth example constitute periodic patterns. Since the coupling matrix depends on both the choice of motor-probe interaction strengths $f_{n}$ and the placement of interaction points $s_{n}$, one could consider the ensemble average over all possible disordered network-probe interactions $\langle\mathbf{F}_{q,w}\rangle_{\text{dis}}$. However, under the assumption of a uniform distribution of possible interaction points $\rho(s_{n})=1/L$, we find $\langle\mathbf{F}_{q,w}\rangle_{\text{dis}}\propto\delta_{q,w}$ due to the orthogonality of the dynamic normal modes. Consequently, the mode correlations will disappear in this ensemble average in spite of non-equilibrium driving. This result holds true, even when the distances between the interaction points $\Delta s_{n}=s_{n+1}-s_{n}$ are drawn from an exponential distribution with some characteristic spatial frequency $1/\ell_{\rm M}$. The ensemble average of the coupling matrix $\langle\mathbf{F}_{q,w}\rangle_{\text{dis}}$ over disordered network configurations with such an exponential distribution of motor interaction points is still diagonal, despite the structured probe-mesh interaction. This is because the relative position of the probe to the network is a uniformly distributed variable and hence the average force geometry is uniform.

In any case, diagonality of $\mathbf{F}_{q,w}$ implies that both broken Onsager’s time reversible symmetry and breaking of detailed balance (see below) will not be visible in normal-mode phase space. In other words, ensemble averaging over trajectories recorded from different probe filaments and thus over different force profiles will tend to conceal the non-equilibrium nature of the dynamics. The same will be true in the limit of high spatial densities of interaction points $\ell_{\rm M}\to 0$, where the sum in Eq. (14) approaches an integral over $s$,

[TABLE]

If force amplitudes don’t spatially vary, the coupling matrix $\mathbf{F}_{q,w}$ therefore also becomes diagonal in the high-density limit Gladrow et al. (2016).

III Brownian dynamics simulations

To numerically test our analytical predictions, we carried out Brownian dynamics simulation of Eq. (1). To this end, we approximated the probe filament as a discrete worm-like chain composed of $N$ beads, connected by stiff elastic springs to ensure local inextensibility.  We discretized the fourth-derivative in equation (1) with a central-stencil scheme. Free-end boundary conditions were implemented by endowing the chain with four ghost-beads, two at each end respectively, to which the fourth-derivative stencil was applied. At each step, we positioned the ghost beads along the tangent at the respective end of the filament. The elastic confinement of the probe filament was simulated by springs attached to every bead of the chain and to a fixed point at a unit distance extending orthogonally to the local filament tangent. Motor-interaction points $s_{n}$ were positioned at an equal distance $\ell_{\text{M}}$ along the filament to apply the assumptions made in the previous section. Motor action was modeled by a telegraph process with a variance as described by Eq. (5). In order to obtain mode amplitude traces $\vec{a}(t)$ from the simulated dynamics, we projected the filament backbone onto the free-rod eigenmodes of the beam equation $y_{q}(s)$ Brangwynne et al. (2007) at each time point $t_{i}$.

IV Active mode dynamics violate detailed balance

Non-equilibrium driving of a network by embedded motors results in an enhancement of fluctuations of an inserted probe filament. We showed above that dynamic normal modes, which are uncorrelated in equilibrium, may start to correlate in such an active environment (Eqs. (12) and (13)). Cross-correlations alone, however, are not a reliable tool to diagnose non-equilibrium or to infer force profiles. Especially lower mode correlations converge rather slowly and tensile effects, which we neglected in our treatment, may result in non-zero cross-correlations in experiments, even in equilibrium.

A more robust measure of non-equilibrium is the breaking of detailed balance. We analyze the breaking of detailed balance by looking for circulating probability currents in normal mode phase space. In particular, we discuss how the structure of such currents depends on the geometry of interactions between our probe filament and the motor-activated network.

While thermal forces drive all modes equally and independently (Eq. (12)), the network-probe interactions $f_{\text{M}}(s,t)$, which are assumed to be distributed heterogenously, excite modes unevenly and in a correlated fashion. In conjunction with different mode relaxation times, this gives rise to a directed probability current in the dynamic normal mode space of the probe filament. We can understand this phenomenon quantitatively in the white-noise limit of the mode equation of motion (see Eq. (6)).

On long timescales $T\gg\tau_{\text{M}}$, the internal motor time scale becomes negligible, and we can approximate the motor forces to appear as white-noise processes. In this limit, we rewrite the correlator of the telegraph process $\langle\mathcal{T}_{n}(t)\mathcal{T}_{m}(t^{\prime})\rangle\to C_{1}+\frac{C_{2}}{2}\delta_{n,m}\tau_{\text{M}}\delta(t-t^{\prime})$. The motor time scale $\tau_{\text{M}}$ here remains only as a proportionality constant of the variance of the active forces to ensure appropriate dimensions. We can now absorb the motor processes $f_{\text{M},q}(t)$ into the thermal force $\xi_{q}(t)$, which yields a new white noise process $\psi_{q}(t)$ with a correlator

[TABLE]

Furthermore, we again adopt the Kelvin-Voigt description of embedding network and solvent $\hat{G}(\omega)=G_{0}+i\omega\eta$, since we now operate in a regime of even lower frequencies. Such a complex shear modulus with constant real part corresponds to a continuum of springs and simple drag penalizing transverse deviations. We can therefore split up the left-hand side of Eq. (1) into an elastic and viscous force. The Langevin equation of motion for the normal modes Eq. (6) then becomes

[TABLE]

The corresponding diffusion matrix is

[TABLE]

with this result, we are in a position to write down the Fokker-Planck equation of the system

[TABLE]

where $\mathbf{K}_{q,w}=-1/\tau_{q}\delta_{q,w}$ denotes the deterministic matrix. The steady-state solution of the above equation is a Gaussian probability density $\rho(\vec{a})=\mathcal{N}^{-1}e^{-\frac{1}{2}\vec{a}^{T}\mathbf{C}^{-1}\mathbf{a}}$ with a normalization constant $\mathcal{N}^{-1}$ and a correlation matrix defined by

[TABLE]

where the superscript indicates that the average is taken in the white-noise limit over $\psi(t)$ defined in Eq. (19). We can analytically obtain the white-noise correlation matrix $\mathbf{C}$ either from Eq. (22) or directly from a power expansion of Eq. (13) up to linear order in $\tau_{\text{M}}$

[TABLE]

For a random collection of motor interaction points $\{s_{n}\}$, drawn from a constant probability density over the filament length $L$, the coupling matrix and thus the correlation matrix may exhibit off-diagonal elements as exemplified in Fig. 6.

Since the correlation matrix $\mathbf{C}_{q,w}$ in Eq. (24) is proportional to the coupling matrix $\mathbf{F}_{q,w}$, the relative values of the entries of this matrix can, in principle, be inferred from the correlation, if mode relaxation times $\tau_{j}$ are known. A direct comparison of Fig. 4 and Fig. 6 reveals how the motor force distribution $f_{\text{M}}(s,t)$ translates into mode cross-correlations.

The quantity in brackets in Eq. (22) can be identified as the probability current $\vec{j}(\vec{a})$, which in the steady-state limit becomes

[TABLE]

with $\mathbf{\Omega}$ denoting a matrix of frequencies, which we discuss below. By inserting the result for $\mathbf{C}$ into Eq. (24), we can directly obtain the current $\vec{j}(\vec{a})$.

The probability currents arising from regularly spaced interaction points $s_{n}$ are plotted in Fig. 7. We calculated these currents in two-dimensional reduced mode spaces $a_{q}\times a_{w}$, where dynamics of other modes were disregarded for simplicity. Due to the symmetry of the motor induced forces $f_{\text{M}}(s)$ in Fig. 5 and Fig. 7, currents develop only in pairs of mode numbers of the same parity, for instance odd-odd, as one would expect from the structure of the coupling matrix (see third panel in Fig. 4).

In our simulations, we measured cumulative probability currents $\vec{J}_{q,w}$ over small rectangles $\square_{n,m}$ in the corresponding $a_{q}\times a_{w}$ subspace. More precisely, we applied the current estimator proposed in Battle et al. (2016) to all simulated two-mode amplitude trajectories $\{a_{q}(t_{i}),a_{w}(t_{i})\}_{i}$. In Fig. 5, coarse-grained probability currents inferred from Brownian dynamics simulations are shown. Following the bootstrapping technique described previously Battle et al. (2016), we calculated error estimates (white ellipses) in order to distinguish significant currents.

V Cycling Frequencies in the White-Noise Approximation

If a system violates detailed balance it does so in any coordinate system; it should not depend on the choice of coordinates whether it is in equilibrium or not. In fact, detailed balance requires the product $\mathbf{K}\mathbf{D}$ to be symmetric Weiss (2003), which is not fulfilled for general $\mathbf{F}_{q,w}$. In a steady-state, the current $\vec{j}(\vec{a})$ must necessarily be purely rotational to fulfill $\vec{\nabla}\cdot\vec{j}(\vec{a})=0$. A rotational probability current manifests itself in a preferred angular direction of the underlying stochastic dynamics of modes $\vec{a}(t)$. Trajectories sampled from such a state must therefore exhibit, on average, a cycling motion. The characteristic frequencies $\omega_{q,w}$ of this cycling motion are, in principle, experimentally accessible and can be considered a scalar measure of irreversible dynamics.

In a two-dimensional space, the stochastic dynamics of the system on average cycles at a single characteristic frequency, which can be calculated analytically as the imaginary part of the first of the two nonzero eigenvalues of the frequency matrix $\mathbf{\Omega}$ Gladrow et al. (2016); Weiss (2003); Ghanta et al. (2017). In higher-dimensional systems, however, the situation is less clear. Without prior knowledge about the structure of the motor-generated forces, we can transform the observed mode traces $\vec{a}(t)$ into so-called correlation-identity coordinates $\widetilde{a}_{q}(t)$, where $\langle\widetilde{a}_{q}(t)\widetilde{a}_{w}(t)\rangle=\widetilde{\mathbf{C}}_{q,w}=\delta_{q,w}$. This is achieved by multiplying the mode vector $\vec{a}(t)$ by the matrix square-root of the inverse-correlation matrix $\sqrt{\mathbf{C}^{-1}}$. We note that the wave-vector indices $q,w,...$ become dimensionless in correlation-identity coordinates.

In this coordinate system, theoretically meaningful cycling frequencies can be estimated by projecting the system onto two-dimensional hyperplanes, disecting mode space. If we followed only $\widetilde{a}_{q}(t)$ and $\tilde{a}_{w}(t)$, we would perceive an apparent frequency $\omega_{q,w}$, which, in general, would be different from any eigenvalue of the frequency matrix $\tilde{\mathbf{\Omega}}$, but can nevertheless be calculated. Indeed, as shown below, the apparent cycling frequency in this hyper-plane spanned by correlation-identity coordinates is given by the matrix elements of $\tilde{\mathbf{\Omega}}_{q,w}$, which here read $(\widetilde{\mathbf{K}}+\widetilde{\mathbf{D}})_{q,w}$.

We begin by showing that integrating out unobserved variables in correlation-identity coordinates effectively reduces the matrix $\tilde{\mathbf{\Omega}}$ to those elements $\tilde{\mathbf{\Omega}}_{q,w}$ corresponding to the two observed degrees of freedom. We indicate $\widetilde{a}_{q}\times\widetilde{a}_{w}$ hyperplane quantities by a superscript $(2)$.

[TABLE]

where the dots indicate integrals over every degree of freedom, except $\widetilde{a}_{q}$ and $\widetilde{a}_{w}$. Here, we used that $\mathbf{\Omega}$ is skew-symmetric in correlation-identity coordinates, i.e. $\widetilde{\mathbf{\Omega}}_{q,w}=-\widetilde{\mathbf{\Omega}}_{w,q}$ Weiss (2003).

Therefore, the current $\widetilde{j}^{(2)}=(\widetilde{j}^{(2)}_{q},\widetilde{j}^{(2)}_{w})^{T}$ on the $\widetilde{a}_{q}\times\widetilde{a}_{w}$ hyperlane is given by the lower dimensional analogue to Eq. (25), $\widetilde{j}^{(2)}=\widetilde{\mathbf{\Omega}}^{(2)}\vec{\widetilde{a}}\widetilde{\rho}^{(2)}(\vec{\widetilde{a}})$ with $\widetilde{\mathbf{\Omega}}^{(2)}$ being an antisymmetric $2\times 2$ matrix with off-diagonal elements $\widetilde{\mathbf{\Omega}}_{q,w}$ and $-\widetilde{\mathbf{\Omega}}_{q,w}$. The term $(2\pi)^{-1}e^{-\frac{\widetilde{a}_{q}^{2}+\widetilde{a}_{w}^{2}}{2}}$ in Eq. (26) is simply the marginal joint distribution $\widetilde{\rho}^{(2)}$ of $\widetilde{a}_{q}$ and $\widetilde{a}_{w}$. The eigenvalues of $\widetilde{\mathbf{\Omega}}$ hence read $\lambda_{\pm}=\pm i\widetilde{\mathbf{\Omega}}_{q,w}$, which implies that $\widetilde{j}^{(2)}(\widetilde{a}_{q},\widetilde{a}_{w})$ is purely rotational. From this analysis, we conclude that the apparent cycling frequency in the $\widetilde{a}_{q}\times\widetilde{a}_{w}$ hyperplane in correlation-identity coordinates is given by

[TABLE]

We note that the matrix elements of $\widetilde{\mathbf{\Omega}}_{q,w}$ must be calculated in the full-dimensional system.

In mode spaces with more than two dimensions, we obtained frequencies numerically (see Fig. 9). To derive analytical predictions for the net cycling frequency, we now consider a two-dimensional system constructed from $a_{q}\times a_{w}$. Here, we do not need to transform the system into correlation-identity coordinates because a 2D system only cycles at a single frequency. Intuitively, this simple 2D model should result in an overestimate of the corresponding frequency of the actual higher-dimensional system since this approach excludes transitions out of the considered plane, which would tend to reduce the in-plane cycling frequencies. Indeed, the two-dimensional frequency $\omega_{q,w}^{\text{2D}}$ appear to be an upper bound for the actual frequencies in the higher-dimensional system, as shown in Fig. 9.

We may, nevertheless, use the approximative calculation to explore under which circumstances we expect detailed balance to be broken. In our two-dimensional subsystem, $\omega_{q,w}^{\text{2D}}$ can be directly calculated as the positive eigenvalue $\lambda_{+}$ of $\mathbf{\Omega}$ in Eq. (25) computed for a system consisting only of modes $a_{q}$ and $a_{w}$ and reads

[TABLE]

this result shows that for (i) diagonal coupling matrices or (ii) equal relaxation times, detailed balance is always maintained. Furthermore, Eq. (28) illustrates how the thermal noise floor, which is constant along the filament, affects the non-equilibrium cycling and may even conceal it in the limit of weak motor action.

As predicted, detailed balance appears to hold, i.e. not be broken, in planes of even and odd modes, for example $a_{2}\times a_{3}$, for the regular geometry of interaction $\{s_{n}\}$. The magnitude and structure of the currents numerically calculated for the pairs $a_{2}$, $a_{6}$ and $a_{3}$ and $a_{5}$ match our analytical predictions in Fig. 5. For higher modes, which correspond to shorter length scales, the finite width of the motor-filament interaction in the simulation results in a coupling of even and odd mode amplitudes.

From an experimental point-of-view, where noise in the data is often limiting, it is desirable to determine the value of a single statistical variable instead of spreading sampled data over a number of variables. The net cycling frequency $\omega_{q,w}$ or the cross-correlations $\langle a_{q}(t)a_{w}(t^{\prime})\rangle$ are examples of such single-variable non-equilibrium indicators. The net cycling frequency in a plane $a_{q}\times a_{w}$ can be inferred from empirical mode traces using the definition $\omega_{q,w}=\langle\dot{\phi}_{q,w}\rangle$, where $\phi_{q,w}(t)=\text{atan}(a_{w}(t)/a_{q}(t))$ is the polar angle. This leads to the observed frequency in the plane

[TABLE]

In Fig. 9 (a), we compare the prediction for the two-dimensional net cycling frequencies $\omega^{\text{2D}}$ based on Eq. (28) with the prediction for the apparent frequencies in a high-dimensional space from Eq. (27) and frequencies obtained from simulations. The distributions of angular displacements $\Delta\theta$ are shown on a log-scale in Fig. 9 (b). As the distributions widen with increasing mode number, it becomes more difficult to precisely estimate the average net cycling frequency.

While it is tempting to interpret $\omega_{q,w}$ as a metric for the distance to equilibrium, the cycling frequencies are actually determined by two features of the system: the natural relaxation times of the respective modes and the strength of driving relative to the thermal background. Thus, the cycling frequencies provide a comprehensive measure of the non-equilibrium dynamics of different scales.

VI Summary

In this paper, we extended existing models describing the stochastic dynamics of semiflexible probe filaments embedded in viscoelastic media to incorporate effects of non-equilibrium force generation in active media. This model would, for example, describe probe filaments such as microtubules or single-walled carbon nanotubes in the cell’s actin cytoskeleton activated by myosin motors. In particular, we derived analytical descriptions of dynamic normal mode cross- and autocorrelations that can be used to characterize the stochastic forces the probe filaments are subject to in the active networks under steady state conditions.

We showed that motor-induced forces may lead to probability currents in phase spaces spanned by dynamic normal mode amplitudes, indicating a breaking of detailed balance. The structure of these currents is closely related to the active force profile along the backbone of the probe. We confirmed these analytical results with currents inferred from Brownian dynamics simulations of semi-flexible filaments driven by molecular motors. In steady-state, we showed that these divergence-free currents form a net cyclic probability flux in the phase space of mode amplitudes with a characteristic frequency $\omega_{q,w}$. We discussed implications of the dimensionality of mode space on the relation between measured frequencies and active force profiles and derived an approximate analytical expression for these cycling frequencies in the white-noise limit of motor action. Furthermore, we indicated special cases in which we expect detailed balance to hold in spite of non-equilibrium driving: In planes $a_{q}\times a_{w}$ spanned by modes with (i) equal relaxation times or (ii) vanishing correlations, detailed balance will hold.

Many predictions of our model, such as the magnitude of probability currents, depend on the geometry of the network around the probe filament, which might be difficult to chart. By contrast, variables like relaxation times or correlation functions of the normal mode amplitudes appear to be more robust against geometrical details of the active forcing. We therefore introduced the breaking of Onsager’s time reversal symmetry as a novel non-equilibrium marker, which only involves temporal aspects of the correlation function. In fact, we show that under non-equilibrium conditions, reciprocity of the correlation function breaks down. This can be used as a powerful tool to detect and quantify non-equilibrium behavior using the dynamics of probe filaments.

In summary, mode cross-correlations and cycling frequencies, may serve in the future to detect non-equilibrium dynamics and characterize the spatial distribution and the temporal behavior of motor-induced forces acting on individual probe filaments. Estimating probability currents from sample trajectories requires rather large data sets, however. The availability of photostable fluorescent filaments, such as single-walled carbon nanotubes, will bring the proposed experiments within practical reach. More generally, our results could be applied to any extended object with non-thermal fluctuations such as chromosomes, membranes, cellular organelles, whole cells or even tissues.

Acknowledgements.

We thank F. MacKintosh, N. Fakhri, F. Mura, F. Gnesotto, and G. Gradziuk for helpful discussion. This research was supported by the National Science Foundation under Grant No. NSF PHY11-25915,by the German Excellence Initiative via the program NanoSystems Initiative Munich (NIM) (C.P.B.), the GRK2062 (C.P.B) and the Deutsche Forschungsgemeinschaft (DFG) Collaborative Research Center SFB 937 (Project A2), the European Research Council Advanced Grant PF7 ERC-2013-AdG, Project 340528 (C.F.S), and the Cluster of Excellence and DFG Research Center Nanoscale Microscopy and Molecular Physiology of the Brain (CNMPB) (C.F.S.).