Computational Analysis of Viscoelastic Properties of Crosslinked Actin Networks

Mechanical force plays an important role in the physiology of eukaryotic cells whose dominant structural constituent is the actin cytoskeleton composed mainly of actin and actin crosslinking proteins (ACPs). Thus, knowledge of rheological properties of actin networks is crucial for understanding the mechanics and processes of cells. We used Brownian dynamics simulations to study the viscoelasticity of crosslinked actin networks. Two methods were employed, bulk rheology and segment-tracking rheology, where the former measures the stress in response to an applied shear strain, and the latter analyzes thermal fluctuations of individual actin segments of the network. It was demonstrated that the storage shear modulus (G′) increases more by the addition of ACPs that form orthogonal crosslinks than by those that form parallel bundles. In networks with orthogonal crosslinks, as crosslink density increases, the power law exponent of G′ as a function of the oscillation frequency decreases from 0.75, which reflects the transverse thermal motion of actin filaments, to near zero at low frequency. Under increasing prestrain, the network becomes more elastic, and three regimes of behavior are observed, each dominated by different mechanisms: bending of actin filaments, bending of ACPs, and at the highest prestrain tested (55%), stretching of actin filaments and ACPs. In the last case, only a small portion of actin filaments connected via highly stressed ACPs support the strain. We thus introduce the concept of a ‘supportive framework,’ as a subset of the full network, which is responsible for high elasticity. Notably, entropic effects due to thermal fluctuations appear to be important only at relatively low prestrains and when the average crosslinking distance is comparable to or greater than the persistence length of the filament. Taken together, our results suggest that viscoelasticity of the actin network is attributable to different mechanisms depending on the amount of prestrain.


Introduction
Actin is the most abundant intracellular protein in eukaryotic cells and plays an important role in a wide range of biological and mechanical phenomena [1]. Monomeric actin (G-actin) selfassembles to a filamentous form, F-actin, which is crosslinked into the actin cytoskeleton by various actin crosslinking proteins (ACPs). It has been known that mechanical force plays a crucial role in the physiology of eukaryotic cells [2], and therefore appropriate functions of living cells are attributable to the rigorous control of their rheological properties [3]. Thus, investigating rheological properties of actin networks is indispensable for elucidating the mechanics of cells as well as for understanding a wide variety of cellular processes. Experiments have been conducted to probe viscoelastic properties of cells and reconstituted actin gels using a variety of techniques such as microbead rheology, magnetic bead cytometry, and bulk rheology [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. In experiments, discrepancies have been observed among measurements using dissimilar methodologies, and many of the observed features are not well understood. For example, viscoelastic moduli measured by single-bead passive microbead rheology are much smaller than those determined by 2-point microrheology or bulk rheology [4][5][6]. Also, although distinct power law responses of the storage modulus have often been observed in vivo and in vitro [7][8][9][10][11], their origin is not yet clearly understood.
Concurrently, characteristics of semi-flexible polymer networks have been studied theoretically and computationally [24][25][26][27][28][29][30][31][32][33][34]. Two- [24,25] and 3-dimensional computational models [26] studying affine and nonaffine deformations of semi-flexible networks responding to large shear strain revealed two regimes dominated by bending or stretching of filaments, respectively. Recently, using a microstructure-based continuum mechanics approach, Palmer and Boyce reproduced many of the rheological properties of actin networks observed in experiments [27]. The viscoelastic behavior of semi-flexible networks was also investigated using dissipative particle dynamics and the concept of microbead rheology, where scale-free behavior of the bead displacement was observed [28,29]. To date, however, most of these models neither explicitly take into account ACP mechanics nor systematically account for thermal fluctuations, nor have they been used to explore the effects of finite prestress on viscoelasticity, all of which are potentially important factors governing matrix viscoelasticity.
With the objective of extending these previous works and providing new insights into underlying mechanisms, we develop a Brownian dynamics (BD) model of the actin network that includes features such as steric interaction among filaments, the usage of explicit crosslinkers, a more realistic morphology, and the consideration of crosslinker stiffness. By measuring stress in response to applied oscillatory shear strain (''bulk rheology'') and thermal fluctuations of individual segments in the polymeric chain (''segment-tracking rheology''), we investigate viscoelastic properties of actin-like networks. Throughout this paper, for convenience, the term, ''actin network'' is used to refer to the network being simulated. It should be noted, however, that some of the properties employed in our model, especially for ACPs, were estimated since they are not well-known experimentally. Due to simplifications in the model and parameter uncertainty, the results should therefore be viewed as representative of a generic crosslinked network, but lack a quantitatively precise correspondence to actin networks.
Nevertheless, we found features that semi-quantitatively capture experimentally observed behaviors of actin networks. The storage and loss moduli, G9 and G0, followed power laws as functions of the oscillation frequency. As the prestrain increased, the network became increasingly elastic. Bending and extensional stiffnesses of actin filaments and ACPs played an important role depending on the degree of prestrain. We found that the mechanical response of the network is dominated by a percolating 'supportive framework,' while other actin filaments contribute little to the viscoelastic moduli. Surprisingly, in typical physiological conditions where the distance between crosslinking points along F-actin is much shorter than the actin persistence length, we found that thermal fluctuation plays little role in viscoelasticity, so that the network consisting of crosslinked F-actins can be viewed essentially as a deterministic overdamped system in a viscous medium. In sum, our computational model elucidates how various mechanical responses (thermal forces and the bending and stretching of actin filaments and ACPs) govern viscoelastic properties of the network under different conditions.

Model overview
To create an actin network, we began with a uniform distribution of actin monomers with ACPs dispersed randomly, and allowed the network to polymerize until 99% of the monomers were incorporated into filaments [35]. Motion of the monomers, filaments, and ACPs followed Brownian dynamics, and bonding occurred according to a first order irreversible process. After polymerization, we applied a coarsegraining procedure explained in Methods to cover a larger system size. In the resulting networks, the filament length distribution and the network morphology [35] bear a closer resemblance to reconstituted actin gels ( Figure 1) than those generated by the random placement of equal-length filaments [24][25][26]29,32]. F-actins and ACPs in our model are characterized by introducing bending and extensional stiffnesses ( Table 1). Two types of ACPs were used, depending on whether they form bundles (ACP B ) or orthogonal networks (ACP C ) during the polymerization process with equilibrium crosslinking angles of 0 and p/2, respectively, between the two crosslinked actin filaments.
Viscoelasticity of the network was probed in two ways: 1) In segment-tracking rheology, we analyzed the random thermal motion of many individual actin segments in the network [36]. 2) In bulk rheology, we applied an oscillatory shear strain on the top surface of the system while holding the bottom surface fixed. The measured viscoelastic moduli may depend on detailed geometry and the extent of percolation, especially if the network is small. Therefore, the use of a geometrically identical network for all simulations enables us to systematically control and isolate the effects of a given parameter. We thus used the two representative networks shown in Figure 1 for our measurements. The filament length (L f ) was 1.5 mm60.65 mm (average6standard deviation), and the actin concentration (C A ) was 12.1 mM. We randomly removed ACPs from the networks to change R, the ratio of ACP concentration to C A (parameters are defined in Table 1). Actin filaments longer than the width of the simulation domain (2.8 mm) were severed to reduce finite size effects.

Comparison to experiments
While recognizing the limitations of any direct quantitative comparisons between our model predictions and experiments, we conducted one set of experimental measurements under conditions similar to those of the simulation (mean filament length ,L f . = 1.5 mm . , R = 0.01, and C A = 12.1 mM) and compared viscoelastic moduli. In order to match ,L f ., gelsolin was added to the sample; the length distribution was determined by fluorescence imaging. One difficulty in matching simulation and experimental conditions originates from the fact that R in the experiment corresponds to the total amount of ACPs in the sample, including both those in active and inactive (partially bound or free) states, whereas in our simulation, it indicates the net amount of active ACPs that crosslink or bundle two filaments. In addition, due to computational constraints, the oscillation frequency tested in the simulation overlaps that of the bulk rheology experiment only in a narrow range. Despite these difficulties, as seen in Figure 2, the values of G9 and G0 computed using the model are in reasonable agreement with the experiment, both qualitatively and quantitatively, perhaps better than one might have expected given the range of uncertainty of some of the parameters.

Author Summary
The actin cytoskeleton provides structural integrity to a cell, is highly dynamic, and plays a central role in a wide variety of phenomena such as migration and the sensation of external forces. For years, researchers have studied the mechanics of the cytoskeleton by creating actin gels in the laboratory in combination with proteins that bridge between and reinforce the actin gel found inside cells. These gels, however, failed to replicate many aspects of cell behavior. Recent studies have shown that tension within the cytoskeleton contributes to the observed stiffness of cells. Still, our understanding of cytoskeletal mechanics is incomplete, and many observed phenomena cannot be explained by existing models. Here, we simulate a three-dimensional network containing actin filaments linked together by other proteins. We studied the relative contributions of thermal fluctuations of the network and the stiffness of filaments and linking proteins. Under conditions that replicate those in a cell, properties of the linking proteins are surprisingly significant, as is the stiffness of the actin filament to stretching. Thermal fluctuations are relatively unimportant, but become more so at low levels of resting tension. At high tensions, a small fraction of filaments support a majority of the load.

Effects of the concentration and type of ACPs on G9 and G0
We computed G9 and G0 of structures crosslinked via ACP C (R = 0, 0.01, and 0.021) and those bundled via ACP B (R = 0, 0.01, 0.02, and 0.04) using both bulk rheology and segment-tracking rheology. Networks without ACPs exhibit a slope of G9 close to 0.75, as indicated by the black solid line in Figures 3A and 4A. This value has been observed in various experiments [8,[37][38][39], and it is known to originate from transverse thermal undulations of actin filaments [30]. Interestingly, when repulsive forces between the filaments are eliminated, the slope of G9 estimated via segmenttracking rheology approaches unity (data not shown), implying that the volume exclusion effect of neighboring filaments creates a tube-shaped space that hampers free translation and rotation of the filament [31]. Although the filament confined in this way can perform reptation, it is operative on long time scales that are beyond those attainable in these simulations, so the mean square displacement (MSD) observed here primarily reflects transverse thermal motions, resulting in the 0.75 slope. On the other hand, the plateau in G9 often exhibited by experiments [6,12,40] was not observed within the frequency range of these simulations. Note that the plateau modulus is induced by entanglement effects that become more pronounced at longer time scales. The combination in these simulations of relatively short filaments leading to longer entanglement time [12] and computational constraints precluding simulations for longer times limited our ability to observe a plateau.
Values of viscoelastic moduli attained using bulk rheology and segment-tracking rheology exhibit surprisingly good agreement ( Figure 3) even though segment-tracking rheology (Equation 11) was originally developed for a test particle much larger than the meshwork of filaments [36]. ACP C elevates the magnitude of G9 and reduces its slope ( Figure 3A), implying that the frequency dependence of G9 is reduced as networks incorporate more ACP C . G9 follows a power law, G9,f s 0.3 Figure 1. Two representative networks used in this study. (A) A network bundled via ACP B and (B) crosslinked via ACP C , both of which consist of actin filaments of various lengths (cyan) and ACPs (red). For visualization, VMD was used [53]. Two arms of each ACP drawn with partial red and cyan are connected to filaments, forming crosslinks. In (A), due to the small computational domain, ladder-like structures are predominant rather than long, thick bundles. The linear dimension of the simulation box is 2.8 mm. These two networks form the basis for all simulations; networks with low R were obtained from these by eliminating a portion of active ACPs. The inset of each plot shows the detailed geometry of bundled or crosslinked structures consisting of actin filaments and ACPs. doi:10.1371/journal.pcbi.1000439.g001 Numbers in parentheses are corresponding dimensionless values as defined in the text. The value and notation of other parameters are the same as in [35].
(dashed line), for f s ,100 Hz at the highest crosslink density, R = 0.021, which is within the range of powers observed in cells, 0.15-0.3 [9][10][11]. However, the magnitude of G9 from these simulations was much lower than in vivo values. This is likely due to many factors, notably the absence of prestrain. G0 increases slightly as the amount of ACP C is increased, but the slope remains similar ( Figure 3B). In addition, a decrease in the phase delay, tan 21 (G0/G9) (Equation 10), accompanies the increase in R, indicating that ACP C elevates the elasticity of the network. At f s ,10 3 Hz, the phase delay depends only weakly on R, but as frequency decreases, it decreases more quickly with higher R, implying a greater effect by crosslinking at lower frequencies.
Networks bundled by ACP B exhibit a behavior distinctly different from that of ACP C (Figure 4). Large differences in G9 and G0 were observed between segment-tracking rheology and bulk rheology. This originates from the heterogeneity of the bundled network attributable to the small computational domain, for which viscoelastic moduli measured by bulk rheology depend strongly on whether or not there are bundles that percolate between the top and bottom boundaries. We thus discuss results of only segment-tracking rheology for networks formed by ACP B .
ACP B increases G9 but has little effect on its slope in contrast to ACP C ( Figure 4A), whereas the phase delay is only slightly influenced by R. ACP C is able to form a well-percolated network even at relatively low R, so that the network gels more efficiently [14] and acts like a single-body elastic object in response to shear stress or strain. By contrast, ACP B bundles filaments together, resulting in the relatively low level of percolation for the same value of R as reflected by its low connectivity [35]. Although the diffusivity of bundled filaments is lower, the lack of connectivity leads to the absence of elastic behavior at long time scales. In order for ACP B to increase network elasticity to a similar extent by ACP C , its concentration may have to be much higher.

Effects of prestrain on G9 and G0
Only bulk rheology is employed here since the application of prestrain leads to a highly nonuniform distribution of the load, with a small fraction of highly tensed filaments and a larger number of filaments under little or no stress. Due to such heterogeneity, segment tracking rheology underestimates G9 and G0 as it randomly traces N MSD segments from the entire network. The effect of prestrain on G9 is analogous to that of ACP C . As seen in Figure 5A, G9 increases and produces a weaker dependence on frequency at higher prestrain (c); at c = 0.55, G9 is virtually independent of frequency and is nearly 100-fold larger than that at c = 0 for f s ,10 Hz. This means that large prestrain transforms the network into a highly elastic one that exhibits a phase delay close to 0 at all frequencies and results in G9 comparable to in vivo values [9][10][11], wherein it has been postulated that prestress or prestrain plays a significant role [14]. G0 also exhibits interesting behavior; at high prestrain, it increases slightly at low frequency, similar to in vitro observations using heavy meromyosin (HMM) [18,41], and a similar increase in G0 was also observed in vivo [19,38]. This  suggests that at low frequencies, viscous effects play an important role. Tharmann and coworkers argued that this trend of G0 may be due to the unbinding of HMM. However, since unbinding was not permitted in these simulations, the increase in the lowfrequency G0 with prestrain must originate from a different mechanism. Note that such a tendency may have been more evident if the simulations were capable of reaching even lower frequencies.
Also, the relation between G9 at f s = 3.16 Hz and prestress (t 0 ) was investigated. It remains relatively constant until a threshold prestress (t 0 ,0.1 Pa, Figure 5B) beyond which it increases following a power law, G9,t 0 0.85 . The exponent of 0.85 is close to the value of ,1 found in in vitro experiments under similar conditions [14].

Effects of extensional stiffness of actin filaments, k s,A
To illustrate how prestrain transforms a network into a more elastic one, we display the network using a color scale depending on bond length averaged for duration of 0.1 ms. Only a small number of actin filaments aligned in the x-z direction are highly stretched ( Figure 6A,B). As mentioned above, this heterogeneity precludes using segment-tracking rheology which measures thermal motions of randomly selected segments, many of which are not a part of the highly stretched filaments in prestrained networks. The mean filament length of the entire network increased by only 0.5% with c = 0.55. However, due to the large value of k s,A (Table 1), this results in large spring forces that contribute significantly to the high magnitude of G9.
We also performed simulations with different extensional stiffness (k s,A = 0.0338 and 0.0068 N/m) for networks with c = 0.55 ( Figure 6C). G9 and G0 decreases with lower k s,A , but its phase delay is virtually unchanged (data not shown). When c = 0, however, variation in k s,A has little or no effect since most actin filaments are not highly stretched ( Figure 6D). Therefore, k s,A affects viscoelasticity only under high prestrain.  We also studied the influence of k s,A on G9 (f s = 3.16 Hz) at different levels of prestress ( Figure 6E). Interestingly, when k s,A is reduced to 0.0068 N/m, the previous slope of 0.85 decreases to about 0.7, suggesting that the slope of the curve in the nonlinear regime increases with k s,A . Note that we used a value of k s,A that is 1/40 of the experimentally measured value due to computational efficiency. Our result is thus consistent with a larger slope observed in experiments, ,1 [14].

Effects of bending stiffness and thermal fluctuation of actin filaments
To address factors affecting viscoelasticity at low prestrain, we probed the influence of actin filament bending stiffness, k b,A , by using three different values (k b,A = 1.056610 218 , 1.056610 219 , and 1.056610 220 Nm) at c = 0 ( Figure 7A) and 0.4 ( Figure 7B). At high prestrain, c = 0.4, variations in k b,A have little effect on G9 and only minor effects on G0, suggesting that in the high prestrain regime, actin bending stiffness does not play a major role in viscoelastic properties ( Figure 7B). By contrast, at c = 0, k b,A influences G9 ( Figure 7A), although perhaps not to the extent that one might expect for reductions in k b,A by factors of 10 and 100.
We also studied the effect of thermal fluctuation of actin filaments using separate simulations with or without the Langevin force term,F i B (Equation 1), at various c and l p (,k b,A /k B T, with T = 300 K) ( Figure 7C). Without F i B , the temperature of the system drops to ,10 K, as measured by applying the equipartition theorem, where the nonzero temperature is due to the externally imposed oscillation. Figure 7C shows the ratio of G9 without thermal fluctuation (TF) to that with TF at f s = 10 Hz under each condition. At high prestrain (c$0.4), the elimination of TF does not affect G9 regardless of l p . Surprisingly, for filaments with bending stiffness comparable to that of actin (l p = 10-20 mm), elimination of thermal effects had no significant effect on G9 at any value of c. Only when k b,A was decreased (l p #3 mm) did thermal fluctuations have a noticeable influence on G9. Under these conditions, l p became comparable to the average distance between crosslinks, l c , (0.393 mm at R = 0.021), in which entropic effects would be expected to play a role [15], especially with a low prestrain. However, it is not clear whether this range of conditions can be attained in vivo. This finding that the relation between l p and l c determines the importance of thermal fluctuations is consistent with previous qualitative predictions [24]. Yet, our result is at odds with a previous view that entropic effects due to thermal undulations of individual filaments are responsible for the elasticity of scruin-crosslinked [15,16] and HMM-crosslinked [18] networks even in the case when l c %l p . In previous numerical studies [25,26], thermal fluctuations were applied only before applying shear strain, not during the measurement of stress, and thus they cannot be used as a clear demonstration that thermal fluctuations are important. Further experimental, numerical, and theoretical investigations are necessary to clarify the role of thermal fluctuation and entropic elasticity. For instance, if network elasticity is mainly governed by enthalpy, and if thermal fluctuation plays little role, we expect that G9 will be minimally affected upon adding crowding agents to increase solvent viscosity and to provide steric barrier to conformational motion. However, crowding agents may affect organization of filaments, such as enhanced bundling [42]. Thus, such experiments will need to be carefully interpreted. Also, the vibrational motion of actin filaments depending on l c measured in a computational study can be compared to theoretical predictions in [24].

Effects of bending and extensional stiffnesses of ACPs
Previous studies showed that the detailed structure of ACPs strongly influences viscoelastic moduli of actin networks. In [20], each synthetically constructed crosslinking molecule of a different length produced distinctive macroscopic mechanical behaviors. In addition, Gardel et al. observed significantly different stress-strain relationships between a reconstituted actin network with intact filamin A (FLNa) and one with mutated hingeless FLNa [14].
Since these previous results led us to anticipate substantial effects of ACPs in our actin networks, we investigated the influence of bending stiffness of ACP C , k b,ACP,1 and k b,ACP,2 . As described in Methods, k b,ACP,1 acts between two arms of the ACP, whereas k b,ACP,2 limits bending between one arm of the ACP and the axis of actin filament to which it is attached, with an equilibrium value of p/2. To study their effects, we decreased both k b,ACP,1 and k b,ACP,2 10 fold, 100 fold, and to zero. This had a significant effect on G9 and G0 at all prestrains (Figure 8), but especially so for the magnitude of G9 in the highly prestrained case with c = 0.4. At small prestrain (c = 0), deformation is mostly associated with the bending of actin filaments, and the ACP angles remain close to their equilibrium values. As prestrain increases (to 0.4), the actin filaments are progressively straightened, especially those that support the bulk of the load. Since this is accompanied by changes in the angle between crosslinked actin filaments, bending stiffness of ACPs becomes an important determinant of G9. As prestrain further increases (to 0.55), ACPs are maximally bent, and stretching of actin filaments is the dominant mechanism for resisting deformation since changing the extensional modulus has the greatest influence on G9. It should be noted, however, that though the geometry of ACP C mimics that of FLNa, the large values of k b,ACP are closer to that of stiff scruin. Studies involving mutant ACP C (e.g. FLNa) that differs in bending stiffness would further clarify our observation.
A change in extensional stiffness of ACP C , k s,ACP , results in little change in G9 and G0 (data not shown) at low c, but plays a significant role at high c. This tendency is not surprising as unbending of the V-shaped ACP C should precede stretch of ACP C arms responding to the load.

Significance of each parameter at various prestrains
Based on the effects of bending and extensional stiffnesses of actin filaments and ACPs as well as thermal fluctuations discussed above, we can estimate the relative importance of each factor over a wide range of prestrain and identify regimes where different phenomena dominate the viscoelastic behavior. For this purpose, each stiffness was decreased by 25-fold from the standard case, and the ratio of G9 with the decreased stiffness to that with the normal stiffness was computed at f s = 10 Hz (Figure 9). A fall in extensional stiffness of both actin filaments (k s,A ) and ACPs (k s,ACP ) decreases G9 more as c increases, implying that the stretch of filaments and ACPs plays an important role in strain-stiffening behavior at high c, at least at the exaggerated levels of extensional compliance used in the simulations. Bending stiffness of ACPs, k b,ACP , is significant at all tested c, but has the largest effect on G9 for c,0.3-0.4. The influence of bending stiffness of actin filaments, k b,A , on G9 is interesting, passing through a minimum at c = 0.2. We measured G9 under the same conditions but without thermal fluctuation, and the minimum at c = 0.2 disappeared; the ratios at c = 0-0.2 are similar to each other. Therefore, thermal fluctuation of actin filaments contributes to an increase in G9 at low c when l p is comparable to l c .
Based on these observations, we propose three distinct regimes: i) A low c regime where k b,A is dominant, with thermal fluctuations playing a substantial role for l p #l c . ii) An intermediate c regime in which the effect of k b,ACP becomes dominant. iii) A high c regime where k s,A and k s,ACP are the predominant factors. Transition from one regime to another, however, is not sharply defined. Others [25,26] have also argued that bending stiffness of filaments would dominate in the low c regime, and extensional stiffness of filaments in the high c regime, but they did not consider stiffness of ACPs as a parameter. Also, since these previous simulations lacked thermal fluctuations, their relative effects could not be assessed. Figure 6B shows that only a subset of the entire network supports the dominant portion of the stress in prestrained networks. While there could be many different criteria for identifying such a 'supportive framework,' we considered two, based either on the stretch of actin filaments or on the bending of ACPs. In the first case, we deleted filaments that were not highly stretched. The remaining network, however, when oscillated at the same strain amplitude, produced very low stress levels, suggesting that a significant fraction of the stress-supporting elements had been removed by this process. In the second case, we use bending force on ACPs as a selection criterion since the change in k b,ACP had a strong effect on viscoelastic moduli ( Figure 8). First, the sum of magnitudes of all bending forces applied on each ACP in a prestrained state was calculated during 0.1 ms, the portion of ACPs with the highest bending forces were selected, and the remaining ACPs were deleted. All actin filaments not connecting between a pair of these highly strained and bent ACPs were removed. We compared stresses exerted on networks in which 25%, 50%, and 75% of ACPs remain under the same imposed prestrain. The network containing only 25% of initial ACPs ( Figure 10A) consists of filaments oriented mostly in the diagonal direction on the x-z plane (cf., Figure 6B). Moreover, as the orientational distribution shows ( Figure 10B, inset), there are (less stretched) filaments oriented in other directions that also help to transmit the applied load, by bending of connecting ACPs, which are not selected when a criterion based on stretch of actin filaments is used. The mean stress and oscillation amplitude of the reduced network remain nearly at original levels ( Figure 10B). For example, the network containing only ,28% of the original actin filaments and ,25% of ACPs produced ,70% of the original mean and amplitude of oscillating stress at the same levels of strain. Thus, bending force on ACPs is a major determinant of the elasticity of prestrained networks. We also confirmed that at high prestrain, ACPs experiencing large bending forces tend to displace in a manner consistent with the deformations being affine, as measured by the S parameter used in [43] (data not shown). This means that parts of the supportive framework that bear greater levels of force experience more affine deformation.

The supportive framework governing viscoelastic moduli
It is well known that most cells are in a prestressed state largely due to the actomyosin contractile apparatus [9,44,45]. By analogy, we hypothesize that the large G9 (,10 3 Pa) measured in cells is mostly due to the supportive framework composed of the contractile apparatus while the rest of cytoskeletal structures play comparatively little role in the macroscopic viscoelasticity of cells.

Limitations resulting from permanent crosslinks
Most ACPs have finite binding lifetimes leading to dissociation and reformation of crosslinks, which is known to play an important  role in the viscoelastic moduli and rheology of actin networks [22,46,47]. However, ACPs in our present simulation do not undergo dynamic rearrangement. This artifact renders the level of stress responding to prestrain unreasonably high. While stress catastrophically drops beyond 1,30 Pa in experiments [14,16,20,48], it can increase over 100 Pa in our simulation. The absence of crosslink dynamics causes G9 in Fig. 5A to be much higher and to exhibit a lower power compared to experiments. Nevertheless, our present results provide valuable insights as well as benchmarks whereby stress relaxation or dynamic reorganization of networks can be compared as an extension of this work.

Conclusion
Using Brownian dynamics simulations, we systematically investigated viscoelastic moduli of actin-like networks. First, viscoelastic moduli measured in our model compared favorably with those of experiments conducted using a bulk rheometer under similar conditions. Then, effects of the kind and concentration of ACPs on G9 and G0 were investigated. With no ACP, G9 exhibited a power law slope of 0.75 against frequency, reflecting the transverse thermal motion of actin filamets. ACP C tended to produce networks that were more elastic with a weaker dependence of G9 on frequency and smaller phase delay than ACP B . The influence of prestrain on G9 and G0 was also studied. High prestrain makes the network more elastic mainly via the stretch of a small fraction of aligned actin filaments. In addition, the effects of extensional and bending stiffnesses of actin filaments and ACPs on viscoelasticity were tested individually. We found that for viscoelastic moduli, three regimes are evident, each governed by different mechanisms depending on the amount of prestrain. At low prestrain, bending stiffness of ACPs and actin filaments plays an important role in determining viscoelasticity. At intermediate prestrain, bending stiffness of ACPs affects viscoelas-tic moduli the most. At high prestrain, the extensional stiffness of actin filaments and ACPs becomes dominant. We also found that entropic effects due to thermal fluctuations of filaments are important only when the prestrain is low and the average distance between active ACPs is comparable to the filament persistence length. Last, we identified the supportive framework that largely accounts for the high elasticity of prestrained networks, as that associated with strained ACPs rather than stretched actin filaments. Even after 75% of the network components are deleted, stress remains at 70% of the original value.
Our computational model using the discrete network of crosslinked F-actins provides insights into the microscopic origin of the viscoelastic behavior of the actin cytoskeleton. Importantly, our findings regarding the limited role of thermal fluctuation, as well as the different contributions to viscoelastic responses by the bending or stretching of F-actin or ACPs, require additional experiments for further investigation. The present computational framework can be further developed to study phenomena associated with ACP unbinding, such as stress relaxation, network reorganization, and plastic deformation by shear.

Methods
Our Brownian dynamics (BD) model extends our previous work in which we investigated the morphology of polymerized and crosslinked actin networks [35]. In this simulation, actin monomers, filaments, and ACPs experience thermal motion and interact with each other with defined potentials and binding probabilities. Individual displacements are governed by the Langevin equation: Figure 10. The supportive framework bearing most of stress. (A) Network composed of filamentous actin connected via ACPs that support the highest 25% of ACP bending forces. In contrast to Figure 6B, there are filaments that are almost perpendicular to the diagonal direction on the x-z plane, which are not highly stretched yet transmit load by bending of ACPs connected to them. (B) Stress exerted by prestrained networks (c = 0.4) consisting of a fraction of actin filaments and ACPs. The extent of ACP bending forces is employed as a criterion to retain elements. Each symbol corresponds to a different percentage ratio of the number of ACPs remaining in a rebuilt network to that in the original network: 100% (black circles), 75% (magenta triangles), 50% (blue inverted triangles), and 25% (green diamonds). The fraction of remaining actin segments is, respectively: 100%, 79%, 52%, and 28%. (inset) Orientation angles of actin segments projected onto the x-z plane for the network in Figure 10A. Segments oriented in the z direction have a value of 0u. Most actin segments in the reduced structure are oriented in the (+x)-(+z) direction (45u), but segments with other orientations are also important, presumably because they transmit stress through bending of the ACPs attached to them. doi:10.1371/journal.pcbi.1000439.g010 where m i is the mass of ith element (actin or ACP), r i is the element's location, f i is the friction coefficient, t is time, and v ' is the velocity of the surrounding medium. F i is a net deterministic force, and F B i is a stochastic force satisfying the fluctuation-dissipation theorem. Considering that inertia of all elements is negligible on the length and time scales of interest, the Langevin equation is simplified by setting the acceleration term to zero and cast in nondimensional form using k B T, f C,A , and L C,A as primary variables, where f C,A and L C,A are the friction coefficient and length of a segment of F-actin (as explained later). Positions of all elements are updated using the Euler integration scheme: where dimensionless variables are indicated by the tilde '' , ''. Two interaction potentials describe bond stretch and bending with stiffnessk k, denoted by subscripts ''s'' and ''b'', respectively [35]:Ũ wherer r 12 is the distance between two points constituting a chain, h is the bending angle, and the subscript 0 denotes the equilibrium value.

Coarse-graining scheme using cylindrical segments
In our previous model, we treated a segment of F-actin as a spherical particle representing two G-actins. To simulate larger length and time scales, we introduced coarse-graining in which a cylindrical segment represents several monomers. We kept thermal forces on ACPs and actin segments in a form similar to our previous model, but incorporated a cylindrical geometry for calculating repulsive forces.
As seen in Figure 11, the points of two adjacent elements on a filament correspond to the ends of one cylindrical segment representing N C actin monomers of the previous model. Accordingly, the diameter of the cylindrical segment, s C,A , is the same as the diameter of actin monomer of the previous model, s A = 7 nm, and its length, L C,A , is N C ?s A . By letting each monomer of the polymerized network evolve into a cylindrical segment, a larger network is created, albeit with a lower concentration of ,1-100 mM. In addition, one point on a filament axis and the other point representing the center of an ACP determine the two end points of each arm of the ACP, indicating that ACPs have a structure with two cylindrical arms whose length and diameter are L C,ACP and s C,ACP respectively. In this study, N C was set to 10; this degree of coarse-graining is appropriate since the length of one rigid cylindrical segment, L C,A = 70 nm, is still much shorter than the persistence length of an actin filament, ,10 mm.
One consequence of this coarse-graining technique is that it increases the arm length of ACPs. That is, since the redundant volume of the original monomers is neglected, the arm length of ACPs is extended by (L C,A 2s C,A )/2. This produces a network for which the shortest distance between two bundled filaments is (L C,A 2s C,A )/2 (63 nm in this case) that exceeds the typical spacing formed by many of the bundling ACPs (e.g. fascin, fimbrin, and aactinin). However, it is not expected that this will significantly alter the qualitative behavior of bundled networks since the change in the ACP arm length is much shorter than the length of F-actin. Previous experiments where different bundling ACPs led to distinct macroscopic behaviors [21,22] are more likely due to the different stiffness and dissimilar binding affinities of the ACPs rather than their physical size.
Repulsive forces are computed according to the following harmonic potential depending on the minimum distance,r r 12 , Figure 11. Schematic diagrams explaining the geometry of a network and the calculation of repulsive forces. (A) A coarse-graining scheme using cylindrical segments with N C = 5. Dashed lines show monomers and ACPs used to generate the network using the polymerization model [35]. Once the network is formed, it is coarse-grained by replacing the original spheres by cylinders as shown. wheres s C,i,1 ands s C,i,2 are diameters of the cylindrical segments (i = A or ACP), andk k r is the strength of repulsive effects. Forces calculated from the potential (Equation 4),F F r , are distributed onto the two end points constituting the cylindrical segment, a and b, via the following equations: where y represents the location on the segment at which the minimum distance, r 12 , is measured ( Figure 11B; 0#y#L C,i ).

Polymerization and crosslinking
Monomer assembly and disassembly are important determinants of the morphology of actin networks. Given the rate constants obtained in recent in vitro experiments [49], it would take ,100 s for an actin filament 1.5 mm in length to completely depolymerize, or ,1-10 s to polymerize the same filament with C A = 12.1 mM. Though in vitro depolymerization of F-actin is very slow, various ACPs can accelerate it in vivo. However, we assume here for simplicity that neither polymerization nor depolymerization of actin filaments occurs within the time scale of interest in this study, ,1 s. Unlike our previous model [35], an ACP can bind to any point along an actin filament in any circumferential direction.

Preparation of a network to estimate viscoelastic properties
The measurement of viscoelastic moduli can be influenced by the detailed geometry and the extent of percolation, especially if the network is small. Therefore, the use of a geometrically identical network for all simulations enables us to systematically control and isolate the effect of a given parameter. In other studies, actin networks have been generated by the random placement of equal-length filaments [24][25][26]29,32]. However, we generated more realistic networks using the previous model that incorporates both polymerization and crosslinking. A somewhat heterogeneous network bundled by ACP B ( Figure 1A) and a well-percolated, homogeneous network crosslinked by ACP C ( Figure 1B) were prepared. In Figure 1A, ladder-like structures consisting of two actin filaments and multiple ACP B are evident, in contrast to thick bundles that are often observed in experiments. The relative absence of thicker bundles in these simulations is attributable to the small domain size compared to an average filament length. With N C = 10, the filament length (L f ) is 1.5 mm60.65 mm (average6standard deviation), and the actin concentration (C A ) is 12.1 mM. We randomly removed ACPs to change R (Table 1), while maintaining the overall network geometry. In addition, actin filaments longer than the width of the simulation domain (2.8 mm) were severed to minimize finite size effects.

Adjusting parameter values for coarse-graining
Upon coarse-graining, several parameters need to be adjusted. For the cylindrical geometry of the segments, approximate forms for friction coefficients are [50]: f C,i,k~3 pgs C,i : 4zL C,i =s C,i 5 where g is the viscosity of the surrounding medium. However, since the mostly crosslinked cylindrical segments move predominantly in the transverse direction, and considering that f C,A,\ is only 1.64 times higher than f C,A,k L C,A =s C,A~1 0 ð Þ , f C,A,\ was used in all directions for simplicity. To test this assumption, we compared simulations with and without the anisotropic friction coefficient at zero prestrain. At f s = 10 Hz, we obtained G9 = 5.07 Pa (isotropic) and 5.00 Pa (anisotropic), and G0 = 2.02 Pa (isotropic) and 2.39 Pa (anisotropic). Since motion along the filament axis is more suppressed with prestrain, the effect of anisotropic friction coefficient will be even less. We also assumed a constant friction coefficient regardless of the filament length to which the segment belongs. Hydrodynamic interactions between filaments are expected to play little role and were ignored since the volume fraction of actin is low (,0.1%), and because the filaments have high aspect ratio (long thin rod) [51].
Since the geometry of ACPs changes after coarse-graining, the following equilibrium values for additional harmonic potentials were used for ACP B : k k s,ACP~5 ,000,k k b,ACP, 1~2 50,k k b, ACP,2~1 ,000 whereL L C,ACP ands s C,ACP are the length and diameter of a cylindrical arm of ACP, h 1,eq is the equilibrium angle between two arms of ACP, and h 2,eq is the angle formed by an arm of ACP and the axis of the actin filament to which the ACP is bound.k k s,ACP ,k k b, ACP,1 , andk k b, ACP,2 are stiffness constants related toL L C,ACP , h 1,eq , and h 2,eq , respectively ( Figure 11A). The value of h 2,eq is the same as in the previous model, and h 1,eq of ACP C was adjusted to maintain an equilibrium minimal distance of 70 nm between two crosslinked filaments.k k s,ACP was set to be one fortieth that of an actin filament due to computational efficiency.k k b,ACP, 1 andk k b, ACP,2 of ACP were estimated to be similar to the bending stiffness of an actin filament. Other parameters were also modulated according to the altered scale, as listed in Table 1. In our previous model by which the network was generated [35], a crosslink was allowed only if the torsional angle between two filaments was close to the equilibrium value (0 for ACP B and p/2 for ACP C ), and a finite stiffness was assigned to the torsional angle between crosslinked filaments. However, because the other two bending forces can preclude free torsional rotation, the torsional force is neglected here for simplicity.

Bulk rheology -simulation
The concept of a strain-controlled bulk rheometer used in experiments was adopted in order to measure the viscoelastic moduli of the generated networks. First, all actin filaments were severed at the upper and lower boundaries, and the periodic boundary condition on those surfaces was deactivated. Cylindrical actin segments within 70 nm from the bottom surface were fixed, whereas those within the top 70 nm were forced to move following an imposed strain. For the application of prestrain, the top boundary was translated at a constant strain rate, _ c c x , up to the desired strain. To measure differential viscoelastic moduli as in [14], a small sinusoidal strain (5%) was superposed on top of the finite prestrain. The sum of forces on the ends of filaments attached to the top boundary was calculated, and divided by the area of the top surface to determine stress. Only the force component parallel to the surface was considered. In addition, due to the small dimension of the system, the time scale for water diffusion through the computational domain is of order 10 ms. Since this is smaller than the smallest period of oscillatory strain, we assumed that the imposed shear strain immediately induces a linear velocity profile within the fluid. Consequently, after calculating the stress due to filament forces, we added a shear stress expressed as: where t zx is the shear stress exerted in the x-direction on a plane perpendicular to the z-direction (pointing from the bottom to the top face), and c x is the shear strain applied in the x-direction. The induced velocity of medium affects the movement of elements via the v ' term of Equation 2.
Finally, dividing the measured stress (t zx ) by the differential strain (c x ), viscoelastic moduli, G9 and G0, can be evaluated: where w is phase delay between strain and stress, and t zx j j and c x j j are the amplitude of the differential stress and differential strain, respectively. Note that w is zero for a perfectly elastic material and equals p=2 for a perfectly viscous one.

Bulk rheology -experiment
We also measured the mechanical properties of in vitro F-actin networks with a rheometer (AR-G2, TA Instruments) using a 40 mm parallel plate geometry. Lyophilized actin monomer from rabbit skeletal muscle was purchased from (Denver, CO). To minimize artifacts caused by sample preparations [52], the actin was stored at high concentration (10 mg/ml) at 280uC and thawed rapidly at 37uC before each experiment. Recombinant filamin A was purified from Sf9 cell lysates, and recombinant human gelsolin was produced in Escherichia coli. Solutions of gelsolin, filamin, actin polymerization buffer, and actin were gently mixed. The solutions were then loaded within 10 s into a rheometer to form a crosslinked F-actin network.
After 2 hr of polymerization at room temperature, frequencydependent shear moduli, G9 and G0, were measured in the range of 0.1-10 Hz. To obtain mechanical properties in the linear elastic regime, strain was maintained below 2%.

Segment-tracking rheology
Many experiments have used microbead rheology to probe viscoelastic moduli of actin gels based on the concept that the thermal motion of the bead is reflective of the gel's viscoelastic properties. Here, we used a variation of this approach, and tracked thermal fluctuations of individual actin segments. First, the domain was divided into a cubic lattice comprised of N MSD cells (N MSD = 512) of equal volume. Then, one cylindrical element was randomly selected per cell, and MSDs of the center of these elements were recorded over time. Using a wellknown approximate method [36], G9 and G0 were calculated from the MSDs. Considering that this method was initially designed for a spherical bead, it was appropriately modified for cylindrical elements: where r b is the effective radius of an actin segment, which satisfies f C,A,\~6 pgr b . C 1zx(v) ½ is the gamma function, wherex(v) is the power law exponent describing the logarithmic

Computational domain size and attainable time range
In spite of the coarse-graining using cylindrical segments, the length and time scales that our model can attain are still much smaller than those of usual experiments. For example, L f in this computational model is a few microns at most, and the width of the 3-D domain is less than 5 mm. In such a small domain, it is difficult to investigate the effects of a wide range of L f or C A on viscoelastic moduli since L f cannot be longer than the width of the domain to avoid artifacts associated with self-repulsion, and because an increase in C A was achieved by a decrease in domain size with a constant number of molecules. If we increase the number of molecules, simulation time significantly increases. On the other hand, many in vitro studies have used quite long actin filaments (,20 mm) with relatively large systems of size ,O(1 mm). It takes about 16 days to reach 1 s in typical conditions using an Intel Xeon 2.66GHz CPU, but experimental results span up to 100-1000 s. Thus, an exhaustive comparison between computational and experimental results was not possible.