Simulations of dynamically cross-linked actin networks: Morphology, rheology, and hydrodynamic interactions

Ondrej Maxian; Raúl P. Peláez; Alex Mogilner; Aleksandar Donev

doi:10.1371/journal.pcbi.1009240

Abstract

Cross-linked actin networks are the primary component of the cell cytoskeleton and have been the subject of numerous experimental and modeling studies. While these studies have demonstrated that the networks are viscoelastic materials, evolving from elastic solids on short timescales to viscous fluids on long ones, questions remain about the duration of each asymptotic regime, the role of the surrounding fluid, and the behavior of the networks on intermediate timescales. Here we perform detailed simulations of passively cross-linked non-Brownian actin networks to quantify the principal timescales involved in the elastoviscous behavior, study the role of nonlocal hydrodynamic interactions, and parameterize continuum models from discrete stochastic simulations. To do this, we extend our recent computational framework for semiflexible filament suspensions, which is based on nonlocal slender body theory, to actin networks with dynamic cross linkers and finite filament lifetime. We introduce a model where the cross linkers are elastic springs with sticky ends stochastically binding to and unbinding from the elastic filaments, which randomly turn over at a characteristic rate. We show that, depending on the parameters, the network evolves to a steady state morphology that is either an isotropic actin mesh or a mesh with embedded actin bundles. For different degrees of bundling, we numerically apply small-amplitude oscillatory shear deformation to extract three timescales from networks of hundreds of filaments and cross linkers. We analyze the dependence of these timescales, which range from the order of hundredths of a second to the actin turnover time of several seconds, on the dynamic nature of the links, solvent viscosity, and filament bending stiffness. We show that the network is mostly elastic on the short time scale, with the elasticity coming mainly from the cross links, and viscous on the long time scale, with the effective viscosity originating primarily from stretching and breaking of the cross links. We show that the influence of nonlocal hydrodynamic interactions depends on the network morphology: for homogeneous meshworks, nonlocal hydrodynamics gives only a small correction to the viscous behavior, but for bundled networks it both hinders the formation of bundles and significantly lowers the resistance to shear once bundles are formed. We use our results to construct three-timescale generalized Maxwell models of the networks.

Author summary

The cytoskeleton is composed of semiflexible, inextensible actin filaments, dynamic cross linkers, and molecular motors, and makes the primary contribution to the structural properties of the cell. Despite its being so fundamental to cell biology, the biological complexity of the cytoskeleton hinders our ability to understand its mechanical properties through in vitro experiments. In this paper, we perform microscopic simulations of actin fibers and transient cross linkers to quantify the principle timescales involved in the network, study how these timescales influence the morphology and rheology of the system, and examine the role of hydrodynamic interactions in cytoskeletal networks. We find three principle timescales which we associate with fiber flexibility, cross linker detachment, and network remodeling, respectively. We show that the morphology of the network is more important on longer timescales, where the viscosity of links inside of fiber bundles is enhanced. We also show that hydrodynamic interactions reduce the stress inside of bundles because of entrainment flows. Finally, we propose a continuum model which can be used to coarse-grain our agent-based simulations and enable modeling of larger systems.

Citation: Maxian O, Peláez RP, Mogilner A, Donev A (2021) Simulations of dynamically cross-linked actin networks: Morphology, rheology, and hydrodynamic interactions. PLoS Comput Biol 17(12): e1009240. https://doi.org/10.1371/journal.pcbi.1009240

Editor: Roeland M.H. Merks, Leiden University Faculty of Science: Universiteit Leiden Faculteit der Wiskunde en Natuurwetenschappen, NETHERLANDS

Received: June 29, 2021; Accepted: November 12, 2021; Published: December 6, 2021

Copyright: © 2021 Maxian et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: Code and input files for all of the simulations are available at https://github.com/stochasticHydroTools/SlenderBody.

Funding: This work is supported by the National Science Foundation via GRFP/DGE-1342536 (to O.M.), DMS-2052515 (to A.M. and A.D.), DMS-1953430 (to A.M.), and RTG/DMS-1646339 (to A.D.). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

In most cells, as much as 10% of all protein is actin [1]. The majority of actin is the F-actin cytoskeleton − a gel made of rapidly-turning-over (assembling and disassembling) actin filaments (which we will call fibers), which are dynamically cross linked by a vast host of actin binding proteins (which we will call cross linkers or CLs). The actin cytoskeleton is largely responsible for the cell’s shape, movements, division, and mechanical response to its external environment [2]. Thus, understanding the mechanical and rheological properties of dynamically cross-linked cytoskeletal networks is the first step to understanding the mechanical properties of the cell at large.

Several in vitro experimental techniques, including active poking, parallel plate shearing, and embedded-microbead tracking, have characterized the mechano-rheological behavior of cross-linked actin gels [1]. Depending on the experimental conditions, these studies report viscoelastic moduli in a wide range, from 0.1 to hundreds of Pascals [3–6]. The mixed viscoelastic mechanical response of the densely cross-linked gel is expected: simply speaking, at short time scales the network of elastic fibers and CLs can be considered permanently interconnected and deform elastically, while at long time scales dynamic CLs are expected to connect fibers only transiently, enabling the brief storage of elastic energy of deformations. This elastic energy is dissipated after the CLs detach, which causes effective viscous behavior on long timescales. Often, the measured elastic modulus for actin gels is about an order of magnitude greater than the viscous one [3, 4, 6]; however, some studies measure elastic and viscous moduli of similar magnitude [7], and in some cases the elastic modulus is smaller than the viscous one [8]. Both moduli are increasing functions of actin and CL concentrations [9, 10]; see especially [10] for master curves over a wide range of CL and actin concentrations.

To understand the intrinsic timescales in the gel, many experimental studies examine the dependence of the loss and storage moduli on the frequency of oscillatory shear deformation ω. All experimental studies agree that the behavior of transiently cross-linked networks on long timescales (low frequencies) is qualitatively viscous both in vitro [11, 12] and in vivo [7], with many studies reporting a decay rate of ω^1/2 for both moduli at low frequencies in the absence of actin turnover [11, 13, 14]. As the frequency increases, some studies report a monotonically-increasing viscous modulus [10], while others display a curious local maximum and minimum in the viscous modulus data [13, 15, 16]. Proceeding to the high-frequency limit, there are again contradictions, as some studies report a viscous modulus that scales as ω^0.75 [10] (consistent with the thermal fluctuations of a single semiflexible filament [6, 17]), while others observe a viscous-fluid-like scaling of ω¹ [18, 19], and still others yield a scaling of ω^1/2 [11, 12, 20], which is characteristic of the Rouse model of polymer physics (fluctuating beads joined by harmonic springs) [21]. Thus the short and intermediate timescale behavior, and its dependence on microscopic parameters, is still an open question, as is the exact meaning of “long” and “short” timescales (i.e., when these observed scalings begin to dominate).

Part of the reason for the variations in experimental data is that each actin binding protein produces a unique change in network morphology [22], which in turn uniquely affects the viscoelastic moduli. For compliant CLs such as filamin, Kasza et al. observe a strong frequency-dependence of ω^0.7 for the viscous modulus and a weak power law scaling of the elastic modulus, with both taking values in the range 0.1–1 Pa [23]. As the CLs get shorter and stiffer, the frequency-dependence in the elastic [8, 10, 23, 24] and sometimes even viscous modulus disappears [23], which suggests a change in the network morphology with the type of cross-linker. It is shown in [11] that changes in the moduli with chemical cross linking are only relevant when the cross-linker length is shorter than the mesh size, because these kind of CLs group the fibers into bundles, thereby changing the macroscopic structure of the network.

Wachsstock et al. [5] study the relationship between mechanical response and morphology using actin filaments with α-actinin cross linkers. Their experiments and modeling show a transition from elastic, solid-like behavior to viscous, fluid-like behavior as the α-actinin and actin concentrations increase, which corresponds to the formation of bundles within the network. These bundles are no longer entangled in a complex network and are free to slide past each other, which leads to viscous behavior. This observation is compatible with the work of Tseng et al. [25], which shows that more homogeneous networks are more elastic in nature, and suggests that the cell makes the strongest elastic structures by combining bundling and cross-linking proteins to form a cross-linked network of bundles [26].

In view of this experimental complexity, many modeling studies address scaling of the cross-linked actin gel’s mechanical moduli with kinetic and microscopic mechanical parameters. Some of the rheological behavior of pure actin networks can be explained with semiflexible polymer theory [27, 28]. For instance, strain-stiffening behavior in filament networks can be accounted for by an entropic model, in which each thermally fluctuating polymer resists affine (stretching) deformations, or an enthalpic model, in which the filaments bend before they stretch, allowing them to reorganize along the direction of shear [29]. Semiflexible filament theories of this nature have recently been extended to transiently cross-linked networks [16, 30, 31], and it has been shown that a broad spectrum of relaxation times appears for a cross linker with a single unbinding rate, because larger and larger wavelength bending modes are able to relax as cross linkers unbind [14, 30]. These theories may explain the soft glassy rheology observed in living cells [14] and the low-frequency regime where the moduli scale with the square root of frequency [16, 30].

Speaking more broadly, the existing models of cross-linked networks fall into two categories: continuum and agent-based (discrete) models. There are a number of continuum models available for the cell cytoskeleton, and in fact there has been a lot of work in recent years on extending the immersed boundary method [32] to cytoskeleton-like fluids. For example, Karcher et al. employ a continuum finite element model to measure the stress induced by magnetic-bead-forcing [33] (a similar model was also used in [34]). They model the cytoskeleton as either a Maxwell (spring and dashpot in series) or Kelvin-Voigt (spring and dashpot in parallel) material and report a large sensitivity of the mechanical behavior on the choice of model. Other studies have approximated the cell cytoskeleton as a poroelastic [35], porous viscoelastic [36], or Brinkmann [37] fluid. In either case, as discussed in [1, 33], continuum models are only as good as their fits to the data, and they in general have difficulty relating macroscopic parameters (such as porosities, discrete timescales, and stiffnesses) to microscopic parameters. In fact, it has been suggested that the cell possesses a continuum of relaxation timescales [14, 38], which would invalidate any continuum model with a discrete number of elements.

Discrete models tend to suffer the opposite problem in that they are detailed and realistic, but are harder to analyze and do not readily extrapolate to macroscopic systems. Most of these models have focused on simulating the steady state structure of actin networks, which presumably could help explain some of the observed rheological behavior. Head et al. [39], for instance, solve an energy minimization problem for straight filaments in two dimensions and use their results to characterize the transition from affine (stretching-dominated) to non-affine (bending-dominated) deformation in cross-linked networks. In a more detailed approach, Kim et al. [40] use Brownian dynamics simulations on a bead-spring model of actin to show that bundling is reduced when CLs are forced to bind in the perpendicular direction (as they do when fascin is the cross-linker). They propose that the ability of a network to flow on long timescales is due to the breaking of CLs with shear, followed by their reformation on a timescale determined by network reorganization and not individual CL binding. There have also been a number of computational studies on the structure and contractile behavior of actomyosin networks, which show that a critical concentration of α-actinin cross-linkers can combine with myosin motors to form ordered bundles of actin filaments with varying polarity [41, 42].

Two recent studies [31, 43] combine computational rheology with microscopic Brownian dynamics simulations to measure the rheological properties of transiently cross-linked actin networks without filament turnover and without hydrodynamic interactions between the filaments. In these studies, a local maximum in the viscous modulus is observed [15], which was shown to correspond to a maximum in CL binding and unbinding events [43]. However, the networks considered in these studies are either homogeneous [43](Fig. 2) or highly bundled [31](Fig. 1), and we know of no detailed simulations that have reported the rheological properties over a range of structures and microscopic parameters. Furthermore, due to the absence of turnover, in both of these studies the structure being simulated is not in equilibrium and evolves significantly with the imposed flow, and so the measured low-frequency viscous and elastic moduli depend on history and how the measurement is performed, and are consequently not well-defined.

Despite, and maybe even because of, the large volume of experimental and theoretical studies on cross-linked actin gel mechanics, a number of questions still remain open. We focus on three of them: What role does the morphology of the network play in the mechanical response? Could there be a significant contribution from nonlocal hydrodynamic interactions? And can we inform a simple and practical spring-and-dashpot model by detailed microscopic simulations?

This paper is built around answering these questions. To do so, we first review our formulation for inextensible, semiflexible filaments [44] and extend it to the case of transient cross-linking. We model filament polymerization dynamics by turning over filaments with a characteristic time τ_f and introduce an operator splitting scheme to update the system by turning over the filaments, updating the cross-linked network, and moving the filaments in a sequential order. Since previous experimental and modeling studies show that the mechanical behavior of actin networks is most affected by the presence of short, stiff CLs which induce morphological changes in the network, we will focus our analysis on CLs with spring stiffness 10 pN/μm and a 50 nm rest length, thereby mimicking α-actinin. We show that changing the turnover time τ_f with all other parameters fixed gives two model systems: an isotropically cross-linked network, or homogeneous meshwork, and a composite network of bundles embedded in the meshwork (here we follow the classification scheme of [45]). Unlike previous simulation studies [31, 43], the networks that we consider here are in a dynamic steady state, which allows us to precisely quantify their steady-state mechanical properties. In particular, we use small amplitude oscillatory shear (SAOS) rheology over a range of biologically-relevant frequencies (0.01 Hz to 100 Hz) to show that the networks have three characteristic timescales on which the links go from rigid to flexible, permanent to dynamic, and viscoelastic to purely viscous. It is on this last timescale that remodeling of the network occurs. Here by “biologically relevant,” we mean timescales from ∼0.1 sec (the fastest relevant microscopic changes in the network structure, like unbinding of a CL [7]) to tens of seconds (turnover time, [46, 47]). Because of this, in vivo rheological experiments often focus on the range from 0.1 to 100 s [7]—we follow their example.

This paper is, to our knowledge, the first to incorporate many-body hydrodynamic interactions (HIs) in the study of transiently cross-linked actin networks. Here we once again find results that depend on the timescale under consideration. At high frequency, we show that the viscous modulus scales linearly with ω, and that accounting for nonlocal hydrodynamics simply changes this coefficient by a fixed percentage. In contrast, at low frequency, we show that nonlocal hydrodynamic interactions lead to a significant decrease in the moduli for structures with a significant degree of bundling, and that the decrease can be explained by the disturbance flows inside a bundle. We conclude this paper by using our knowledge of the timescales and role of nonlocal hydrodynamic interactions to suggest coarse-grained continuum models for the transiently cross-linked actin gel.

Materials and methods

We first introduce our simulation method for passively cross-linked actin gels. This begins with a brief review of our recently-developed spectral method for semiflexible, inextensible fibers in the presence of cross linkers [44]. We next introduce our models of dynamic cross-linking and actin turnover, both of which are based on the next reaction stochastic simulation algorithm. We then discuss the three approximations to the mobility (force-velocity) relationship: (non-isotropic) local drag, intra-fiber hydrodynamics, and full nonlocal hydrodynamics, as well as how we evaluate them. Most of the details of this can be found in [44]; here we only elaborate on the modifications we have made for efficient simulation of a bundled many-fiber gel. We conclude this section with our temporal integration (network evolution) strategy, which uses operator splitting to turn over the fiber positions, update the cross-linked network, and update the fiber positions in a sequential order.

Semiflexible, inextensible fibers

Since actin fibers are slender and semiflexible, we can represent them by smooth Chebyshev interpolants X(s), where s ∈ [0, L] is an arclength parameterization of the fiber centerline and L is the fiber length. We denote the tangent vector by τ(s) = X_s(s). For inextensible filaments such as actin, we have the constraint (1) for all s and t. To enforce the constraint, we introduce a Lagrange multiplier force density λ(s, t) on each fiber.

At every instant in time, each fiber resists bending with bending force density (per unit length) (2) where the constant linear operator gives f^κ with the “free fiber” boundary conditions [48] (3)

In [44](Sec. 4.1.3), we discuss how to discretize the operator to spectral accuracy in a manner consistent with the boundary conditions Eq (3) using a method called rectangular spectral collocation [49, 50].

Beginning with a single fiber, let us introduce the mobility operator and a spatially and temporally varying background flow u₀(x, t). Then, accounting for the constraint force, bending force, and any other external forces f^(CL) (e.g., those coming from attached CLs), the evolution equation on the fiber is given by (4) subject to the boundary conditions Eq (3) and inextensibility constraint Eq (1). In [44], we show how to obtain λ and enforce inextensibility via the principle of virtual work. We discuss different approximations to in a later section.

Dynamic cross-linking model

We now specify a method to compute the external or cross-linking force density f^(CL) in Eq (4). Here we extend [44] to account for transient cross linking, where the links can come on and off through the course of a simulation. To model this, we assume that the density of cross-links is sufficiently high that CLs do not get locally depleted by binding, and that the links diffuse rapidly. By “rapid” diffusion, we mean that the Brownian motion of the links is fast relative to the residence time of a single link, which is on the order of one second (see Table 1). To justify this, we consider the translational diffusion of a rod of length ℓ = 50 nm and radius a = 4 nm. The translational diffusion coefficient of the rigid rod can be approximated by Broersma’s relations as [51] (5)

Download:

Table 1. Simulation parameters.

https://doi.org/10.1371/journal.pcbi.1009240.t001

For our parameters in Table 1, we obtain ζ = 2.53, γ_∥ = 1.25, γ_⊥ = 0.19, and D_t = 0.15 μm²/s. This implies a center of mass displacement of 〈r²〉 = 6D_tt ≈ 0.9t. Thus the time for a CL to travel a mesh size of r = 0.2 μm is about 0.04 seconds, which is much faster than the one second lifetime of each link. These assumptions considerably simplify our simulations, as we do not have to keep track of explicit CLs, but only the number of links bound to each fiber, which can be tuned by setting the ratio of the binding and unbinding rates.

There are four reactions in the model (parameters are summarized in Table 1):

The attachment of one end of a free CL to a fiber, which occurs with rate k_on (units 1/(length × time)). Here k_on accounts for the rate of diffusion of the CL until one end gets sufficiently close to a fiber, and then the reaction for the CL to bind to the site.
The detachment of a singly-bound CL end to become a free-floating CL, which occurs with rate k_off (1/time). This is the reverse of reaction 1.
The attachment of the second end of a singly-bound CL to another (distinct) fiber, which occurs with rate k_on,s (units 1/(length × time)). This necessarily follows after reaction 1. The rate k_on,s accounts for the diffusion of the second end to find a fiber within a sufficient distance, and then the binding reaction for the end to latch on to the fiber. We model the thermal fluctuations in the CL length by allowing a singly-bound CL to bind its free end to another attachment site on a different fiber that is within a distance interval (6) and kT ≈ 4 × 10⁻³ pN ⋅ μm. In this paper, we will fix K_c = 10 pN/μm [58], so that the distance a link can stretch is δℓ = 0.02 μm. This is 40% of the rest length of ℓ = 0.05 μm.
The detachment of a single end of the doubly-bound CL, so that it becomes singly-bound, which occurs with rate k_off,s (1/time). This is the reverse of reaction 2.

A schematic of the four reactions is shown in Fig 1a. In our model, we do not account for any dependence of the binding/unbinding rates on the CL stretch.

Download:

Fig 1. Dynamic of model and resulting bundling behavior.

(a) Our model of dynamic cross-linking with fiber turnover. We coarse-grain the dynamics of individual CLs into a rate for each end to bind to a fiber. The first end (purple) can bind to one fiber at any binding site. Once bound, we account for thermal fluctuations of the CL length by allowing the second end (red) to bind to any other fiber within a distance of the first end, where ℓ is the CL rest length and K_c is the CL stiffness. To model actin turnover, we allow each fiber to disassemble with rate 1/τ_f and for a new (straight) fiber to assemble in a random place (nascent fiber is shown in green) at the same time. (b) Consecutive simulation snapshots illustrate how the model reproduces the bundling behavior characteristic of an actin mesh cross-linked with α-actinin. A pair of fibers that are close enough to be crosslinked by a thermally stretched CL are pulled together when this CL relaxes to its rest length. This brings the fibers closer together, promoting binding of additional CLs. When multiple CLs along the inter-fiber overlap relax to their rest length, cross-linked pairs of fibers align and stay close together, making a bundle.

https://doi.org/10.1371/journal.pcbi.1009240.g001

To simulate these reactions stochastically, we employ a variant of the Stochastic simulation/ Gillespie algorithm [62, 63]. To advance the system by a time step Δt, we calculate the rate of each reaction and use the rate to sample an exponentially-distributed time for each reaction to occur. We choose the minimum of these times, advance the system to that time, and then recalculate the rates for the reactions that were affected by the previous step. We organize the event queue efficiently using a heap data structure [64], so that we can compute the next reaction in time.

We break each filament into N_s possible uniformly-spaced binding sites, where each of the sites represents a binding region of size L/(N_s − 1) = Δs_u. To properly resolve the possible connections between fibers, we require . We then compute the rates of each reaction in the following way:

Single end binding: the rate of a single end binding to one of the N_s attachment points on one of the F fibers is k_onΔs_u. We implement this as a single event with rate N_sFk_onΔs_u, where the specific site is chosen uniformly at random once the event is selected. Multiple CLs are allowed to bind to a site, since each site represents a fiber segment of length Δs_u that is typically longer than the minimal possible distance between two CLs.
Single end unbinding: we keep track of the number of bound ends on each of the N_sF sites. If a site has b ends bound to it, we schedule a single end unbinding event from that site with rate k_offb. We do this for each site separately, so that this reaction contributes N_sF events to the event queue.
Second end binding: we first construct a neighbor list of all possible pairs of sites (i, j) that are on distinct fibers and separated by a distance ℓ + δℓ or less, using a linked list cell data structure [65]. If site i has b free ends bound to it, we schedule a binding event for that pair of sites with rate bk_on,sΔs_u. Notice that this naturally implies a zero binding rate if there are no free ends bound to site i. We treat the (j, i) connection separate from the (i, j) one for simplicity at the cost of increasing the number of events in the queue.
Second end unbinding: if a link exists between sites i and j, it can unbind from one of the sites with rate k_off,s. Letting β be the number of links in the system, the total rate for which bound links unbind from one end is 2βk_off,s. We therefore schedule a single event with rate 2βk_off,s, and, once it is chosen, select the link and end to unbind randomly and uniformly.

In the rest of this paper, we will denote the connectivity of the network (number of single ends bound to each site, list of sites connected by CLs) as Y(t). The CL force in Eq (4) is then a function of the fiber positions and network configuration, which we denote as f^(CL)(X;Y).

We use our previous work [44](Sec. 6.1) to define the cross-linking force density between two fibers X⁽ⁱ⁾ and X^(j) in a manner that can preserve the problem smoothness and therefore the spectral accuracy of our numerical method. Suppose that a CL connects two fibers by attaching to arclength coordinate on fiber i and on fiber j. Then the force density on fiber i due to the CL is (7) where K_c is the spring constant for the CL (units force/length), ℓ is the rest length, and δ_h is a Gaussian density with standard deviation σ. As discussed in [44](Sec. 6.1), the Gaussian width σ depends on how many points are used to discretize the fibers. In this paper, we use N = 16 points in the fiber discretization, and so we also use the required σ/L = 0.1 to preserve smoothness. In the limit σ → 0, the Gaussian density becomes a Dirac delta function, and the force density Eq (7) becomes the expected expression for a linear spring connecting the points and . This model is of course an approximation; as shown in [58], the stress-strain relationship for α-actinin, which is the model CL we consider here, is nonlinear and exhibits hysteresis for loading and unloading. We follow a number of other modeling studies [43, 66] in approximating α-actinin as a linear spring, although we do not include an angular stiffness, since recent experimental results [67] indicate that the dynamics of the α-actinin-actin bond are insensitive to rotation.

Fiber turnover

We account for the turnover of actin filaments by removing a filament (and all of its attached links) from the system and rebirthing a new, straight filament with a random location and orientation elsewhere in the system (see schematic in Fig 1a). Our reasoning for this is two-fold: first, it is simple to implement, and second, our simplifying assumptions are supported by experimental data. Experiments have shown that pointed-end depolymerization from one end of an actin filament is too slow to explain the observed rates of filament turnover in vivo [68], and that depolymerization actually occurs in bursts as filaments abruptly sever into smaller pieces [69]. In this sense, our model of depolymerization as an instantaneous event is supported by the data. For polymerization, the process is normally linear, with monomers being added at the plus end until it is capped. It has been shown experimentally [70], however, that the combined time of growth, capping, and disassembly is still significantly smaller than the time in which the filament is at a finite length, so we assume the former to be instantaneous relative to the latter. We refer to [71](Sec. 2.5) for a mathematical treatment of linear polymerization.

Having chosen this model, fiber turnover can be added as another reaction in our stochastic simulation algorithm. Defining the mean turnover time of each fiber as τ_f, we add a fiber turnover reaction to our reaction list with rate F/τ_f. If this reaction is chosen, we select the fiber to turnover randomly and uniformly. We emphasize that our time steps are at most 10⁻³ seconds while the turnover times are of the order seconds, so we almost always see zero or one fiber turnover events per time step.

Mobility evaluation

We complete our description of the evolution equation Eq (4) by specifying the mobility operator . Our approach is fully described in [44](Sec. 2, Appendix A) and is based on an improved version of the traditional slender body theory for Stokes flow [48, 72, 73]. For a single fiber, the self-mobility gives the fiber velocity in Eq (4) as (8) where r(s, s′) = X(s) − X(s′), r = ‖r‖, and . Thus the velocity on the fiber can be written as a leading order local drag term (the first line), which accounts for forcing localized around X(s), plus a “finite part” integral remainder term (the second line), which accounts for the intra-fiber nonlocal hydrodynamic interactions. In Eq (8), c(s) is a local drag coefficient which has a logarithmic dependence on the fiber radius a as (9) away from the fiber endpoints. Near the endpoints, Eq (9) as written is singular, and we regularize it over a distance δ = 0.1L as discussed in [44](Sec. 2.1). Throughout this paper, we will fix a = 4 nm [52].

When there are multiple fibers, the velocity Eq (8) has to be modified to account for the flows induced by fibers j ≠ i on fiber i. Denoting this flow by U_NL and indexing the ith fiber with a superscript (i), we use a line integral of the Rotne-Prager-Yamakawa (RPY) hydrodynamic tensor over all other fibers to obtain the total disturbance flow (10)

The RPY tensor S_RPY, which involves the Stokeslet and doublet singularities of Stokes flow, is a specific form of a symmetrically regularized Stokeslet, as explained in more detail in [44](Sec. 2) and [74]. The constant a* is proportional to the fiber radius a, and is chosen to match the slender body theory Eq (8) with the RPY integral for a single fiber [44](Eq (15), note that in that equation a = ϵL).

Our strategy for evaluating the total velocity on each fiber is documented in [44](Sec. 4). The most computationally intensive part of this calculation is the evaluation of the nonlocal integrals Eq (10). These integrals present challenges because they naively cost operations to evaluate, where F is the number of fibers, and because they can be nearly singular when fibers are close together. For triply periodic systems, which we consider in this paper, we reduce the cost to using the positively split Ewald method of [75, 76]. For a summary of how this method is adapted to sheared domains, see [44](Sec. 4.3).

In our previous work, we dealt with the nearly singular nature of the integrals for nearly contacting fibers by correcting the result from Ewald splitting with the special quadrature method of [77]. While this is the most accurate way to evaluate Eq (10), our recent experiments show that the fastest way to compute the integrals Eq (10) with sufficient accuracy is via oversampled direct quadrature on a GPU. Specifically, we discretize all integrals for j ≠ i using Clenshaw-Curtis quadrature with weights w, (11)

We then perform the summation with periodic boundary conditions using the positively split Ewald method [44, 75, 76], which we implement on a GPU in the UAMMD library [78], with N_u sufficiently large to resolve all of the integrals in Eq (10).

Following [44], to investigate the role of nonlocal hydrodynamic interactions, we vary the model used to compute the fluid velocity on the fibers. Our first option, which we refer to as the local drag model, is to include only the first line of Eq (8). Note that this mobility is non-isotropic, since the c(s) term in Eq (8) gives approximately twice the velocity for a force in the tangential direction, so even the local drag model we use is a significant improvement over formulations based on an isotropic local drag coefficient [43, 79], which are also missing the logarithmic dependence in Eq (9). Another option is to include only the U_L term in each fiber’s velocity so that (nonlocal) hydrodynamic interactions between a fiber and itself (which for dilute suspensions are the most important) are accounted for, but not interactions between a fiber and the others. We term this the intra-fiber mobility. Finally, when all terms U_L + U_NL are included on each fiber, we refer to the mobility as the full (nonlocal) hydrodynamic mobility.

SAOS rheology

To quantify the viscoelastic behavior of the network, we measure the linear elastic (also called storage) modulus G′ and viscous (loss) modulus G′′. We employ small-amplitude oscillatory shear (SAOS) rheology to do this [44, 80], so that we impose a shear flow of the form (12)

The nonzero component of the rate of strain tensor for this flow is given by (13) while the nonzero component of the strain tensor is (14) where is the maximum strain in the system. The bulk elastic (G′) and viscous modulus (G′′) relate the stress to the strain (for the elastic modulus G′) and rate of strain (viscous modulus G′′) [80] via (15)

As discussed in [44](Sec. 6.2), the stress tensor for our system can be decomposed into a part coming from the background fluid and a part coming from the internal fiber stresses, (16)

Substituting the background flow Eq (12), we obtain the stress for the background (pure viscous) fluid (17)

A pure viscous fluid therefore has viscous modulus scaling linearly with ω as G′′(ω) = 2πωμ. This will be useful to us when we examine the behavior of the viscous modulus of the actin gel.

The contribution of the fibers to the stress tensor is the integral ∫X(s)f(s) ds, where f is the total force density on the fibers [48, 81], including the elastic forces exerted by the cross-linkers. This must be handled carefully in periodic boundary conditions, as the correct periodic copy of the fiber must be used to ensure the stress is symmetric. The formula we use is given in [44](Eqs. (124,125)).

Temporal integration and splitting

We have three different events to process at each time step: fiber turnover, link turnover, and fiber movement. We will treat them sequentially using a first order Lie splitting scheme [82]. For each time step n of duration Δt, from time nΔt to (n + 1)Δt, we

Turnover filaments to update configuration .
Process binding/unbinding events to evolve Y_n → Y_n+1.
Use the method of [44] to evolve .

The evolution of the fiber network in step 3 is performed using the method of [44](Sec. 4.5). This method obtains second-order accuracy and unconditional stability for fibers with permanent CLs using a combination of a linearly-implicit trapezoidal discretization for the fiber bending force and extrapolations for arguments of nonlinear terms (e.g., we use X_n+1/2,* = 3/2X_n − 1/2X_n−1 as the argument for the mobility and apply the nonlocal part of the mobility to an extrapolated constraint force λ_n+1/2,* = 2λ_n−1/2 − λ_n−3/2). We use this temporal discretization wherever it applies and is sufficiently stable.

However, extrapolations such as X_n+1/2,* are nonsensical for fibers that are turned over in the previous time step. For those fibers, we set and compute a guess constraint force λ_n+1/2,* by solving Eq (4) with the local drag mobility on each fiber that is respawned. Additionally, in our temporal discretization [44](Eq. (102)), the forcing inside the nonlocal mobility is treated entirely explicitly, since there we assume that the local drag part dominates the fiber motion. There are two ways this could go wrong: first, if a uniform, isotopic suspension is sufficiently concentrated, the temporal integrator becomes unstable. This was the case we already dealt with in [44], where we switched to implicit treatment of the bending force in all parts of the mobility and used an iterative solver. When the fibers are in bundles, however, we find that the approach of [44] is prohibitively expensive due to the ill-conditioning of the elastic force operator (which involves fourth derivatives) combined with the extrapolations for X_n+1/2,* and λ_n+1/2,*. For this reason, we switch to a first-order backward Euler discretization of fiber elasticity, where X_n+1/2,* is replaced by and the elastic force is computed on X_n+1. The backward Euler discretization is our preferred one for simulations with significant bundling and full hydrodynamics.

Finally, we discretize the stress tensor in a manner consistent with the numerical method, either as or as . The moduli G′ and G′′ are then evaluated by integrating the stress against the orthogonal sine and cosine functions [44](Eq. (126)).

Results and discussion

This section presents our results for the viscoelasticity of dynamically cross-linked actin networks. We use the parameters in Table 1 to perform the simulations, and show that changing the fiber turnover time gives a pair of dynamic steady states with varied degrees of bundling. We then use a stress relaxation test to show that the networks relax to a state of zero stress on a timescale of a few seconds. This is quantified more precisely with our rheological data, where we establish three discrete relaxation timescales and discuss the behavior on short, long, and intermediate timescales in more detail. Our discussion of the role of nonlocal hydrodynamic interactions (HIs) necessarily follows this, as the role of HIs is more significant on long timescales than short ones. We conclude this section by developing a continuum Maxwell-type model that is informed by our knowledge of the three timescales.

Our simulations use a background fluid with viscosity μ = 0.1 Pa ⋅ s, mimicking the larger viscosity of the cytoplasm [56]. We will, however, assume the background fluid is Newtonian, in contrast to actual cytoplasm, so that we can isolate the contribution of the fibers to the viscoelastic moduli.

Fiber turnover creates dynamic steady states with varied degrees of bundling

We first show that varying the turnover time τ_f of the F fibers creates a set of dynamic steady states with varying degrees of bundling. We also verify that the statistics are independent of the system size by considering three possible cubic periodic unit cells of different edge length L_d, all with the same mesh size (computed assuming that the fibers are distributed isotropically) μm. The fibers have length L = 1 μm throughout, so we will consider domains of size L_d = 2 μm (F = 200 fibers), L_d = 3 μm (F = 675 fibers), and L_d = 4 μm (F = 1600 fibers). We will show that simulations of a homogeneous meshwork (τ_f = 5 seconds) are suitably carried out in the smallest of the three systems, while runs where bundling is more considerable require the next-largest system (L_d = 3 μm).

As shown in Fig 1b, dynamic cross-linking leads to bundling of fibers, a phenomenon that is well documented experimentally for a variety of CL types [9, 83, 84]. To define a bundle, we use the network Y to map the fibers to a connected graph [85]. Unlike [85], we say that two fibers are connected by an edge if they are connected by two links a distance L/4 apart, the idea being that the fibers are sufficiently aligned in that case for the links to constrain their orientations. “Bundles” are then the connected components of this graph (two or more fibers per bundle). To better understand the network morphology, we define a bundle density (units bundles/μm³), where B is the total number of bundles, and also quantify the percentage of fibers in bundles. We define as the average number of links bound to each fiber (the total number of links is ).

Homogeneous meshwork morphology.

Beginning with turnover time τ_f = 5 s, we perform a steady state run to 10τ_f = 50 s using the three domain sizes. As shown in Fig 2a and S1 Video, the network is made primarily of fibers oriented isotropically (as can be verified by computing an order parameter of the type [85](Eq. (21))), with a few bundles of at most two to three aligned fibers. We quantify this more precisely in Table 2, where we see that there is on average one bundle per μm³, and that less than 10% of the fibers are in bundles. Because there are only a small number of bundles and no permanent structures over long timescales in this system, we classify this system as a homogeneous meshwork. We will report results for it using a domain size of L_d = 2 μm (we have verified that using L_d = 3 μm does not change the results significantly; see S1 Fig).

Download:

Table 2. Systems considered and their mean link and bundle densities.

https://doi.org/10.1371/journal.pcbi.1009240.t002

Download:

Fig 2. Morphology of the dynamic steady states.

The actin gel for (a) τ_f = 5 seconds and (b) τ_f = 10 seconds, both shown on a domain with edge length L_d = 3 μm. Fibers in the same bundle are colored with the same color. In (a), we observe a more homogeneous mesh. In this case, we use L_d = 2 μm in most simulations. In (b), we observe multiple bundles embedded into a mesh, and we use L_d = 3 μm in most simulations.

https://doi.org/10.1371/journal.pcbi.1009240.g002

Bundles embedded in meshwork morphology (B-In-M).

In Fig 2b, we show snapshots from the dynamic steady state with a (doubled) turnover time of τ_f = 10 s. We observe significantly more bundling and fiber alignment, as well as bundles with several (four to six) fibers in them. In Table 2, we see that the steady-state link density has more than doubled from the 5 s turnover case, and that the percentage of fibers in bundles has gone from about 10% for τ_f = 5 s to 25% for τ_f = 10 s. The fluctuations of link density in the smaller system (L_d = 2) are quite large (about 20%, see S2 Fig), which makes averaging too inaccurate. For this reason, for τ_f = 10 seconds we run a larger system with L_d = 3 μm, which we show in Fig 2b and S2 Video. We use the L_d = 4 system with 1600 fibers to verify that the finite system size has little effect; see, for example, Table 3.

Download:

Table 3. Network structure statistics for the runs over 120 seconds for the B-In-M network with ω = 1 Hz.

https://doi.org/10.1371/journal.pcbi.1009240.t003

Because this system has a significant number, but not more than half, of the fibers in bundles, and because the maximum bundle size is still at most ten fibers, we refer to this system as a bundles embedded in meshwork (B-In-M) morphology, where there are small bundles of a few fibers embedded in an otherwise homogeneous mesh.

Other systems (parameter variations).

In addition to varying the morphology through the turnover time, we also look at systems with the same morphologies, but different microscopic parameters. Table 2 gives our description of these systems. For a homogeneous meshwork, we consider τ_f = 5 s again, but this time with ten times smaller bending stiffness, so that κ = 0.007 pN ⋅ μm². Table 2 shows that the morphology in this case is the same as the homogeneous meshwork with κ = 0.07 pN ⋅ μm², so changing the bending stiffness has a minimal impact on the morphology. Indeed, the fibers remain relatively straight even with the smaller bending stiffness (see S3 Fig), just as they are straight in both the B-In-M and homogeneous meshwork systems in Fig 2. This implies that the degree of cross-linking/network connectedness is insufficient to induce filament bending, even in B-In-M geometries.

We also consider B-In-M systems with larger viscosity, μ = 1 Pa ⋅ s instead of 0.1, and another system where we double the rates of link binding and unbinding. Table 2 shows how we vary the turnover time in these systems to obtain a B-In-M morphology. When the viscosity is increased by a factor of 10, we only need to increase the turnover time by a factor of 1.5 to obtain a similar morphology (match as closely as possible). Meanwhile, doubling the rates forces us to cut the turnover time almost in half to obtain the same morphology. This provides our first indication that the morphology, or the timescale on which the network organizes itself into bundles, is strongly dependent on the link binding and unbinding rates and weakly dependent on the underlying fluid viscosity. This combination of behavior suggests we are operating in a regime where the CLs are attached for long enough to move the fibers into a quasi-steady state (since using twice the link turnover rate speeds up the process by a factor of two and making the dynamics slower by changing viscosity has little effect).

The linear response regime (LRR): Viscoelastic moduli and stress spectra.

Our goal for the rheology experiments is to measure the shear moduli of the network in its dynamic steady state, and in the linear response regime (LRR). In this section we briefly demonstrate how to find the LRR and what happens when we go beyond it, using the B-In-M geometry as an example.

Table 3 shows the steady state statistics and mean viscoelastic moduli for the B-In-M system with varying strains γ. In Table 3 and throughout this paper, we report the viscous modulus G′′ without the contribution from the viscosity of the Newtonian solvent . We observe from the table that the LRR is γ ≲ 0.05, since using γ > 0.05 leads to significantly larger moduli. In particular, Table 3 shows that larger strains of γ ≥ 0.10 induce higher bundle densities and the formation of more links. Furthermore, the links that do form in the larger-strain systems are themselves more strained, as the mean square link strain at γ = 0.1 and γ = 0.15 is 10% and 25% higher, respectively, than in the linear regime. The combination of these two factors gives large increases in the moduli when we reach the nonlinear regime. The viscoelastic moduli for larger strains are also subject to significantly larger uncertainties (e.g., the uncertainity for γ = 0.15 is at least double the uncertainity for γ = 0.1). The larger uncertainties come from fluctuations in the underlying microsctructure, as the fluctuations in the link density increase by 50% from γ = 0.01 to γ = 0.15.

Another way we can confirm that we are in the LRR is by looking at the Fourier spectrum of the stress. Specifically, we will write the stress as (18) and use the discrete Fourier transform to look at the amplitude of the coefficients for various γ. In S4 Fig, we show the spectrum for integer values of k with τ_f = 10 s, where we observe a peak at the obvious location of k = 1. It is the emergence of higher harmonics for larger strains, however, that tells when we leave the LRR. From S4 Fig, we see that γ = 0.1 is definitely not in the linear regime for the B-In-M geometry, since it has a k = 3 Fourier coefficient that is only about ten times smaller than the one for k = 1. Increasing γ to 0.15, we see an even stronger response from the k = 3 harmonic, and the emergence of the k = 2 harmonic as well.

Timescales of stress relaxation

We first seek to gain some understanding of the timescales in the system with a stress relaxation test, snapshots from which are shown in Fig 3, with full movies shown in S3 Video (for the homogeneous meshwork) and S4 Video (for the B-In-M geometry). We simulate for 0.25 seconds using ω = 1 Hz, so that the system ends at its maximum strain γ (the results are largely independent of both the strain and shear rate used to obtain it). In Fig 4, we show the decay of the stress, with t = 0 denoting the point at which the system reaches the maximum strain. In addition to our dynamic network geometries, for which we average over 30 trials, we consider permanently linked, interconnected networks of the type considered in [44], where ensures the network is well above the rigidity percolation threshold (over 95% of the fibers are in one connected component, where two fibers are connected if they are linked by at least one CL). These networks, for which we do not include fiber turnover, give smoother stress profiles, and so we only average over 5 trials.

Download:

Fig 3. Snapshots from the stress relaxation test in a homogeneous meshwork (top) and B-In-M geometry (bottom).

We begin with an unsheared unit cell at left, then shear the network until it reaches a maximum strain (20% in this case, shown at right), after which we turn off the shear and measure the relaxation of the stress. As in Fig 2, the colored fibers are in bundles, and the CLs are shown in black. For the B-In-M geometry, these snapshots are from a smaller domain (L_d = 2 μm) than we typically use so that we can also see the CLs.

https://doi.org/10.1371/journal.pcbi.1009240.g003

Download:

Fig 4. Normalized stress profiles over time in the stress relaxation test.

We consider four different systems: a permanent network (blue), for which the stress relaxes to a nonzero value (g₀ ≈ 0.07 Pa after accounting for normalization), and three dynamic networks, for which the stress relaxes to a value of zero (orange is the homogeneous meshwork, yellow the B-In-M morphology, and purple the B-In-M morphology with ten times larger viscosity). To illustrate the point that there are multiple intrinsic relaxation timescales in the system, we show a double-exponential curve which approximately matches the decay of stress for both permanent and transient CLs.

https://doi.org/10.1371/journal.pcbi.1009240.g004

Fig 4 gives us two pieces of information that are important for our analysis going forward. First, we observe that statically-linked networks can store elastic energy, since g₀ ≔ σ₂₁(t → ∞) > 0, while all of the energy dissipates in transiently-linked networks. Second, both permanent and dynamic networks have multiple relaxation timescales, as most of the stress (≈ 80% for dynamic networks) relaxes over the first 0.2 s or so, with the remaining part taking seconds to relax. In Fig 4, we show a two-timescale decay curve that approximately, but not exactly, matches the decay of the stress in both statically- and transiently-linked networks. Given that one of the time constants is on the order 0.1 s and the other is on the order 1 s, Fig 4 demonstrates that both slow and fast timescales exist in cross-linked networks. Our rheological experiments will give more precise estimates for the individual timescales, but it is worth noting that the dissipation of all of the stress over a timescale implies that the longest timescale in the system cannot be much longer than five seconds. This is not surprising since the filament turnover time is of this order, but is in sharp contrast to modeling studies without turnover, where stress takes hundreds of seconds to dissipate [31](Fig. 4).

One observation we can make from Fig 4 is that the long-time behavior is similar for all dynamic networks, so the long timescale behavior is roughly independent of viscosity. We made this latter observation once before when we saw that a similar turnover time can be used to generate a B-In-M morphology for larger viscosity systems. We will look at the short-timescale behavior more closely in our rheological tests, which we discuss next.

Viscoelastic moduli

In S5 and S6 Videos, we show animations of the oscillatory shear experiment for the homogeneous meshwork and B-In-M geometries, respectively. In Fig 5, we plot the elastic and viscous moduli as a function of frequency for the five different systems defined in Table 2. We also include the data for a permanent, interconnected network of the type considered in [44], which we show with and without the constant elastic element g₀ = G′(ω → 0) determined in Fig 4. For the system with ten times larger viscosity, we show the data with time/frequency rescaled by a factor of ten, so that the value for ω = 1 Hz is mapped to ω = 10 Hz. These data are for the local drag model, since it is faster to simulate; we will consider the effect of nonlocal hydrodynamics in a later section.

Download:

Fig 5. Elastic (left) and viscous (right) modulus for the five systems we study.

The system parameters are described in Tables 1 and 2 (see Table 2 for what parameters are varied in each system). For the B-In-M morphology with ten times the viscosity (μ = 1 Pa ⋅ s), we show the data with time rescaled by a factor of ten. We also show the results for permanently-cross-linked networks, and for the elastic modulus the dashed light blue curve is the remaining elastic modulus when g₀ ≈ 0.07 Pa (measured in Fig 4) is subtracted off.

https://doi.org/10.1371/journal.pcbi.1009240.g005

Examining the data, we make the following observations:

On very short timescales (shorter than 0.05 s, or frequencies larger than 20 Hz), the viscous moduli all approach the viscous-fluid scaling of ≈ (0.13 Pa ⋅ s) ω. The timescale on which this occurs is directly proportional to the viscosity, as the rescaled G′′ curve for the larger viscosity system makes the transition at the same time as the other B-In-M curves. In addition, the elastic moduli for B-In-M networks are about twice as large as those of homogeneous meshwork curves at high frequency, so the elastic modulus at short times is proportional to the link density.
On very long timescales (several seconds, low frequencies), the viscous modulus again scales linearly with frequency, but with a larger slope (G′′ ≈ (0.3 Pa ⋅ s) ω), which indicates that the links have become viscous on those timescales. The elastic modulus is much smaller than the viscous modulus on these timescales, and both of the moduli scale nonlinearly with the link density, as B-In-M morphologies have moduli about 4 times as large as homogeneous meshworks.
There is an intermediate timescale of about 1 second on which the moduli for the B-In-M meshwork (yellow) diverge from those for the system with twice the binding/unbinding rate (purple), and on which both of these curves diverge from the light blue permanent network curve. In particular, for ω ≤ 2 Hz the elastic and viscous modulus are both smaller for the system with faster link turnover. This timescale is related to (in fact, it is exactly equal to) 1/k_off, the timescale on which individual links bind and unbind. The rescaled curve with larger viscosity also diverges from the other B-In-M curves on this intermediate timescale, which indicates that a timescale has been introduced (by the dynamic linking) that does not scale with viscosity.

We will explore each of these fast, slow, and intermediate timescales in subsequent sections. The picture we sketch is of a rigid network for timescales shorter than the fastest intrinsic timescale τ₁ ∼ 0.1 s; on these timescales the links are predominantly elastic and the viscous modulus is the same as if the fibers were in a fluid without CLs. On timescales longer than τ₁ but shorter than the intermediate timescale τ₂ ∼ 1 s, the links still appear permanent, but they provide additional viscosity by deforming the network to its steady state in response to flow deformations. On timescales longer than τ₂, but shorter than the longest timescale τ₃ ∼ 5 s, the links come off before they are able to fully respond to the deformations induced by the background flow, and the slope of the elastic and viscous modulus changes relative to the one observed in permanent networks. Finally, on timescales longer than τ₃, the network is almost completely viscous, and the viscosity is controlled by the morphology.

Short timescales: Elastic links and viscous fibers

Our hypothesis is that the shortest timescale in the system, τ₁, dictates when the links make a negligible contribution to the viscous modulus. For timescales τ < τ₁, the network is essentially frozen: the fibers are rigid, the links are static, and the fibers and links contribute exclusively to the viscous and elastic modulus, respectively. Specifically, as ω → ∞, we expect that the links will make a constant contribution to G′ and that G′′ will scale like η₀ ω, i.e., as a viscous fluid, with the viscous coefficient η₀ independent of the number of links in the system. In Fig 5, we show that η₀ ≈ 0.13 Pa ⋅ s if ω is given in Hz. We have confirmed that this scaling continues for frequencies all the way up to 1000 Hz.

To verify that the viscous modulus at short timescales is dominated by the fibers only (independent of the number and nature of the links), we compare to the theoretical shear viscosity enhancement for an isotropic suspension of rigid cylindrical fibers, where the mobility is computed by the local drag approximation (first line in Eq (8)). For dilute suspensions, the enhanced viscosity is given to order (log ϵ)⁻² in terms of the fiber density f and aspect ratio ϵ as [86–88] (19)

For our parameters (ϵ = 0.004, ), we obtain for cylindrical fibers. Recalling that G′′ = 2πωη₀ for a viscous fluid, our scaling in Fig 5 shows that we obtain a similar value of η₀/μ = 0.13/(2π × 0.1) ≈ 0.21. Note that, to obtain Eq (19), we dropped the finite part integral in [86](Eq. (7.1)), which gives 1/2 instead of 3/2 in the denominator of [86](Eq. (8.13)). If we include intra-fiber hydrodynamics in the mobility, the 0.5 in Eq (19) changes back to 1.5. Substituting our parameters and recomputing, we get with intra-fiber hydrodynamics, which means the local drag approximation gives only 80% of the correct viscosity.

It is important to separate the viscous contribution from the links, which approaches zero at short timescales, from that due to the fibers, which is infinite as ω → ∞. For this reason, we define (20) as the viscous modulus coming from the fibers. In Fig 6a, we examine the behavior of the moduli due to the links, G′ and , rescaled by the average number of links in the system . After rescaling time in the system with larger viscosity, we see that the data for ω > 10 Hz can be collapsed onto a single curve for both the elastic and viscous modulus. Thus on timescales shorter than τ₁ ∼ 0.1 s, each link behaves as an elastic element with strength independent of morphology. Note that the elastic modulus is more than 2–3 times larger than the viscous modulus (coming from the links) at high frequencies.

Download:

Fig 6. Renormalizing the elastic and viscous modulus on short and long timescles.

We show the elastic (left) and viscous (right) modulus due to the CLs for the five systems we study. Systems are described in Tables 1 and 2 (see Table 2 for what parameters are varied in each system). To obtain a viscous modulus due to the links alone, we subtract the component due to the fibers defined in Eq (20). (a) We normalize by the link density and see that the curves all match at short timescales on the order τ₁ ≈ 0.02 s (after we also rescale time by viscosity). (b) We normalize by the link density multiplied by the bundle density. For the system with ten times larger viscosity, we also include the raw data (not rescaled) as a dotted green line. In the viscous modulus plot on the right, we show a linear slope as a dashed black line and define τ₃ ≈ 5 s as the end of the low-frequency linear regime in G′′. This timescale is roughly independent of viscosity.

https://doi.org/10.1371/journal.pcbi.1009240.g006

Long timescales: Viscosity depends on morphology and link unbinding rate

Let us now transition to the opposite limit of long timescales. At long timescales, Fig 6a shows that each system is about five times more viscous than elastic (e.g., compare the values at ω = 0.2 Hz), and that normalizing by does not place the curves on top of each other. We say that these “long timescales” are longer than τ₃, the longest intrinsic timescale in the system, or the timescale on which the network totally “remodels” itself. Our goal in this section is two-fold: first, we would like to find a new rescaling that better matches the curves at low frequencies, and second, we would like to estimate τ₃, which we do by finding the regime in which G′′ is a linear function of ω for each system.

Fig 6b shows our attempt to accomplish both of these goals. For both the viscous and elastic modulus, we rescale by the mean link density multiplied by the mean number of bundles , the idea being that bundles make a stronger contribution to the moduli at low frequency. We see a better match at low frequencies than in Fig 6a, which implies that bundles have a stronger effect on the mechanical behavior at long timescales than at short ones. In particular, the homogeneous meshworks (blue and orange) and B-In-M (yellow) curves all scale onto the same curve for frequencies ω ≤ ω₃ ≈ 0.2 Hz = 1/τ₃, but the curves with larger viscosity (green) and faster link (un)binding (purple) do not scale as well. In both of these cases, there is a change in a timescale shorter than τ₃, as we discuss in the next section.

To estimate τ₃, we determine the regime where the viscous modulus G′′ scales linearly with ω and equate τ₃ with the end of the linear region. From Fig 6b, we see that τ₃ ≈ 5 seconds, and that the linear regime endures longer for homogeneous meshworks, which implies that they start becoming purely viscous at shorter timescales than systems with bundles (or, equivalently, these networks “remodel” themselves faster).

To understand exactly how much viscosity is being provided by the links and bundles in the different systems, we use the slope in Fig 6b to obtain (21)

Since has units links per fiber and has units bundles per μm³, 0.04 in this equation has units Pa ⋅ s/μm³. To extract a viscosity, we write , so that in the B-In-M system with and μm⁻³ we obtain μ_CL/μ = (0.04 × 12.5 × 2)/(2π × 0.1) = 1.6. This is eight times more viscosity than the fibers alone (0.2) and 1.6 times the background fluid viscosity. For the homogeneous meshwork with and μm⁻³, the additional viscosity is μ_CL/μ = (0.04 × 5.75 × 1)/(2π × 0.1) = 0.4. This is twice the viscosity of the fibers but 2.5 times smaller than that of the background fluid. Comparing the two morphologies confirms some of the prior intuition that bundles embedded in meshworks provide stronger resistance to flow than homogeneous meshworks [26].

Intermediate timescales

So far, we have addressed the two extreme limits of the suspension behavior. Over short timescales, τ < τ₁, the entire network is rigid, the CLs are purely elastic, and the additional viscosity is the same as that due to the fibers alone. Over long timescales, τ > τ₃, all of the energy in the network is dissipated. Bundles break up, there is no elastic modulus, and the morphology of the network, and the amount by which links can deform it, determines the viscous behavior.

There is also an intermediate timescale τ₂ at which some of the links, but not the entire network, start to behave viscously. While we expect τ₂ ∼ 1/k_off, this could change based on the network parameters, as the same value of k_off acts differently depending on the speed of network deformation. To better understand the timescale τ₂, we compare the elastic modulus with dynamic links to the one for statically-linked networks (by starting from the same equilibrium structure and keeping the links fixed). Similar to [44], we do this for five seconds or five cycles, and measure the modulus over the last three cycles.

Fig 7a shows the results for our B-In-M structures. The solid lines show the elastic modulus with dynamic linking, and the dotted lines show the elastic modulus with permanent linking. The timescale τ₂ is when the dynamic linking elastic modulus diverges significantly from the permanent link modulus, so that the links come off before they can provide their maximum elastic response. From Fig 7a, we see that τ₂ is obviously smaller when we increase the rate of link turnover (the purple curve diverges from the dotted yellow one faster than the solid yellow one does); this is expected from τ₂ ∼ 1/k_off, which is 1 s for the yellow curve and 0.5 s for the purple one.

Download:

Fig 7. Viscoelastic behavior on medium timescales.

(a) Elastic modulus at medium frequencies. For each network type (indicated in the legend), we compare the solid line, which has the elastic modulus with dynamic links, to the dotted line, which is the elastic modulus with permanent links. The intermediate timescale τ₂≈ 0.5 s is the inverse of the frequency where the data start to diverge. We consider only bundle-in-mesh morphologies using our standard parameters (yellow), twice the link turnover rate (purple, but the static link reference is still the dotted yellow), and ten times the viscosity (green). (b) Rescaling time to get all of the viscous modulus curves on the same plot. This is the same plot as Fig 6b (right panel), but now we rescale the time for the larger viscosity green curve by a factor of 4 instead of 10 (ω → 4ω), and we rescale the time for the faster link turnover purple curve by 1.5 (ω → ω/1.5). This demonstrates that the data do not scale simply with the parameters at medium and low frequencies, when multiple timescales are involved.

https://doi.org/10.1371/journal.pcbi.1009240.g007

Morphology has no impact on the timescale τ₂, as we obtain the same characteristic time τ₂ ≈ 1 s for the homogeneous meshwork (not shown, because the results are indistinguishable) as the B-In-M morphology. Thus the morphology only appears to have a strong influence on the longest timescale τ₃.

It is harder to determine the effect of viscosity on τ₂. On the one hand, Fig 7a shows the true data for G′ vs. ω (without scaling time by μ), so it appears that the dynamic curve diverges from the permanent one at about the same timescale regardless of viscosity. On the other hand, for a fixed k_off the difference between the solid and dashed curves at ω = 1 Hz is larger for larger viscosity (green). Thus the timescale τ₂ depends on the interaction of dynamic linking and deformation in a nontrivial way.

To illustrate this complicated dependency, we show that different rescalings are required when the viscosity (dynamics) or link turnover rates change from the base parameters. We recall from Fig 6b that the larger viscosity (green) and faster link turnover (purple) B-In-M systems do not follow the same normalization as the other systems. Fig 7b shows that scaling time by a factor of 4 (and not the expected 10) puts the green larger viscosity curve onto the rest of the curves at low frequencies, while scaling time by a factor of 1.5 (and not the expected 2) puts the curve with twice the link turnover rate on top of the others as well. The reason for these particular scalings is not obvious to us, other than that they fall somewhere between the expected scaling and unity.

Role of nonlocal hydrodynamic interactions

We have seen that the network behaves very differently on short timescales, where each link behaves elastically independent of the network morphology, than on long timescales, where the system is purely viscous and the morphology exerts a significant influence on the resistance to flow. In this section, we show that nonlocal hydrodynamic interactions significantly decrease the elastic and viscous moduli in B-In-M morphologies. We show that the decrease is most significant at low frequencies, i.e., the frequencies where we already know morphology has a strong influence on the behavior. We explain the decrease in two ways: first, nonlocal hydrodynamic interactions slow down the bundling process, causing less bundles to form for a fixed turnover time, and second, they create entrainment flows that reduce the stress inside of the bundles that do form.

We consider only the two characteristic morphologies in this section without varying any other parameters from those given in Table 1. For each of the morphologies, we compute the viscoelastic moduli with nonlocal hydrodynamics, then compute them again using local drag and intra-fiber hydrodynamics. In Fig 8, we plot the error/fraction of the moduli recovered with the various mobility operators. For a homogeneous meshwork system (blue lines), we see that the elastic modulus is the same within 10–20% whether we use local drag, intra-fiber, or full hydrodynamics. At low frequencies, the viscous modulus is also the same within 10% regardless of the hydrodynamic model used, but there is an obvious error in the viscous modulus when the local drag model is used at high frequencies. This discrepancy of about 20% is also present in permanent networks [44](Fig. 9(b)) and in B-In-M morphologies, and in all cases it can be mostly recovered by adding only intra-fiber hydrodynamics to the mobility. The 20% difference between local drag and intra-fiber hydrodynamics is similar to the theoretical estimate discussed after Eq (19). There appears to be a small increase (at most 5%) in the viscous modulus when we switch from intra-fiber hydrodynamics to fully nonlocal hydrodynamics, which would be more significant at smaller mesh sizes [87].

Download:

Fig 8. Proportion of the (left) elastic and (right) viscous modulus that is recovered using various mobility approximations.

We compute the modulus using full hydrodynamics, then plot the fraction of it recovered using local drag or intra-fiber hydrodynamics. We show the results for the homogeneous meshwork in blue and the B-In-M morphology in orange. For each line color, a solid line shows the results for local drag and a dotted line shows the results for intra-fiber hydrodynamics. Intra-fiber hydrodynamics cannot explain the deviations in the elastic and viscous modulus at low frequency for the B-In-M geometry.

https://doi.org/10.1371/journal.pcbi.1009240.g008

The more interesting changes with nonlocal hydrodynamics come at low frequencies in the B-In-M geometries. Unlike with the homogeneous meshwork, the local drag model makes about a 50% overestimation in the elastic modulus and a 30% overshoot in the viscous modulus at low frequencies for B-In-M morphologies, with the largest change coming at ω = 0.5 Hz, or a timescale of 2 seconds. We have previously seen that morphology has a strong influence on the moduli at these timescales, so we explore the impact of nonlocal hydrodynamics on morphology next.

Rescaling of time cannot explain long-timescale moduli.

The first possible explanation for the decrease in G′ and G′′ with nonlocal hydrodynamics at low frequencies is a change in the network structure. In this section, we show that, while there are less bundles and links in the dynamic steady state with full (intra- and inter-fiber) hydrodynamics for a given turnover time, this by itself cannot entirely explain the decrease in the moduli.

To do this, we focus on the frequency where full hydrodynamics matters the most as a percentage. For ω = 0.5 Hz, we measure the moduli with intra-fiber hydrodynamics as Pa, Pa, while the moduli with full hydrodynamics are Pa, Pa. As shown in Table 4, the steady state link and bundle density are 12.3 and 2.2 μm⁻³, respectively, with intra-fiber hydrodynamics, while with full hydrodynamics they are are 11.7 and 2.0 μm⁻³. Thus at its dynamic steady state, the network has on average fewer links and fewer bundles when we simulate with full hydrodynamics than when we simulate with intra-fiber hydrodynamics. This is because disturbance flows in bundling fibers (fibers moving towards each other) oppose the direction of motion, which slows down the bundling process. Thus when we fix a turnover time, there are on average fewer bundles (and fewer links) when the bundling process is slower.

Download:

Table 4. Statistics for system with ω = 0.5 Hz and varying turnover times and mobility models.

https://doi.org/10.1371/journal.pcbi.1009240.t004

To compensate for the changes in structure, we drop the turnover time or increase the viscosity for intra-fiber hydrodynamics simulations so that the steady state morphology matches the morphology with hydrodynamics as closely as possible. We then measure the moduli for these new steady states. In Table 4, we show our attempt to match the statistics for intra-fiber hydrodynamics with smaller turnover times with those for hydrodynamics with τ_f = 10 s. Comparing τ_f = 10 s with full hydrodynamics with τ_f = 9.3 s with intra-fiber (IF) hydrodynamics, we see the IF simulations have larger moduli, even when we closely match the steady state link density and bundle density . The same is true when we attempt to rescale time by increasing the background fluid viscosity (second to last line in Table 4), as in that case the larger values of G′ and G′′ persist despite a decrease in link and bundle density. So nonlocal hydrodynamics must both change the bundle morphology and lower the stress for a fixed morphology.

Bundled fibers: When flow reduces stress.

The reason for the decrease in moduli with inter-fiber hydrodynamics for a fixed morphology has to do with the flows inside a bundle. If the fibers are packed tightly within a bundle, then we expect their disturbance flows to have a strong influence on each other. When a fiber in a bundle moves with the bulk fluid, the disturbance flow it creates is in the same direction as the bulk motion, so the other fibers in the bundle will naturally move as well. This effect, which is not present when we use local drag or intra-fiber mobility, explains the reduction in stress for a fixed rate of strain.

Fig 9 establishes this more rigorously. We consider a set of nine fibers with L = 1 μm without cross-links. The fibers are arranged in a regular octagon with side length ℓ with another fiber centered at the origin. We apply a constant straining flow with strength s⁻¹ and measure the stress over the first second in Fig 9. The solid orange line shows the stress with full hydrodynamics and the simulation parameter of ℓ = 0.05 μm, so that the fibers are very close. In this case, we see that the stress with full hydrodynamics is slightly more than half of that with local drag (which is independent of ℓ). Similar to our simulation result for cross-linked gels, the stress decrease is not explained by intra-fiber hydrodynamics, which increases the stress from local drag by the theoretically predicted 20–25%. In fact, as discussed in below Eq (19) and in [44](Sec. 6.3.2), disturbance flows from intra-fiber hydrodynamics (which try to stretch the fiber) make the constraint forces on each fiber larger, which causes an increase in stress. Thus the decrease in stress with nonlocal hydrodynamics comes from the nonlocal flows induced by the fibers on each other. We demonstrate this in Fig 9 by spacing the fibers farther apart (ℓ = 0.2 μm) and showing that the stress with full hydrodynamics increases towards the intra-fiber curve.

Download:

Fig 9. Reduction of stress in bundles explains smaller moduli with hydrodynamics.

(Left) We manufacture a bundle geometry without CLs by placing nine fibers of length L = 1 μm (red) in and around an octagon with side length ℓ and straining with constant rate s⁻¹ until t = 1 s (blue). (Right) The resulting stress evolution for different hydrodynamic models. The local drag (blue) and intra-fiber (yellow) results are independent of ℓ, while the stress for full hydrodynamics (orange) depends strongly on ℓ. For ℓ = 0.05 μm (solid orange, the simulation parameters), there is a significant decrease in stress which comes from the entrainment of the fibers in each other’s flow fields. For ℓ = 0.20 μm (dashed orange), the decrease is minimal; note that for full hydrodynamics with ℓ ≫ L we would recover the “intra-fiber” curve.

https://doi.org/10.1371/journal.pcbi.1009240.g009

A continuum framework: The generalized Maxwell model

In order to enable whole-cell modeling, it is important to introduce a continuum model for the passively cross-linked actin gel. We are motivated by the behavior of the network: at short timescales, it scales as a viscous fluid with a constant elastic modulus (a dashpot in parallel with a spring), while at long timescales it is purely viscous. Combining these two, we introduce a generalized Maxwell model of the type shown in Fig 10. We have the viscous dashpot of strength η₀ in parallel with three Maxwell elements, each of which has a strength g_i and an associated relaxation timescale τ_i, on which the element goes from elastic (for timescales shorter than τ_i) to viscous (longer than τ_i). The total elastic and viscous (in excess of the solvent viscosity) modulus are given by sums of the viscous dashpot and each Maxwell element [80] (22) where we have denoted the 7 parameters by a vector p. As discussed in [89, 90], we can obtain the fitting parameters (g_i, τ_i, η₀) by maximizing the log-likelihood function (23) where and are the uncertainties in G′(ω_k) and G′′(ω_k), respectively, and K is the number of frequencies studied. Weighting the observations by uncertainty will cause the fit to give more weight to more certain measurements. The uncertainty in the fit is estimated by the square root of the diagonal entries of the inverse of the Fisher information matrix, (24) so that p_i±δp_i is a 95% confidence interval for p_i.

Download:

Fig 10. Our continuum model, informed by the timescales discussed in previous sections.

We use three Maxwell elements with timescales τ₁, τ₂, and τ₃, all in parallel with a viscous dashpot to describe the network. The viscous dashpot η₀ represents the high frequency viscosity of the permanently cross-linked fiber suspension. The first Maxwell element has timescale τ₁ ≈ 0.02 seconds associated with it, and represents the relaxation of the fibers to a transient elastic equilibrium (the networks before and after relaxation are shown to the left and right of this Maxwell element; the relaxing fibers are shown in blue); on this timescale, the links are effectively static. The second Maxwell element, with timescale τ₂ ≈ 0.5 s, represents the unbinding of some links (shown more transparent than the others) and the appearance of new links (orange) − compare the networks to the left and right of this Maxwell element. The third Maxwell element with timescale τ₃ ≈ 5 s represents network remodeling (compare the networks to the left and right of this element); for timescales larger than τ₃, some of the fibers (shown in green) and links (orange) turn over and the network completely remodels from the initial state.

https://doi.org/10.1371/journal.pcbi.1009240.g010

In S5 Fig, we verify that using three (as opposed to two or four) Maxwell elements is indeed the best choice for the bundle-embedded meshwork with τ_f = 10 s and full hydrodynamics (the other systems are similar). As detailed in [90, 91], we determine this by increasing the number of Maxwell elements until the fit stops improving substantially and we start to overfit the data, leading to ill-conditioning. In addition to comparing the fit with three timescales to that with two or four, in S5 Fig we also compare to a fit using a continuum of timescales, which is suggested by physical theories [14, 30] and used in prior studies [7]. In this approach, the discrete values of g_i in Eq (22) are replaced by a continuous function g(τ), which is integrated (instead of summed) on 0 ≤ τ ≤ τ_max. If we take the common choice of g(τ) = g₀τ^−α [92], we obtain (25) where (g₀, α, τ_max, η₀) are the fitting parameters. In S5 Fig, we show that the fit to the data with this continuum spectrum of timescales is not substantially better than the three timescale generalized Maxwell fit. Our preference is the three timescale fit, since in this case we can assign physical to meaning to each of the timescales, as we did in earlier sections. Given that adding more discrete timescales leads to ill-conditioning, the continuum approach, where there are infinitely many timescales (and in fact infinitely many for choices g(τ)), is even more poorly conditioned, and in that case it is difficult to assign a physical meaning to the fitting parameters.

The three-timescale fit that we obtain for the B-In-M system with full hydrodynamics, along with the contribution of each of the separate Maxwell elements, is shown in Fig 11. Admittedly, the elastic modulus data G′(ω) at small frequencies do not fit the generalized Maxwell model since they do not decay as ω². This suggests that the material response is likely more complicated than just three Maxwell elements. At the same time, we have already shown that G′′(ω) scales linearly with ω at small frequencies, and thus that the viscous modulus, which dominates the elastic modulus at low frequencies, does fit the generalized Maxwell assumption. Our choice is therefore to fit only G′′, and not G′, for frequencies less than 0.1 Hz (10 s timescales).

Download:

Fig 11. Fitting the model of Fig 10 (the moduli from Eq (22)) to the data from the B-In-M system.

The data are shown in blue, with the total fit shown in red. The dotted lines show the contribution of each element in Fig 10 to the total fit. Here τ₁ = 0.04 s is the fastest timescale shown in purple, τ₂ = 0.41 s is the intermediate timescale shown in green, τ₃ = 4.5 s is the longest timescale shown in light blue. The yellow dotted line shows the contribution of the pure viscous element.

https://doi.org/10.1371/journal.pcbi.1009240.g011

Table 5 gives the fitting parameters for all of the systems we have considered with multiple options for the mobility. For almost all systems, we have τ₁ ≈ 0.03 s, τ₂ ≈ 0.3 s, and τ₃ ≈ 3 s, which is in line with our estimates of these timescales in Figs 6 and 7a. The only exception is the system with higher viscosity; in this case we have already shown that the timescale τ₁ scales with viscosity to become about 0.3 s, while the timescale τ₂ remains about 0.5 s. The result is that τ₁ and τ₂ blend into a single timescale, and a two-timescale model is the best choice for the data, as we show in S6 Fig. Comparing B-In-M and meshwork geometries, we observe that the long timescale τ₃ is either longer (comparing B-In-M and meshwork with full hydrodynamics) or its contribution in Eq (22) is stronger (comparing B-In-M and meshwork with local drag), which is consistent with our observation that B-In-M geometries have longer remodeling times and a higher resistance to deformation. Finally, we notice how the viscosity η₀ scales directly with μ and increases by ≈ 30% when we account for intra-fiber or nonlocal hydrodynamics, in accordance with our observations in previous sections.

Download:

Table 5. Parameters in the three-timescale generalized Maxwell model shown in Fig 10.

https://doi.org/10.1371/journal.pcbi.1009240.t005

Conclusion

We experimented numerically with dynamic actin gels made of micron-long actin fibers at concentrations characteristic of in vitro experiments. In the simulations, the fibers turned over rapidly on a 5 to 10 second scale, similar to some experimental measurements [46, 47]. The gel was crosslinked dynamically with flexible spring-like cross linkers (CLs) characterized by mechanical and kinetic parameters similar to those for α-actinin, a ubiquitous actin CL. We observed that gels with faster fiber turnover evolved into relatively homogeneous meshworks, while gels with slower turnovers morphed into bundle-in-mesh (B-in-M) networks with small actin bundles immersed into cross-linked meshworks. This result agrees with experimental observations [93].

Our simulations revealed three principal time scales characterizing the mechanics of dynamic cross-linked actin gels. On the fastest time scale, ∼0.03 seconds, actin fibers relax viscously, locally, and rapidly to the transient elastic equilibrium generated by the network configuration. On intermediate timescales, ∼0.5 seconds, the CLs turn over, generating new transient elastic equilibria, and finally on the slow time scale, ∼3–5 seconds, the network undergoes global remodeling through a combination of filament and CL turnover. One of the most useful practical results of the simulations is that very complex cross-linked gel mechanics can be approximated with the generalized Maxwell mechanical circuit shown in Fig 10, in which three Maxwell elements in parallel correspond to elastoviscous gel deformations on the three characteristic time scales, in addition to the effective viscosity of the actin fiber suspension in the background fluid.

Our simulations predict that at small frequencies (time scales longer than a few seconds), the effective elastic and viscous behaviors are controlled by the CL mechanics (as opposed to fiber mechanics). In this regime, an effectively viscous mechanical response to deformation dominates relatively weak elastic behavior. At low frequencies, the dominant viscous response originates from the CLs stretching and deforming the fibers (in addition to the significant viscosity of the background fluid); the fibers by themselves contribute little to the net viscosity. The mechanical moduli scale nonlinearly with the CL density, because the CLs within actin bundles respond to shear differently than the CLs in the actin mesh. Thus, the gel’s morphology exerts a significant influence on the actin resistance to flow at long timescales.

At high frequencies, or timescales shorter than one tenth of a second, we found that the effective elastic modulus of the network is higher than (but on the same order of magnitude as) the effective viscous modulus (once we removed the viscous scaling generated by the fibers). The elasticity simply originates from multiple elastic springs of CLs that can be considered static on the short time scale. The greatest contribution to the net viscosity is that of the background fluid, with the fiber suspension contributing about 25% additional viscosity. At moderate frequencies around 1 Hz (intermediate timescale, between 0.1 and 3 seconds), the mechanical behavior of the crosslinked network is complex, combining comparable contributions of all mechanical factors described for low and high frequencies.

Our model allowed us, for the first time, to estimate quantitatively the role of nonlocal inter-fiber hydrodynamic interactions on the actin gel rheology. We demonstrated that, while these interactions have little effect in a homogeneous cross-linked actin meshwork, they significantly decrease the elastic and viscous moduli of heterogeneous B-In-M gels. We showed that this decrease is most significant at low frequencies, i.e., on timescales of a few seconds, when the network morphology has a strong influence on its mechanical behavior. At these frequencies, in the presence of bundles, the hydrodynamic interactions cause 30–50% downward corrections to the viscoelastic moduli. Two factors are responsible for this effect: first, hydrodynamic interactions slow down the bundling process, so fewer bundles form, and second, nonlocal hydrodynamics creates entrainment flows that reduce the stress inside of the bundles that do form.

While almost all of our simulations fell in the linear viscoelastic regime, we did observe a strain hardening effect, with the network resistance increasing beyond the linear response for higher-amplitude deformations. This effect has been observed before experimentally [9], which has prompted a variety of explanations and computational studies. Both Mulla et al. [14] and Kim et al. [94], for instance, propose that strain hardening is due to the mechanics of individual filaments, which stiffen nonlinearly in response to applied stress. Because our goal was to measure the viscoelastic moduli of the networks in their dynamic steady state, we did not systematically explore the viscoelastic behavior at large strains or perform measurements at zero frequency. That said, our simulations do offer two interesting insights: first, the mechanical behavior becomes nonlinear at just 10% or greater strains, which is a typical result for in vitro networks [95], and second, there are more bundles (which are more resistive to strain) at larger strain rates, and so at least one reason for the nonlinearity is the mechanosensitive character of the actin gel morphology. Our conclusion is identical to the modeling study [96], which found that strain combined with dynamic cross-linking induces signficant morphological changes in the network, leading to strain hardening [96](Fig. 3). In fact, it was recently shown that actin filament mobility is a necessary condition for strain hardening and mechanical hysteresis [97].

In the linear response regime, our results are in line with most of the prior experimental and modeling studies. We found viscoelastic moduli on the order 0.1–10 Pa [10, 23], with the links becoming viscous on long timescales and elastic on short ones [7]. Our measurements showed that decreasing the fiber bending stiffness decreased the viscoelastic moduli, but by a slower than linear rate [66] (in our study, decreasing κ by a factor of ten reduced the moduli by at most 30%). We also found that the elastic stress in the gel (from the CLs exerting force on the fibers) is much larger than the viscous drag stress over long times, which agrees with a previous computational study [94].

We failed, however, to observe a local maximum or minimum in the viscous modulus at intermediate frequencies, as was reported in [15, 43]. One explanation for this is a difference in parameters. In the computational study [43], the unloaded CL unbinding rate is 0.1 s⁻¹, and the CL stiffness is 10 pN/nm (10,000 pN/μm). In this parameter regime, the CLs always make the most of their mechanical contribution to the network mechanics while bound, since they cross-link the fibers for seconds and are extremely stiff. Indeed, [43](Fig. 6) shows that the peak in G′′, when the timescale of the driving frequency matches the link unbinding rate, becomes less sharp and almost disappears when k_off increases to the order of one second, as it is in our study here. This suggests that the local maximum is only observed when the relaxation of the network is many orders of magnitude faster than the link binding rate. This issue clearly requires more investigation, and in the future we will explore how using a force-dependent unbinding rate, as is done in [43], impacts the results, and whether this assumption is responsible for the peak.

What all modeling studies, including ours, suffer from is that only few features of the actin gel are simulated, with many factors, forces, and processes either ignored or approximated crudely. There is, of course, a good reason for this: even a very limited caricature of the gel exhibits complex dynamics; the fact that many experimental studies paint very different pictures of the gel mechanics likely stems from the same biological complexity. We propose that to understand the full complexity, one has to keep adding dynamic processes to previous simpler benchmark models and examine the changes that new dynamics bring. Thus, our next step will be to add thermal forces responsible for random bending of the fibers. For a one micron long fiber with persistence length close to 20 microns, the thermal bending causes lateral fluctuations of no more than 30 to 40 nm, which is significantly less than the mesh size of the networks we investigated, so we do not expect that our results will change drastically. Indeed, a previous computational study in permanent networks showed that thermal fluctuations only impact the mechanical moduli when the average distance between cross links is greater than the persistence length of the fibers [66](Fig. 7(c)), which is a regime we are far outside of here.

Despite this speculation, the fact remains that our study was unable to reproduce a number of experimental observations in transiently cross-linked actin networks that have been traced to the underlying thermal fluctuations of the filaments. For instance, we report a viscous modulus that decays as for a viscous fluid (G′′ ∼ ω at low frequencies), whereas a number of recent experimental, modeling, and theoretical studies [13, 16, 30, 31] report an ω^1/2 dependence of both moduli at low frequencies. It is important to note, however, that fiber turnover, which was not accounted for in these previous studies, changes the behavior on timescales larger than the turnover time, as does our use of shorter filaments relative to [31]. In the high-frequency range, our results show viscous-fluid scaling of G′′ at large ω, but experiments [9] and models [43] have shown an ω^3/4 dependence at high frequencies. Now that we have performed a benchmark study which shows that a non-Brownian cross-linked actin network exhibits this behavior, it will be interesting to see if adding thermal bending modes [30, 31] reproduces the scaling relations observed previously, and how these relations are affected by hydrodynamic interactions. In concert with this, we will consider more bundled morphologies, where the flexibility of the filaments, and therefore their transverse bending fluctuations, may become more important. Similarly, we will test how explicit simulation of translational and rotational diffusion and (un)binding of individual CLs [31], as well as force dependence of CL kinetics, affect these results.

Other mechanical features of actin that will be interesting to explore include elastic twisting of actin filaments, which have been posited to contribute to elastic energy storage [85] and to the emergence of chirality in cross-linked gels [98]. Last, but not least, in vivo, the gels are very heterogeneous in the sense that there is a distribution of fiber lengths and turnover times, and a diversity of CL types co-existing with myosin motors [70, 99, 100]. How network heterogeneity and active molecular force generation affect the gel mechanics is another important question.

Supporting information

S1 Fig. Dynamic steady state for τ_f = 5 seconds.

Number of links per fiber (top left), bundle density (top right, number of bundles / ), % of fibers in bundles (bottom left), and maximum bundle size (bottom right) for simulations from an initialized isotropic state run to the dynamic steady state for τ_f = 5 seconds using three different system sizes: L_d = 2 (200 fibers, blue), L_d = 3 (675 fibers, red), and L_d = 4 (1600 fibers, yellow). Error bars are two standard errors over five trials.

https://doi.org/10.1371/journal.pcbi.1009240.s001

(EPS)

S2 Fig. Dynamic steady state for τ_f = 10 seconds.

Number of links per fiber (top left), bundle density (top right, number of bundles / ), % of fibers in bundles (bottom left), and maximum bundle size (bottom right) for simulations from an initialized isotropic state run to the dynamic steady state for τ_f = 10 seconds using three different system sizes: L_d = 2 (200 fibers, blue), L_d = 3 (675 fibers, red), and L_d = 4 (1600 fibers, yellow). Error bars are two standard errors over five trials.

https://doi.org/10.1371/journal.pcbi.1009240.s002

(EPS)

S3 Fig. Homogeneous meshwork for smaller bending stiffness κ.

Side (left) and top (right) view of the homogeneous meshwork with κ = 0.007 pN ⋅ μm², which is 1/10 of the stiffness of actin. Despite the smaller bending stiffness, the fibers still appear relatively straight, with the exception of some of the fibers in bundles.

https://doi.org/10.1371/journal.pcbi.1009240.s003

(TIFF)

S4 Fig. Spectrum of the stress for the B-In-M network.

The coefficients are defined in (18) in the main text. We consider various maximum strain values γ = 0.01 (green), 0.025 (purple), 0.05 (yellow), 0.1 (orange), and 0.15 (blue). We see the emergence of a k = 3 nonlinear harmonic for larger strains, suggesting a nonlinear response.

https://doi.org/10.1371/journal.pcbi.1009240.s004

(EPS)

S5 Fig. Maxwell model fit for the B-In-M system with full hydrodynamics and a varying number of elements.

We consider 1 (red), 2 (yellow), 3 (purple), and 4 (green) Maxwell modes in parallel with a viscous dashpot. The elastic modulus is only fit for ω ≥ 0.1 Hz. The B-In-M system with μ = 0.1 Pa ⋅ s is best fit with three elements. The dashed line shows to a continuous spectrum of timescales, (25) in the main text, with g₀ = 0.22, α = 1.4, τ_max = 3.5, and η₀ = 0.16, which is not a better fit than the 3 timescale generalized Maxwell model.

https://doi.org/10.1371/journal.pcbi.1009240.s005

(EPS)

S6 Fig. Maxwell model fit for the system with ten times larger viscosity and local drag mobility.

We consider 1 (red), 2 (yellow), 3 (purple), and 4 (green) Maxwell modes in parallel with a viscous dashpot. The elastic modulus is only fit for ω ≥ 0.1 Hz. This higher viscosity system (μ = 1 Pa ⋅ s) is best fit with only two elements, since the three and four element Maxwell models have a timescale larger than 100 seconds, indicating overfitting.

https://doi.org/10.1371/journal.pcbi.1009240.s006

(EPS)

S1 Video. Homogeneous meshwork system with L_d = 2 μm in its dynamic steady state.

https://doi.org/10.1371/journal.pcbi.1009240.s007

(MP4)

S2 Video. Bundle-in-mesh geometry with L_d = 3 μm in its dynamic steady state.

https://doi.org/10.1371/journal.pcbi.1009240.s008

(MP4)

S3 Video. Stress relaxation test for the homogeneous meshwork geometry.

This shows five seconds of simulation time in a domain with L_d = 2 μm.

https://doi.org/10.1371/journal.pcbi.1009240.s009

(MP4)

S4 Video. Stress relaxation test for the bundle-in-mesh geometry.

This shows five seconds of simulation time in a domain with L_d = 2 μm, which is smaller than the L_d = 3 μm domain that we use for data collection in this paper.

https://doi.org/10.1371/journal.pcbi.1009240.s010

(MP4)

S5 Video. Oscillatory shear rheology experiment for the homogeneous meshwork geometry.

This shows 5 seconds of simulation time with frequency ω = 1 Hz on a domain of L_d = 2 μm microns with maximum strain γ = 0.2. This strain is well outside the linear regime, and is used for visualization purposes only.

https://doi.org/10.1371/journal.pcbi.1009240.s011

(MP4)

S6 Video. Oscillatory shear rheology experiment for the bundle-in-mesh geometry.

This shows 5 seconds of simulation time with frequency ω = 1 Hz on a domain of L_d = 2 μm microns with maximum strain γ = 0.2. In addition to a larger strain, we also use a smaller domain here for visualization purposes (typically we use L_d = 3 μm for the B-In-M geometry).

https://doi.org/10.1371/journal.pcbi.1009240.s012

(MP4)

Acknowledgments

We thank the NYU High Performance Computing Center for access to the Greene supercomputer cluster, which we used for the simulations in this paper.

References

1. Mofrad MR. Rheology of the cytoskeleton. Annual Review of Fluid Mechanics. 2009;41:433–453.
- View Article
- Google Scholar
2. Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P. Molecular biology of the cell. Garland Science; 2002.
3. Hinner B, Tempel M, Sackmann E, Kroy K, Frey E. Entanglement, elasticity, and viscous relaxation of actin solutions. Physical Review Letters. 1998;81(12):2614.
- View Article
- Google Scholar
4. Janmey PA, Hvidt S, Käs J, Lerche D, Maggs A, Sackmann E, et al. The mechanical properties of actin gels. Elastic modulus and filament motions. Journal of Biological Chemistry. 1994;269(51):32503–32513. pmid:7798252
- View Article
- PubMed/NCBI
- Google Scholar
5. Wachsstock DH, Schwartz W, Pollard TD. Affinity of alpha-actinin for actin determines the structure and mechanical properties of actin filament gels. Biophysical journal. 1993;65(1):205.
- View Article
- Google Scholar
6. Lieleg O, Claessens MMAE, Luan Y, Bausch A. Transient binding and dissipation in cross-linked actin networks. Physical review letters. 2008;101(10):108101. pmid:18851260
- View Article
- PubMed/NCBI
- Google Scholar
7. Chaubet L, Chaudhary AR, Heris HK, Ehrlicher AJ, Hendricks AG. Dynamic actin cross-linking governs the cytoplasm’s transition to fluid-like behavior. Molecular biology of the cell. 2020;31(16):1744–1752. pmid:32579489
- View Article
- PubMed/NCBI
- Google Scholar
8. Wachsstock DH, Schwarz W, Pollard T. Cross-linker dynamics determine the mechanical properties of actin gels. Biophysical journal. 1994;66(3):801–809. pmid:8011912
- View Article
- PubMed/NCBI
- Google Scholar
9. Gardel M, Shin JH, MacKintosh F, Mahadevan L, Matsudaira P, Weitz D. Elastic behavior of cross-linked and bundled actin networks. Science. 2004;304(5675):1301–1305. pmid:15166374
- View Article
- PubMed/NCBI
- Google Scholar
10. Gardel M, Shin JH, MacKintosh F, Mahadevan L, Matsudaira P, Weitz D. Scaling of F-actin network rheology to probe single filament elasticity and dynamics. Physical review letters. 2004;93(18):188102. pmid:15525211
- View Article
- PubMed/NCBI
- Google Scholar
11. Müller O, Gaub H, Bärmann M, Sackmann E. Viscoelastic moduli of sterically and chemically cross-linked actin networks in the dilute to semidilute regime: measurements by oscillating disk rheometer. Macromolecules. 1991;24(11):3111–3120.
- View Article
- Google Scholar
12. Ruddies R, Goldmann W, Isenberg G, Sackmann E. The viscoelasticity of entangled actin networks: the influence of defects and modulation by talin and vinculin. European biophysics journal. 1993;22(5):309–321. pmid:8112218
- View Article
- PubMed/NCBI
- Google Scholar
13. Yao NY, Becker DJ, Broedersz CP, Depken M, MacKintosh FC, Pollak MR, et al. Nonlinear viscoelasticity of actin transiently cross-linked with mutant α-actinin-4. Journal of molecular biology. 2011;411(5):1062–1071. pmid:21762701
- View Article
- PubMed/NCBI
- Google Scholar
14. Mulla Y, MacKintosh F, Koenderink GH. Origin of slow stress relaxation in the cytoskeleton. Physical review letters. 2019;122(21):218102. pmid:31283330
- View Article
- PubMed/NCBI
- Google Scholar
15. Lieleg O, Schmoller K, Claessens MMAE, Bausch A. Cytoskeletal polymer networks: viscoelastic properties are determined by the microscopic interaction potential of cross-links. Biophysical journal. 2009;96(11):4725–4732. pmid:19486695
- View Article
- PubMed/NCBI
- Google Scholar
16. Yao NY, Broedersz CP, Depken M, Becker DJ, Pollak MR, MacKintosh FC, et al. Stress-enhanced gelation: A dynamic nonlinearity of elasticity. Physical review letters. 2013;110(1):018103. pmid:23383843
- View Article
- PubMed/NCBI
- Google Scholar
17. Morse DC. Viscoelasticity of concentrated isotropic solutions of semiflexible polymers. 1. Model and stress tensor. Macromolecules. 1998;31(20):7030–7043.
- View Article
- Google Scholar
18. Fabry B, Maksym GN, Butler JP, Glogauer M, Navajas D, Fredberg JJ. Scaling the microrheology of living cells. Physical review letters. 2001;87(14):148102. pmid:11580676
- View Article
- PubMed/NCBI
- Google Scholar
19. Fabry B, Maksym GN, Butler JP, Glogauer M, Navajas D, Taback NA, et al. Time scale and other invariants of integrative mechanical behavior in living cells. Physical Review E. 2003;68(4):041914. pmid:14682980
- View Article
- PubMed/NCBI
- Google Scholar
20. Ziemann F, Rädler J, Sackmann E. Local measurements of viscoelastic moduli of entangled actin networks using an oscillating magnetic bead micro-rheometer. Biophysical journal. 1994;66(6):2210–2216. pmid:8075354
- View Article
- PubMed/NCBI
- Google Scholar
21. Rouse PE Jr. A theory of the linear viscoelastic properties of dilute solutions of coiling polymers. The Journal of Chemical Physics. 1953;21(7):1272–1280.
- View Article
- Google Scholar
22. Freedman SL, Suarez C, Winkelman JD, Kovar DR, Voth GA, Dinner AR, et al. Mechanical and kinetic factors drive sorting of F-actin cross-linkers on bundles. Proceedings of the National Academy of Sciences. 2019;116(33):16192–16197. pmid:31346091
- View Article
- PubMed/NCBI
- Google Scholar
23. Kasza K, Broedersz C, Koenderink G, Lin Y, Messner W, Millman E, et al. Actin filament length tunes elasticity of flexibly cross-linked actin networks. Biophysical journal. 2010;99(4):1091–1100. pmid:20712992
- View Article
- PubMed/NCBI
- Google Scholar
24. Wagner B, Tharmann R, Haase I, Fischer M, Bausch AR. Cytoskeletal polymer networks: the molecular structure of cross-linkers determines macroscopic properties. Proceedings of the National Academy of Sciences. 2006;103(38):13974–13978. pmid:16963567
- View Article
- PubMed/NCBI
- Google Scholar
25. Tseng Y, An KM, Wirtz D. Microheterogeneity controls the rate of gelation of actin filament networks. Journal of Biological Chemistry. 2002;277(20):18143–18150. pmid:11889122
- View Article
- PubMed/NCBI
- Google Scholar
26. Tseng Y, Schafer BW, Almo SC, Wirtz D. Functional synergy of actin filament cross-linking proteins. Journal of Biological Chemistry. 2002;277(28):25609–25616. pmid:12006593
- View Article
- PubMed/NCBI
- Google Scholar
27. MacKintosh F, Käs J, Janmey P. Elasticity of semiflexible biopolymer networks. Physical review letters. 1995;75(24):4425. pmid:10059905
- View Article
- PubMed/NCBI
- Google Scholar
28. Isambert H, Maggs A. Dynamics and rheology of actin solutions. Macromolecules. 1996;29(3):1036–1040.
- View Article
- Google Scholar
29. Pegoraro AF, Janmey P, Weitz DA. Mechanical properties of the cytoskeleton and cells. Cold Spring Harbor perspectives in biology. 2017;9(11):a022038. pmid:29092896
- View Article
- PubMed/NCBI
- Google Scholar
30. Broedersz CP, Depken M, Yao NY, Pollak MR, Weitz DA, MacKintosh FC. Cross-link-governed dynamics of biopolymer networks. Physical review letters. 2010;105(23):238101. pmid:21231506
- View Article
- PubMed/NCBI
- Google Scholar
31. Müller KW, Bruinsma RF, Lieleg O, Bausch AR, Wall WA, Levine AJ. Rheology of semiflexible bundle networks with transient linkers. Physical review letters. 2014;112(23):238102. pmid:24972229
- View Article
- PubMed/NCBI
- Google Scholar
32. Peskin CS. The immersed boundary method. Acta Numer. 2002;11:479–517.
- View Article
- Google Scholar
33. Karcher H, Lammerding J, Huang H, Lee RT, Kamm RD, Kaazempur-Mofrad MR. A three-dimensional viscoelastic model for cell deformation with experimental verification. Biophysical journal. 2003;85(5):3336–3349. pmid:14581235
- View Article
- PubMed/NCBI
- Google Scholar
34. Mijailovich SM, Kojic M, Zivkovic M, Fabry B, Fredberg JJ. A finite element model of cell deformation during magnetic bead twisting. Journal of Applied Physiology. 2002;93(4):1429–1436. pmid:12235044
- View Article
- PubMed/NCBI
- Google Scholar
35. Strychalski W, Copos CA, Lewis OL, Guy RD. A poroelastic immersed boundary method with applications to cell biology. Journal of Computational Physics. 2015;282:77–97.
- View Article
- Google Scholar
36. Copos CA, Guy RD. A porous viscoelastic model for the cell cytoskeleton. The ANZIAM Journal. 2018;59(4):472–498.
- View Article
- Google Scholar
37. Guy RD, Nakagaki T, Wright GB. Flow-induced channel formation in the cytoplasm of motile cells. Physical Review E. 2011;84(1):016310. pmid:21867307
- View Article
- PubMed/NCBI
- Google Scholar
38. Desprat N, Richert A, Simeon J, Asnacios A. Creep function of a single living cell. Biophysical journal. 2005;88(3):2224–2233. pmid:15596508
- View Article
- PubMed/NCBI
- Google Scholar
39. Head DA, Levine AJ, MacKintosh F. Deformation of cross-linked semiflexible polymer networks. Physical review letters. 2003;91(10):108102. pmid:14525510
- View Article
- PubMed/NCBI
- Google Scholar
40. Kim T, Hwang W, Kamm R. Computational analysis of a cross-linked actin-like network. Experimental Mechanics. 2009;49(1):91–104.
- View Article
- Google Scholar
41. Popov K, Komianos J, Papoian GA. MEDYAN: Mechanochemical simulations of contraction and polarity alignment in actomyosin networks. PLoS computational biology. 2016;12(4):e1004877. pmid:27120189
- View Article
- PubMed/NCBI
- Google Scholar
42. Li X, Ni Q, He X, Kong J, Lim SM, Papoian GA, et al. Tensile force-induced cytoskeletal remodeling: Mechanics before chemistry. PLoS computational biology. 2020;16(6):e1007693. pmid:32520928
- View Article
- PubMed/NCBI
- Google Scholar
43. Wei X, Fang C, Gong B, Yao J, Qian J, Lin Y. Viscoelasticity of 3D actin networks dictated by the mechanochemical characteristics of cross-linkers. Soft Matter. 2021;. pmid:33646227
- View Article
- PubMed/NCBI
- Google Scholar
44. Maxian O, Mogilner A, Donev A. Integral-based spectral method for inextensible slender fibers in Stokes flow. Physical Review Fluids. 2021;6(1):014102.
- View Article
- Google Scholar
45. Lieleg O, Claessens MM, Bausch AR. Structure and dynamics of cross-linked actin networks. Soft Matter. 2010;6(2):218–225.
- View Article
- Google Scholar
46. Guha M, Zhou M, Wang Yl. Cortical actin turnover during cytokinesis requires myosin II. Current biology. 2005;15(8):732–736. pmid:15854905
- View Article
- PubMed/NCBI
- Google Scholar
47. Theriot JA, Mitchison TJ. Actin microfilament dynamics in locomoting cells. Nature. 1991;352(6331):126–131. pmid:2067574
- View Article
- PubMed/NCBI
- Google Scholar
48. Tornberg AK, Shelley MJ. Simulating the dynamics and interactions of flexible fibers in Stokes flows. Journal of Computational Physics. 2004;196(1):8–40.
- View Article
- Google Scholar
49. Aurentz JL, Trefethen LN. Block operators and spectral discretizations. SIAM Review. 2017;59(2):423–446.
- View Article
- Google Scholar
50. Driscoll TA, Hale N. Rectangular spectral collocation. IMA Journal of Numerical Analysis. 2015;36(1):108–132.
- View Article
- Google Scholar
51. Zero K, Pecora R. Rotational and translational diffusion in semidilute solutions of rigid-rod macromolecules. Macromolecules. 1982;15(1):87–93.
- View Article
- Google Scholar
52. Grazi E. What is the diameter of the actin filament? FEBS letters. 1997;405(3):249–252. pmid:9108298
- View Article
- PubMed/NCBI
- Google Scholar
53. McGrath JL, Osborn EA, Tardy YS, Dewey CF, Hartwig JH. Regulation of the actin cycle in vivo by actin filament severing. Proceedings of the National Academy of Sciences. 2000;97(12):6532–6537. pmid:10823888
- View Article
- PubMed/NCBI
- Google Scholar
54. Yoshida A, Sakai N, Uekusa Y, Deguchi K, Gilmore JL, Kumeta M, et al. Probing in vivo dynamics of mitochondria and cortical actin networks using high-speed atomic force/fluorescence microscopy. Genes to Cells. 2015;20(2):85–94. pmid:25440894
- View Article
- PubMed/NCBI
- Google Scholar
55. Zhang Y, Yoshida A, Sakai N, Uekusa Y, Kumeta M, Yoshimura SH. In vivo dynamics of the cortical actin network revealed by fast-scanning atomic force microscopy. Microscopy. 2017;66(4):272–282. pmid:28531263
- View Article
- PubMed/NCBI
- Google Scholar
56. Luby-Phelps K. Cytoarchitecture and physical properties of cytoplasm: volume, viscosity, diffusion, intracellular surface area. In: International review of cytology. vol. 192. Elsevier; 1999. p. 189–221.
57. Gittes F, Mickey B, Nettleton J, Howard J. Flexural rigidity of microtubules and actin filaments measured from thermal fluctuations in shape. The Journal of cell biology. 1993;120(4):923–934. pmid:8432732
- View Article
- PubMed/NCBI
- Google Scholar
58. Le S, Hu X, Yao M, Chen H, Yu M, Xu X, et al. Mechanotransmission and mechanosensing of human alpha-actinin 1. Cell reports. 2017;21(10):2714–2723. pmid:29212020
- View Article
- PubMed/NCBI
- Google Scholar
59. Meyer RK, Aebi U. Bundling of actin filaments by alpha-actinin depends on its molecular length. The Journal of cell biology. 1990;110(6):2013–2024. pmid:2351691
- View Article
- PubMed/NCBI
- Google Scholar
60. Kuhlman PA, Ellis J, Critchley DR, Bagshaw CR. The kinetics of the interaction between the actin-binding domain of α-actinin and F-actin. FEBS letters. 1994;339(3):297–301. pmid:8112470
- View Article
- PubMed/NCBI
- Google Scholar
61. Xu J, Wirtz D, Pollard TD. Dynamic cross-linking by α-actinin determines the mechanical properties of actin filament networks. Journal of Biological Chemistry. 1998;273(16):9570–9576. pmid:9545287
- View Article
- PubMed/NCBI
- Google Scholar
62. Gillespie DT. Stochastic simulation of chemical kinetics. Annu Rev Phys Chem. 2007;58:35–55. pmid:17037977
- View Article
- PubMed/NCBI
- Google Scholar
63. Anderson DF. A modified next reaction method for simulating chemical systems with time dependent propensities and delays. The Journal of chemical physics. 2007;127(21):214107. pmid:18067349
- View Article
- PubMed/NCBI
- Google Scholar
64. Donev A, Yang CY, Kim C. Efficient reactive Brownian dynamics. The Journal of chemical physics. 2018;148(3):034103. pmid:29352800
- View Article
- PubMed/NCBI
- Google Scholar
65. Allen MP, Tildesley DJ. Computer simulation of liquids. Oxford university press; 2017.
66. Kim T, Hwang W, Lee H, Kamm RD. Computational analysis of viscoelastic properties of crosslinked actin networks. PLoS Comput Biol. 2009;5(7):e1000439. pmid:19609348
- View Article
- PubMed/NCBI
- Google Scholar
67. Courson DS, Rock RS. Actin cross-link assembly and disassembly mechanics for α-actinin and fascin. Journal of Biological Chemistry. 2010;285(34):26350–26357. pmid:20551315
- View Article
- PubMed/NCBI
- Google Scholar
68. Johnston AB, Collins A, Goode BL. High-speed depolymerization at actin filament ends jointly catalysed by Twinfilin and Srv2/CAP. Nature cell biology. 2015;17(11):1504–1511. pmid:26458246
- View Article
- PubMed/NCBI
- Google Scholar
69. Kueh HY, Charras GT, Mitchison TJ, Brieher WM. Actin disassembly by cofilin, coronin, and Aip1 occurs in bursts and is inhibited by barbed-end cappers. The Journal of cell biology. 2008;182(2):341–353. pmid:18663144
- View Article
- PubMed/NCBI
- Google Scholar
70. Manhart A, Icheva A, Guerin C, Klar T, Boujemaa-Paterski R, Thery M, et al. Reconstitution of the equilibrium state of dynamic actin networks. bioRxiv. 2018; p. 437806.
- View Article
- Google Scholar
71. Nazockdast E, Rahimian A, Zorin D, Shelley M. A fast platform for simulating semi-flexible fiber suspensions applied to cell mechanics. J Comput Phys. 2017;329:173–209.
- View Article
- Google Scholar
72. Keller JB, Rubinow SI. Slender-body theory for slow viscous flow. Journal of Fluid Mechanics. 1976;75(4):705–714.
- View Article
- Google Scholar
73. Johnson RE. An improved slender-body theory for Stokes flow. Journal of Fluid Mechanics. 1980;99(2):411–431.
- View Article
- Google Scholar
74. Wajnryb E, Mizerski KA, Zuk PJ, Szymczak P. Generalization of the Rotne–Prager–Yamakawa mobility and shear disturbance tensors. Journal of Fluid Mechanics. 2013;731.
- View Article
- Google Scholar
75. Fiore AM, Balboa Usabiaga F, Donev A, Swan JW. Rapid sampling of stochastic displacements in Brownian dynamics simulations. The Journal of chemical physics. 2017;146(12):124116. pmid:28388117
- View Article
- PubMed/NCBI
- Google Scholar
76. Fiore AM, Swan JW. Rapid sampling of stochastic displacements in Brownian dynamics simulations with stresslet constraints. The Journal of chemical physics. 2018;148(4):044114. pmid:29390810
- View Article
- PubMed/NCBI
- Google Scholar
77. af Klinteberg L, Barnett AH. Accurate quadrature of nearly singular line integrals in two and three dimensions by singularity swapping. BIT Numerical Mathematics. 2020; p. 1–36.
- View Article
- Google Scholar
78. Perez RP. Universally Adaptable Multiscale Molecular Dynamics (UAMMD); 2021. https://github.com/RaulPPelaez/UAMMD.
79. Stam S, Alberts J, Gardel ML, Munro E. Isoforms confer characteristic force generation and mechanosensation by myosin II filaments. Biophysical journal. 2015;108(8):1997–2006. pmid:25902439
- View Article
- PubMed/NCBI
- Google Scholar
80. Morrison FA, et al. Understanding rheology. Oxford University Press, USA; 2001.
81. Batchelor G. The stress system in a suspension of force-free particles. Journal of fluid mechanics. 1970;41(3):545–570.
- View Article
- Google Scholar
82. Holden H, Karlsen KH, Lie KA. Splitting methods for partial differential equations with rough solutions: Analysis and MATLAB programs. vol. 11. European Mathematical Society; 2010.
83. Lieleg O, Schmoller KM, Cyron CJ, Luan Y, Wall WA, Bausch AR. Structural polymorphism in heterogeneous cytoskeletal networks. Soft Matter. 2009;5(9):1796–1803.
- View Article
- Google Scholar
84. Schmoller K, Lieleg O, Bausch A. Structural and viscoelastic properties of actin/filamin networks: cross-linked versus bundled networks. Biophysical journal. 2009;97(1):83–89. pmid:19580746
- View Article
- PubMed/NCBI
- Google Scholar
85. Ma R, Berro J. Structural organization and energy storage in crosslinked actin assemblies. PLoS computational biology. 2018;14(5):e1006150. pmid:29813051
- View Article
- PubMed/NCBI
- Google Scholar
86. Batchelor G. Slender-body theory for particles of arbitrary cross-section in Stokes flow. Journal of Fluid Mechanics. 1970;44(3):419–440.
- View Article
- Google Scholar
87. Shaqfeh ES, Fredrickson GH. The hydrodynamic stress in a suspension of rods. Physics of Fluids A: Fluid Dynamics. 1990;2(1):7–24.
- View Article
- Google Scholar
88. Mackaplow MB, Shaqfeh ES. A numerical study of the rheological properties of suspensions of rigid, non-Brownian fibres. Journal of Fluid Mechanics. 1996;329:155–186.
- View Article
- Google Scholar
89. Winter HH. Analysis of dynamic mechanical data: inversion into a relaxation time spectrum and consistency check. Journal of Non-Newtonian Fluid Mechanics. 1997;68(2-3):225–239.
- View Article
- Google Scholar
90. McDougall I, Orbey N, Dealy JM. Inferring meaningful relaxation spectra from experimental data. Journal of Rheology. 2014;58(3):779–797.
- View Article
- Google Scholar
91. Baumgaertel M, Winter H. Determination of discrete relaxation and retardation time spectra from dynamic mechanical data. Rheologica Acta. 1989;28(6):511–519.
- View Article
- Google Scholar
92. Malkin AY. The use of a continuous relaxation spectrum for describing the viscoelastic properties of polymers. Polymer Science Series A. 2006;48(1):39–45.
- View Article
- Google Scholar
93. Henty JL, Bledsoe SW, Khurana P, Meagher RB, Day B, Blanchoin L, et al. Arabidopsis actin depolymerizing factor4 modulates the stochastic dynamic behavior of actin filaments in the cortical array of epidermal cells. The Plant Cell. 2011;23(10):3711–3726. pmid:22010035
- View Article
- PubMed/NCBI
- Google Scholar
94. Kim T, Gardel ML, Munro E. Determinants of fluidlike behavior and effective viscosity in cross-linked actin networks. Biophysical journal. 2014;106(3):526–534. pmid:24507593
- View Article
- PubMed/NCBI
- Google Scholar
95. Stricker J, Falzone T, Gardel ML. Mechanics of the F-actin cytoskeleton. Journal of biomechanics. 2010;43(1):9–14. pmid:19913792
- View Article
- PubMed/NCBI
- Google Scholar
96. Åström JA, Kumar PS, Vattulainen I, Karttunen M. Strain hardening, avalanches, and strain softening in dense cross-linked actin networks. Physical Review E. 2008;77(5):051913. pmid:18643108
- View Article
- PubMed/NCBI
- Google Scholar
97. Scheff DR, Redford SA, Lorpaiboon C, Majumdar S, Dinner AR, Gardel ML. Actin filament alignment causes mechanical hysteresis in cross-linked networks. Soft Matter. 2021;17(22):5499–5507. pmid:33989373
- View Article
- PubMed/NCBI
- Google Scholar
98. Tee YH, Shemesh T, Thiagarajan V, Hariadi RF, Anderson KL, Page C, et al. Cellular chirality arising from the self-organization of the actin cytoskeleton. Nature cell biology. 2015;17(4):445–457. pmid:25799062
- View Article
- PubMed/NCBI
- Google Scholar
99. Weirich KL, Banerjee S, Dasbiswas K, Witten TA, Vaikuntanathan S, Gardel ML. Liquid behavior of cross-linked actin bundles. Proceedings of the National Academy of Sciences. 2017;114(9):2131–2136. pmid:28202730
- View Article
- PubMed/NCBI
- Google Scholar
100. Weirich KL, Stam S, Munro E, Gardel ML. Actin bundle architecture and mechanics regulate myosin II force generation. Biophysical Journal. 2021;120(10):1957–1970. pmid:33798565
- View Article
- PubMed/NCBI
- Google Scholar

[ref1] 1. Mofrad MR. Rheology of the cytoskeleton. Annual Review of Fluid Mechanics. 2009;41:433–453.
View Article
Google Scholar

[2] View Article

[3] Google Scholar

[ref2] 2. Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P. Molecular biology of the cell. Garland Science; 2002.

[ref3] 3. Hinner B, Tempel M, Sackmann E, Kroy K, Frey E. Entanglement, elasticity, and viscous relaxation of actin solutions. Physical Review Letters. 1998;81(12):2614.
View Article
Google Scholar

[6] View Article

[7] Google Scholar

[ref4] 4. Janmey PA, Hvidt S, Käs J, Lerche D, Maggs A, Sackmann E, et al. The mechanical properties of actin gels. Elastic modulus and filament motions. Journal of Biological Chemistry. 1994;269(51):32503–32513. pmid:7798252
View Article
PubMed/NCBI
Google Scholar

[9] View Article

[10] PubMed/NCBI

[11] Google Scholar

[ref5] 5. Wachsstock DH, Schwartz W, Pollard TD. Affinity of alpha-actinin for actin determines the structure and mechanical properties of actin filament gels. Biophysical journal. 1993;65(1):205.
View Article
Google Scholar

[13] View Article

[14] Google Scholar

[ref6] 6. Lieleg O, Claessens MMAE, Luan Y, Bausch A. Transient binding and dissipation in cross-linked actin networks. Physical review letters. 2008;101(10):108101. pmid:18851260
View Article
PubMed/NCBI
Google Scholar

[16] View Article

[17] PubMed/NCBI

[18] Google Scholar

[ref7] 7. Chaubet L, Chaudhary AR, Heris HK, Ehrlicher AJ, Hendricks AG. Dynamic actin cross-linking governs the cytoplasm’s transition to fluid-like behavior. Molecular biology of the cell. 2020;31(16):1744–1752. pmid:32579489
View Article
PubMed/NCBI
Google Scholar

[20] View Article

[21] PubMed/NCBI

[22] Google Scholar

[ref8] 8. Wachsstock DH, Schwarz W, Pollard T. Cross-linker dynamics determine the mechanical properties of actin gels. Biophysical journal. 1994;66(3):801–809. pmid:8011912
View Article
PubMed/NCBI
Google Scholar

[24] View Article

[25] PubMed/NCBI

[26] Google Scholar

[ref9] 9. Gardel M, Shin JH, MacKintosh F, Mahadevan L, Matsudaira P, Weitz D. Elastic behavior of cross-linked and bundled actin networks. Science. 2004;304(5675):1301–1305. pmid:15166374
View Article
PubMed/NCBI
Google Scholar

[28] View Article

[29] PubMed/NCBI

[30] Google Scholar

[ref10] 10. Gardel M, Shin JH, MacKintosh F, Mahadevan L, Matsudaira P, Weitz D. Scaling of F-actin network rheology to probe single filament elasticity and dynamics. Physical review letters. 2004;93(18):188102. pmid:15525211
View Article
PubMed/NCBI
Google Scholar

[32] View Article

[33] PubMed/NCBI

[34] Google Scholar

[ref11] 11. Müller O, Gaub H, Bärmann M, Sackmann E. Viscoelastic moduli of sterically and chemically cross-linked actin networks in the dilute to semidilute regime: measurements by oscillating disk rheometer. Macromolecules. 1991;24(11):3111–3120.
View Article
Google Scholar

[36] View Article

[37] Google Scholar

[ref12] 12. Ruddies R, Goldmann W, Isenberg G, Sackmann E. The viscoelasticity of entangled actin networks: the influence of defects and modulation by talin and vinculin. European biophysics journal. 1993;22(5):309–321. pmid:8112218
View Article
PubMed/NCBI
Google Scholar

[39] View Article

[40] PubMed/NCBI

[41] Google Scholar

[ref13] 13. Yao NY, Becker DJ, Broedersz CP, Depken M, MacKintosh FC, Pollak MR, et al. Nonlinear viscoelasticity of actin transiently cross-linked with mutant α-actinin-4. Journal of molecular biology. 2011;411(5):1062–1071. pmid:21762701
View Article
PubMed/NCBI
Google Scholar

[43] View Article

[44] PubMed/NCBI

[45] Google Scholar

[ref14] 14. Mulla Y, MacKintosh F, Koenderink GH. Origin of slow stress relaxation in the cytoskeleton. Physical review letters. 2019;122(21):218102. pmid:31283330
View Article
PubMed/NCBI
Google Scholar

[47] View Article

[48] PubMed/NCBI

[49] Google Scholar

[ref15] 15. Lieleg O, Schmoller K, Claessens MMAE, Bausch A. Cytoskeletal polymer networks: viscoelastic properties are determined by the microscopic interaction potential of cross-links. Biophysical journal. 2009;96(11):4725–4732. pmid:19486695
View Article
PubMed/NCBI
Google Scholar

[51] View Article

[52] PubMed/NCBI

[53] Google Scholar

[ref16] 16. Yao NY, Broedersz CP, Depken M, Becker DJ, Pollak MR, MacKintosh FC, et al. Stress-enhanced gelation: A dynamic nonlinearity of elasticity. Physical review letters. 2013;110(1):018103. pmid:23383843
View Article
PubMed/NCBI
Google Scholar

[55] View Article

[56] PubMed/NCBI

[57] Google Scholar

[ref17] 17. Morse DC. Viscoelasticity of concentrated isotropic solutions of semiflexible polymers. 1. Model and stress tensor. Macromolecules. 1998;31(20):7030–7043.
View Article
Google Scholar

[59] View Article

[60] Google Scholar

[ref18] 18. Fabry B, Maksym GN, Butler JP, Glogauer M, Navajas D, Fredberg JJ. Scaling the microrheology of living cells. Physical review letters. 2001;87(14):148102. pmid:11580676
View Article
PubMed/NCBI
Google Scholar

[62] View Article

[63] PubMed/NCBI

[64] Google Scholar

[ref19] 19. Fabry B, Maksym GN, Butler JP, Glogauer M, Navajas D, Taback NA, et al. Time scale and other invariants of integrative mechanical behavior in living cells. Physical Review E. 2003;68(4):041914. pmid:14682980
View Article
PubMed/NCBI
Google Scholar

[66] View Article

[67] PubMed/NCBI

[68] Google Scholar

[ref20] 20. Ziemann F, Rädler J, Sackmann E. Local measurements of viscoelastic moduli of entangled actin networks using an oscillating magnetic bead micro-rheometer. Biophysical journal. 1994;66(6):2210–2216. pmid:8075354
View Article
PubMed/NCBI
Google Scholar

[70] View Article

[71] PubMed/NCBI

[72] Google Scholar

[ref21] 21. Rouse PE Jr. A theory of the linear viscoelastic properties of dilute solutions of coiling polymers. The Journal of Chemical Physics. 1953;21(7):1272–1280.
View Article
Google Scholar

[74] View Article

[75] Google Scholar

[ref22] 22. Freedman SL, Suarez C, Winkelman JD, Kovar DR, Voth GA, Dinner AR, et al. Mechanical and kinetic factors drive sorting of F-actin cross-linkers on bundles. Proceedings of the National Academy of Sciences. 2019;116(33):16192–16197. pmid:31346091
View Article
PubMed/NCBI
Google Scholar

[77] View Article

[78] PubMed/NCBI

[79] Google Scholar

[ref23] 23. Kasza K, Broedersz C, Koenderink G, Lin Y, Messner W, Millman E, et al. Actin filament length tunes elasticity of flexibly cross-linked actin networks. Biophysical journal. 2010;99(4):1091–1100. pmid:20712992
View Article
PubMed/NCBI
Google Scholar

[81] View Article

[82] PubMed/NCBI

[83] Google Scholar

[ref24] 24. Wagner B, Tharmann R, Haase I, Fischer M, Bausch AR. Cytoskeletal polymer networks: the molecular structure of cross-linkers determines macroscopic properties. Proceedings of the National Academy of Sciences. 2006;103(38):13974–13978. pmid:16963567
View Article
PubMed/NCBI
Google Scholar

[85] View Article

[86] PubMed/NCBI

[87] Google Scholar

[ref25] 25. Tseng Y, An KM, Wirtz D. Microheterogeneity controls the rate of gelation of actin filament networks. Journal of Biological Chemistry. 2002;277(20):18143–18150. pmid:11889122
View Article
PubMed/NCBI
Google Scholar

[89] View Article

[90] PubMed/NCBI

[91] Google Scholar

[ref26] 26. Tseng Y, Schafer BW, Almo SC, Wirtz D. Functional synergy of actin filament cross-linking proteins. Journal of Biological Chemistry. 2002;277(28):25609–25616. pmid:12006593
View Article
PubMed/NCBI
Google Scholar

[93] View Article

[94] PubMed/NCBI

[95] Google Scholar

[ref27] 27. MacKintosh F, Käs J, Janmey P. Elasticity of semiflexible biopolymer networks. Physical review letters. 1995;75(24):4425. pmid:10059905
View Article
PubMed/NCBI
Google Scholar

[97] View Article

[98] PubMed/NCBI

[99] Google Scholar

[ref28] 28. Isambert H, Maggs A. Dynamics and rheology of actin solutions. Macromolecules. 1996;29(3):1036–1040.
View Article
Google Scholar

[101] View Article

[102] Google Scholar

[ref29] 29. Pegoraro AF, Janmey P, Weitz DA. Mechanical properties of the cytoskeleton and cells. Cold Spring Harbor perspectives in biology. 2017;9(11):a022038. pmid:29092896
View Article
PubMed/NCBI
Google Scholar

[104] View Article

[105] PubMed/NCBI

[106] Google Scholar

[ref30] 30. Broedersz CP, Depken M, Yao NY, Pollak MR, Weitz DA, MacKintosh FC. Cross-link-governed dynamics of biopolymer networks. Physical review letters. 2010;105(23):238101. pmid:21231506
View Article
PubMed/NCBI
Google Scholar

[108] View Article

[109] PubMed/NCBI

[110] Google Scholar

[ref31] 31. Müller KW, Bruinsma RF, Lieleg O, Bausch AR, Wall WA, Levine AJ. Rheology of semiflexible bundle networks with transient linkers. Physical review letters. 2014;112(23):238102. pmid:24972229
View Article
PubMed/NCBI
Google Scholar

[112] View Article

[113] PubMed/NCBI

[114] Google Scholar

[ref32] 32. Peskin CS. The immersed boundary method. Acta Numer. 2002;11:479–517.
View Article
Google Scholar

[116] View Article

[117] Google Scholar

[ref33] 33. Karcher H, Lammerding J, Huang H, Lee RT, Kamm RD, Kaazempur-Mofrad MR. A three-dimensional viscoelastic model for cell deformation with experimental verification. Biophysical journal. 2003;85(5):3336–3349. pmid:14581235
View Article
PubMed/NCBI
Google Scholar

[119] View Article

[120] PubMed/NCBI

[121] Google Scholar

[ref34] 34. Mijailovich SM, Kojic M, Zivkovic M, Fabry B, Fredberg JJ. A finite element model of cell deformation during magnetic bead twisting. Journal of Applied Physiology. 2002;93(4):1429–1436. pmid:12235044
View Article
PubMed/NCBI
Google Scholar

[123] View Article

[124] PubMed/NCBI

[125] Google Scholar

[ref35] 35. Strychalski W, Copos CA, Lewis OL, Guy RD. A poroelastic immersed boundary method with applications to cell biology. Journal of Computational Physics. 2015;282:77–97.
View Article
Google Scholar

[127] View Article

[128] Google Scholar

[ref36] 36. Copos CA, Guy RD. A porous viscoelastic model for the cell cytoskeleton. The ANZIAM Journal. 2018;59(4):472–498.
View Article
Google Scholar

[130] View Article

[131] Google Scholar

[ref37] 37. Guy RD, Nakagaki T, Wright GB. Flow-induced channel formation in the cytoplasm of motile cells. Physical Review E. 2011;84(1):016310. pmid:21867307
View Article
PubMed/NCBI
Google Scholar

[133] View Article

[134] PubMed/NCBI

[135] Google Scholar

[ref38] 38. Desprat N, Richert A, Simeon J, Asnacios A. Creep function of a single living cell. Biophysical journal. 2005;88(3):2224–2233. pmid:15596508
View Article
PubMed/NCBI
Google Scholar

[137] View Article

[138] PubMed/NCBI

[139] Google Scholar

[ref39] 39. Head DA, Levine AJ, MacKintosh F. Deformation of cross-linked semiflexible polymer networks. Physical review letters. 2003;91(10):108102. pmid:14525510
View Article
PubMed/NCBI
Google Scholar

[141] View Article

[142] PubMed/NCBI

[143] Google Scholar

[ref40] 40. Kim T, Hwang W, Kamm R. Computational analysis of a cross-linked actin-like network. Experimental Mechanics. 2009;49(1):91–104.
View Article
Google Scholar

[145] View Article

[146] Google Scholar

[ref41] 41. Popov K, Komianos J, Papoian GA. MEDYAN: Mechanochemical simulations of contraction and polarity alignment in actomyosin networks. PLoS computational biology. 2016;12(4):e1004877. pmid:27120189
View Article
PubMed/NCBI
Google Scholar

[148] View Article

[149] PubMed/NCBI

[150] Google Scholar

[ref42] 42. Li X, Ni Q, He X, Kong J, Lim SM, Papoian GA, et al. Tensile force-induced cytoskeletal remodeling: Mechanics before chemistry. PLoS computational biology. 2020;16(6):e1007693. pmid:32520928
View Article
PubMed/NCBI
Google Scholar

[152] View Article

[153] PubMed/NCBI

[154] Google Scholar

[ref43] 43. Wei X, Fang C, Gong B, Yao J, Qian J, Lin Y. Viscoelasticity of 3D actin networks dictated by the mechanochemical characteristics of cross-linkers. Soft Matter. 2021;. pmid:33646227
View Article
PubMed/NCBI
Google Scholar

[156] View Article

[157] PubMed/NCBI

[158] Google Scholar

[ref44] 44. Maxian O, Mogilner A, Donev A. Integral-based spectral method for inextensible slender fibers in Stokes flow. Physical Review Fluids. 2021;6(1):014102.
View Article
Google Scholar

[160] View Article

[161] Google Scholar

[ref45] 45. Lieleg O, Claessens MM, Bausch AR. Structure and dynamics of cross-linked actin networks. Soft Matter. 2010;6(2):218–225.
View Article
Google Scholar

[163] View Article

[164] Google Scholar

[ref46] 46. Guha M, Zhou M, Wang Yl. Cortical actin turnover during cytokinesis requires myosin II. Current biology. 2005;15(8):732–736. pmid:15854905
View Article
PubMed/NCBI
Google Scholar

[166] View Article

[167] PubMed/NCBI

[168] Google Scholar

[ref47] 47. Theriot JA, Mitchison TJ. Actin microfilament dynamics in locomoting cells. Nature. 1991;352(6331):126–131. pmid:2067574
View Article
PubMed/NCBI
Google Scholar

[170] View Article

[171] PubMed/NCBI

[172] Google Scholar

[ref48] 48. Tornberg AK, Shelley MJ. Simulating the dynamics and interactions of flexible fibers in Stokes flows. Journal of Computational Physics. 2004;196(1):8–40.
View Article
Google Scholar

[174] View Article

[175] Google Scholar

[ref49] 49. Aurentz JL, Trefethen LN. Block operators and spectral discretizations. SIAM Review. 2017;59(2):423–446.
View Article
Google Scholar

[177] View Article

[178] Google Scholar

[ref50] 50. Driscoll TA, Hale N. Rectangular spectral collocation. IMA Journal of Numerical Analysis. 2015;36(1):108–132.
View Article
Google Scholar

[180] View Article

[181] Google Scholar

[ref51] 51. Zero K, Pecora R. Rotational and translational diffusion in semidilute solutions of rigid-rod macromolecules. Macromolecules. 1982;15(1):87–93.
View Article
Google Scholar

[183] View Article

[184] Google Scholar

[ref52] 52. Grazi E. What is the diameter of the actin filament? FEBS letters. 1997;405(3):249–252. pmid:9108298
View Article
PubMed/NCBI
Google Scholar

[186] View Article

[187] PubMed/NCBI

[188] Google Scholar

[ref53] 53. McGrath JL, Osborn EA, Tardy YS, Dewey CF, Hartwig JH. Regulation of the actin cycle in vivo by actin filament severing. Proceedings of the National Academy of Sciences. 2000;97(12):6532–6537. pmid:10823888
View Article
PubMed/NCBI
Google Scholar

[190] View Article

[191] PubMed/NCBI

[192] Google Scholar

[ref54] 54. Yoshida A, Sakai N, Uekusa Y, Deguchi K, Gilmore JL, Kumeta M, et al. Probing in vivo dynamics of mitochondria and cortical actin networks using high-speed atomic force/fluorescence microscopy. Genes to Cells. 2015;20(2):85–94. pmid:25440894
View Article
PubMed/NCBI
Google Scholar

[194] View Article

[195] PubMed/NCBI

[196] Google Scholar

[ref55] 55. Zhang Y, Yoshida A, Sakai N, Uekusa Y, Kumeta M, Yoshimura SH. In vivo dynamics of the cortical actin network revealed by fast-scanning atomic force microscopy. Microscopy. 2017;66(4):272–282. pmid:28531263
View Article
PubMed/NCBI
Google Scholar

[198] View Article

[199] PubMed/NCBI

[200] Google Scholar

[ref56] 56. Luby-Phelps K. Cytoarchitecture and physical properties of cytoplasm: volume, viscosity, diffusion, intracellular surface area. In: International review of cytology. vol. 192. Elsevier; 1999. p. 189–221.

[ref57] 57. Gittes F, Mickey B, Nettleton J, Howard J. Flexural rigidity of microtubules and actin filaments measured from thermal fluctuations in shape. The Journal of cell biology. 1993;120(4):923–934. pmid:8432732
View Article
PubMed/NCBI
Google Scholar

[203] View Article

[204] PubMed/NCBI

[205] Google Scholar

[ref58] 58. Le S, Hu X, Yao M, Chen H, Yu M, Xu X, et al. Mechanotransmission and mechanosensing of human alpha-actinin 1. Cell reports. 2017;21(10):2714–2723. pmid:29212020
View Article
PubMed/NCBI
Google Scholar

[207] View Article

[208] PubMed/NCBI

[209] Google Scholar

[ref59] 59. Meyer RK, Aebi U. Bundling of actin filaments by alpha-actinin depends on its molecular length. The Journal of cell biology. 1990;110(6):2013–2024. pmid:2351691
View Article
PubMed/NCBI
Google Scholar

[211] View Article

[212] PubMed/NCBI

[213] Google Scholar

[ref60] 60. Kuhlman PA, Ellis J, Critchley DR, Bagshaw CR. The kinetics of the interaction between the actin-binding domain of α-actinin and F-actin. FEBS letters. 1994;339(3):297–301. pmid:8112470
View Article
PubMed/NCBI
Google Scholar

[215] View Article

[216] PubMed/NCBI

[217] Google Scholar

[ref61] 61. Xu J, Wirtz D, Pollard TD. Dynamic cross-linking by α-actinin determines the mechanical properties of actin filament networks. Journal of Biological Chemistry. 1998;273(16):9570–9576. pmid:9545287
View Article
PubMed/NCBI
Google Scholar

[219] View Article

[220] PubMed/NCBI

[221] Google Scholar

[ref62] 62. Gillespie DT. Stochastic simulation of chemical kinetics. Annu Rev Phys Chem. 2007;58:35–55. pmid:17037977
View Article
PubMed/NCBI
Google Scholar

[223] View Article

[224] PubMed/NCBI

[225] Google Scholar

[ref63] 63. Anderson DF. A modified next reaction method for simulating chemical systems with time dependent propensities and delays. The Journal of chemical physics. 2007;127(21):214107. pmid:18067349
View Article
PubMed/NCBI
Google Scholar

[227] View Article

[228] PubMed/NCBI

[229] Google Scholar

[ref64] 64. Donev A, Yang CY, Kim C. Efficient reactive Brownian dynamics. The Journal of chemical physics. 2018;148(3):034103. pmid:29352800
View Article
PubMed/NCBI
Google Scholar

[231] View Article

[232] PubMed/NCBI

[233] Google Scholar

[ref65] 65. Allen MP, Tildesley DJ. Computer simulation of liquids. Oxford university press; 2017.

[ref66] 66. Kim T, Hwang W, Lee H, Kamm RD. Computational analysis of viscoelastic properties of crosslinked actin networks. PLoS Comput Biol. 2009;5(7):e1000439. pmid:19609348
View Article
PubMed/NCBI
Google Scholar

[236] View Article

[237] PubMed/NCBI

[238] Google Scholar

[ref67] 67. Courson DS, Rock RS. Actin cross-link assembly and disassembly mechanics for α-actinin and fascin. Journal of Biological Chemistry. 2010;285(34):26350–26357. pmid:20551315
View Article
PubMed/NCBI
Google Scholar

[240] View Article

[241] PubMed/NCBI

[242] Google Scholar

[ref68] 68. Johnston AB, Collins A, Goode BL. High-speed depolymerization at actin filament ends jointly catalysed by Twinfilin and Srv2/CAP. Nature cell biology. 2015;17(11):1504–1511. pmid:26458246
View Article
PubMed/NCBI
Google Scholar

[244] View Article

[245] PubMed/NCBI

[246] Google Scholar

[ref69] 69. Kueh HY, Charras GT, Mitchison TJ, Brieher WM. Actin disassembly by cofilin, coronin, and Aip1 occurs in bursts and is inhibited by barbed-end cappers. The Journal of cell biology. 2008;182(2):341–353. pmid:18663144
View Article
PubMed/NCBI
Google Scholar

[248] View Article

[249] PubMed/NCBI

[250] Google Scholar

[ref70] 70. Manhart A, Icheva A, Guerin C, Klar T, Boujemaa-Paterski R, Thery M, et al. Reconstitution of the equilibrium state of dynamic actin networks. bioRxiv. 2018; p. 437806.
View Article
Google Scholar

[252] View Article

[253] Google Scholar

[ref71] 71. Nazockdast E, Rahimian A, Zorin D, Shelley M. A fast platform for simulating semi-flexible fiber suspensions applied to cell mechanics. J Comput Phys. 2017;329:173–209.
View Article
Google Scholar

[255] View Article

[256] Google Scholar

[ref72] 72. Keller JB, Rubinow SI. Slender-body theory for slow viscous flow. Journal of Fluid Mechanics. 1976;75(4):705–714.
View Article
Google Scholar

[258] View Article

[259] Google Scholar

[ref73] 73. Johnson RE. An improved slender-body theory for Stokes flow. Journal of Fluid Mechanics. 1980;99(2):411–431.
View Article
Google Scholar

[261] View Article

[262] Google Scholar

[ref74] 74. Wajnryb E, Mizerski KA, Zuk PJ, Szymczak P. Generalization of the Rotne–Prager–Yamakawa mobility and shear disturbance tensors. Journal of Fluid Mechanics. 2013;731.
View Article
Google Scholar

[264] View Article

[265] Google Scholar

[ref75] 75. Fiore AM, Balboa Usabiaga F, Donev A, Swan JW. Rapid sampling of stochastic displacements in Brownian dynamics simulations. The Journal of chemical physics. 2017;146(12):124116. pmid:28388117
View Article
PubMed/NCBI
Google Scholar

[267] View Article

[268] PubMed/NCBI

[269] Google Scholar

[ref76] 76. Fiore AM, Swan JW. Rapid sampling of stochastic displacements in Brownian dynamics simulations with stresslet constraints. The Journal of chemical physics. 2018;148(4):044114. pmid:29390810
View Article
PubMed/NCBI
Google Scholar

[271] View Article

[272] PubMed/NCBI

[273] Google Scholar

[ref77] 77. af Klinteberg L, Barnett AH. Accurate quadrature of nearly singular line integrals in two and three dimensions by singularity swapping. BIT Numerical Mathematics. 2020; p. 1–36.
View Article
Google Scholar

[275] View Article

[276] Google Scholar

[ref78] 78. Perez RP. Universally Adaptable Multiscale Molecular Dynamics (UAMMD); 2021. https://github.com/RaulPPelaez/UAMMD.

[ref79] 79. Stam S, Alberts J, Gardel ML, Munro E. Isoforms confer characteristic force generation and mechanosensation by myosin II filaments. Biophysical journal. 2015;108(8):1997–2006. pmid:25902439
View Article
PubMed/NCBI
Google Scholar

[279] View Article

[280] PubMed/NCBI

[281] Google Scholar

[ref80] 80. Morrison FA, et al. Understanding rheology. Oxford University Press, USA; 2001.

[ref81] 81. Batchelor G. The stress system in a suspension of force-free particles. Journal of fluid mechanics. 1970;41(3):545–570.
View Article
Google Scholar

[284] View Article

[285] Google Scholar

[ref82] 82. Holden H, Karlsen KH, Lie KA. Splitting methods for partial differential equations with rough solutions: Analysis and MATLAB programs. vol. 11. European Mathematical Society; 2010.

[ref83] 83. Lieleg O, Schmoller KM, Cyron CJ, Luan Y, Wall WA, Bausch AR. Structural polymorphism in heterogeneous cytoskeletal networks. Soft Matter. 2009;5(9):1796–1803.
View Article
Google Scholar

[288] View Article

[289] Google Scholar

[ref84] 84. Schmoller K, Lieleg O, Bausch A. Structural and viscoelastic properties of actin/filamin networks: cross-linked versus bundled networks. Biophysical journal. 2009;97(1):83–89. pmid:19580746
View Article
PubMed/NCBI
Google Scholar

[291] View Article

[292] PubMed/NCBI

[293] Google Scholar

[ref85] 85. Ma R, Berro J. Structural organization and energy storage in crosslinked actin assemblies. PLoS computational biology. 2018;14(5):e1006150. pmid:29813051
View Article
PubMed/NCBI
Google Scholar

[295] View Article

[296] PubMed/NCBI

[297] Google Scholar

[ref86] 86. Batchelor G. Slender-body theory for particles of arbitrary cross-section in Stokes flow. Journal of Fluid Mechanics. 1970;44(3):419–440.
View Article
Google Scholar

[299] View Article

[300] Google Scholar

[ref87] 87. Shaqfeh ES, Fredrickson GH. The hydrodynamic stress in a suspension of rods. Physics of Fluids A: Fluid Dynamics. 1990;2(1):7–24.
View Article
Google Scholar

[302] View Article

[303] Google Scholar

[ref88] 88. Mackaplow MB, Shaqfeh ES. A numerical study of the rheological properties of suspensions of rigid, non-Brownian fibres. Journal of Fluid Mechanics. 1996;329:155–186.
View Article
Google Scholar

[305] View Article

[306] Google Scholar

[ref89] 89. Winter HH. Analysis of dynamic mechanical data: inversion into a relaxation time spectrum and consistency check. Journal of Non-Newtonian Fluid Mechanics. 1997;68(2-3):225–239.
View Article
Google Scholar

[308] View Article

[309] Google Scholar

[ref90] 90. McDougall I, Orbey N, Dealy JM. Inferring meaningful relaxation spectra from experimental data. Journal of Rheology. 2014;58(3):779–797.
View Article
Google Scholar

[311] View Article

[312] Google Scholar

[ref91] 91. Baumgaertel M, Winter H. Determination of discrete relaxation and retardation time spectra from dynamic mechanical data. Rheologica Acta. 1989;28(6):511–519.
View Article
Google Scholar

[314] View Article

[315] Google Scholar

[ref92] 92. Malkin AY. The use of a continuous relaxation spectrum for describing the viscoelastic properties of polymers. Polymer Science Series A. 2006;48(1):39–45.
View Article
Google Scholar

[317] View Article

[318] Google Scholar

[ref93] 93. Henty JL, Bledsoe SW, Khurana P, Meagher RB, Day B, Blanchoin L, et al. Arabidopsis actin depolymerizing factor4 modulates the stochastic dynamic behavior of actin filaments in the cortical array of epidermal cells. The Plant Cell. 2011;23(10):3711–3726. pmid:22010035
View Article
PubMed/NCBI
Google Scholar

[320] View Article

[321] PubMed/NCBI

[322] Google Scholar

[ref94] 94. Kim T, Gardel ML, Munro E. Determinants of fluidlike behavior and effective viscosity in cross-linked actin networks. Biophysical journal. 2014;106(3):526–534. pmid:24507593
View Article
PubMed/NCBI
Google Scholar

[324] View Article

[325] PubMed/NCBI

[326] Google Scholar

[ref95] 95. Stricker J, Falzone T, Gardel ML. Mechanics of the F-actin cytoskeleton. Journal of biomechanics. 2010;43(1):9–14. pmid:19913792
View Article
PubMed/NCBI
Google Scholar

[328] View Article

[329] PubMed/NCBI

[330] Google Scholar

[ref96] 96. Åström JA, Kumar PS, Vattulainen I, Karttunen M. Strain hardening, avalanches, and strain softening in dense cross-linked actin networks. Physical Review E. 2008;77(5):051913. pmid:18643108
View Article
PubMed/NCBI
Google Scholar

[332] View Article

[333] PubMed/NCBI

[334] Google Scholar

[ref97] 97. Scheff DR, Redford SA, Lorpaiboon C, Majumdar S, Dinner AR, Gardel ML. Actin filament alignment causes mechanical hysteresis in cross-linked networks. Soft Matter. 2021;17(22):5499–5507. pmid:33989373
View Article
PubMed/NCBI
Google Scholar

[336] View Article

[337] PubMed/NCBI

[338] Google Scholar

[ref98] 98. Tee YH, Shemesh T, Thiagarajan V, Hariadi RF, Anderson KL, Page C, et al. Cellular chirality arising from the self-organization of the actin cytoskeleton. Nature cell biology. 2015;17(4):445–457. pmid:25799062
View Article
PubMed/NCBI
Google Scholar

[340] View Article

[341] PubMed/NCBI

[342] Google Scholar

[ref99] 99. Weirich KL, Banerjee S, Dasbiswas K, Witten TA, Vaikuntanathan S, Gardel ML. Liquid behavior of cross-linked actin bundles. Proceedings of the National Academy of Sciences. 2017;114(9):2131–2136. pmid:28202730
View Article
PubMed/NCBI
Google Scholar

[344] View Article

[345] PubMed/NCBI

[346] Google Scholar

[ref100] 100. Weirich KL, Stam S, Munro E, Gardel ML. Actin bundle architecture and mechanics regulate myosin II force generation. Biophysical Journal. 2021;120(10):1957–1970. pmid:33798565
View Article
PubMed/NCBI
Google Scholar

[348] View Article

[349] PubMed/NCBI

[350] Google Scholar

Figures

Abstract

Author summary

Introduction

Materials and methods

Semiflexible, inextensible fibers

Dynamic cross-linking model

Fiber turnover

Mobility evaluation

SAOS rheology

Temporal integration and splitting

Results and discussion

Fiber turnover creates dynamic steady states with varied degrees of bundling

Homogeneous meshwork morphology.

Bundles embedded in meshwork morphology (B-In-M).

Other systems (parameter variations).

The linear response regime (LRR): Viscoelastic moduli and stress spectra.

Timescales of stress relaxation

Viscoelastic moduli

Short timescales: Elastic links and viscous fibers

Long timescales: Viscosity depends on morphology and link unbinding rate

Intermediate timescales

Role of nonlocal hydrodynamic interactions

Rescaling of time cannot explain long-timescale moduli.

Bundled fibers: When flow reduces stress.

A continuum framework: The generalized Maxwell model

Conclusion

Supporting information

S1 Fig. Dynamic steady state for τf = 5 seconds.

S2 Fig. Dynamic steady state for τf = 10 seconds.

S3 Fig. Homogeneous meshwork for smaller bending stiffness κ.

S4 Fig. Spectrum of the stress for the B-In-M network.

S5 Fig. Maxwell model fit for the B-In-M system with full hydrodynamics and a varying number of elements.

S6 Fig. Maxwell model fit for the system with ten times larger viscosity and local drag mobility.

S1 Video. Homogeneous meshwork system with Ld = 2 μm in its dynamic steady state.

S2 Video. Bundle-in-mesh geometry with Ld = 3 μm in its dynamic steady state.

S3 Video. Stress relaxation test for the homogeneous meshwork geometry.

S4 Video. Stress relaxation test for the bundle-in-mesh geometry.

S5 Video. Oscillatory shear rheology experiment for the homogeneous meshwork geometry.

S6 Video. Oscillatory shear rheology experiment for the bundle-in-mesh geometry.

Acknowledgments

References

S1 Fig. Dynamic steady state for τ_f = 5 seconds.

S2 Fig. Dynamic steady state for τ_f = 10 seconds.

S1 Video. Homogeneous meshwork system with L_d = 2 μm in its dynamic steady state.

S2 Video. Bundle-in-mesh geometry with L_d = 3 μm in its dynamic steady state.