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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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):

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.
2. 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.
3. 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.
4. 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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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:

1. 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.
2. 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.
3. 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.
4. 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

1. Turnover filaments to update configuration .
2. Process binding/unbinding events to evolve Y_n → Y_n+1.
3. 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)).

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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

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