PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Comput Methods Biomech Biomed Engin. Author manuscript; available in PMC 2010 May 21.
Published in final edited form as:
PMCID: PMC2874252
NIHMSID: NIHMS79295

A microstructurally informed model for the mechanical response of three-dimensional actin networks

R.Y. KWON,*,**## A.J. LEW,## and C.R. JACOBS**##

Abstract

We propose a class of microstructurally informed models for the linear elastic mechanical behavior of cross-linked polymer networks such as the actin cytoskeleton. Salient features of the models include the possibility to represent anisotropic mechanical behavior resulting from anisotropic filament distributions, and a power-law scaling of the mechanical properties with the filament density. Mechanical models within the class are parameterized by seven different constants. We demonstrate a procedure for determining these constants using finite element models of three-dimensional actin networks. Actin filaments and cross-links were modeled as elastic rods, and the networks were constructed at physiological volume fractions and at the scale of an image voxel. We show the performance of the model in estimating the mechanical behavior of the networks over a wide range of filament densities and degrees of anisotropy.

Keywords: Homogenization, Constitutive model, Finite element model, Actin cytoskeleton, Cell mechanics

1 Introduction

Numerous experiments have shown mechanical loading to be an important factor in the development and/or maintenance of a wide variety of tissues such as muscle, cartilage, tendon, and bone, and organs such as the heart and lung. The deformations of the cells within these tissues and organs are dictated by their mechanical behavior under loading. Thus, it comes as no surprise that cellular mechanical behavior has been implicated as an important factor in the pathology of many diseases such as osteoporosis, osteoarthritis, cancer, heart failure, and several pulmonary disorders [1].

Our understanding of the mechanical regulation of the pathologic processes involved in these diseases would be greatly enhanced if it were possible to predict the mechanical behavior of a particular cell from microscopically obtained observations. A critical component governing the mechanical behavior of adherent cells is the actin cytoskeleton, a three-dimensional network of cross-linked actin filaments (figure 1). The microstructure of the actin cytoskeleton is highly dynamic and can change dramatically in response to mechanical loading. A growing body of evidence suggests that the ability of cells to convert mechanical signals into biochemical signals depends on the actin cytoskeletal microstructure [2, 3]. Microstructurally-based models of the actin cytoskeleton would be ideal for investigating the mechanical implications of actin microstructural organisation, since representative cytoskeletal networks observed in vitro could be utilized. Important microstructural features, such as spatial and angular heterogeneity, could be directly accounted for, allowing investigation of underlying mechanical ‘principles’ that may be governing cytoskeletal microarchitecture.

Figure 1
Fluorescent image of the actin cytoskeleton of a MC3T3-E1 osteoblastic (bone) cell. Microstructurally, the actin cytoskeleton is highly heterogeneous in both the number and orientation of filaments at each point. The inset contains a brightfield image ...

Homogenized models for the actin cytoskeleton have yet to be obtained. A homogenized model in this context is a constitutive model for a continuum capable, in some appropriate sense, of approximating the mechanical behavior of the network. The advantage of a homogenized model is that the details of the network microstructure are used to generate the model, but are not needed thereafter. Although some advanced results are available for structured networks [4], rigorous mathematical results on general network homogenization problems remain elusive. Recent and classical theoretical investigations of ‘stiff’ or ‘semi-flexible’ polymer networks have yielded important insight into the mechanics of this class of networks and have generally identified geometric properties of random networks [5, 6, 7], different elastic regimes [8], scaling behaviors, and methods for explicit calculation of the macroscopic network elastic moduli from microscopic properties [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. All of these results but [24] have been obtained for two-dimensional, isotropic networks. However, anisotropy is highly relevant for actin networks, which forms aligned bundles of actin within the cytoskeleton (both in cultured cells and in vivo [27]) as well as cross-linked gels [28, 29]. The direct consequence of anisotropy is that, as opposed to previous investigations that needed to solely focus on the tensile and/or shear modulus that characterize isotropic materials, a suitable homogenized model will require the calculation of the entire elastic tensor (21 independent components) to fully specify the mechanical behavior of the network.

Incomplete knowledge of microscopic network properties is a unique challenge which makes constitutive modeling of many biopolymer networks difficult. In the case of the actin cytoskeleton, although the dimensions and material properties of individual actin filaments have been measured [30], many important microscopic properties of actin cytoskeletal networks have yet to be elucidated. For example, there are a wide variety of cytoskeletal cross-linking proteins in vivo whose mechanical behavior need to be characterized. In addition, although it is well accepted that the cytoskeletal network can be subject to a prestress, the degree to which each filament is prestressed is not known.

We propose here a novel class of models for the homogenized linear elastic response of cross-linked polymer networks such as the actin cytoskeleton. The proposed models can be constructed based on the filament angular distribution and spatial density. We deliberately avoided making specific assumptions on whether the elasticity of the network is the result of entropic [11] or enthalpic [8, 15] contributions, the nature of the cross-links between filaments in the network, or the existence of a prestress in the network, since these are still largely under discussion. Instead, our goal was to formulate a class of models that account for some features of the microstructure of the network, and that, through a suitable validation procedure, could be tailored to represent its homogenized linear elastic response under any or many of these conditions.

The 21 elastic moduli of the model are determined by postulating an ansatz or functional form inspired from the exact expression for affinely deformed networks with anisotropic filament distributions, see e.g., [31, 32]. The model accounts for the possibly different-than-linear exponents observed in the power law dependence of the elastic moduli with filament density (see, e.g., [8, 15]), and the effect of cross-links in the Poisson ratio. There are only 7 independent parameters, which need to be calibrated from a relatively small number of simulations of fully resolved and explicitly represented networks. We showcase the performance of the model by predicting the mechanical response of finite element models of three-dimensional, anisotropic networks of elastic rods with semi-flexible cross-links. The calibrated model shows good performance over a wide range of angular distributions and spatial densities away from the vicinity of the point of calibration. The particular type of networks chosen for this example was motivated by two-dimensional analogs that have been previously adopted as possible descriptions for the actin cytoskeleton [13, 15]. We expect, however, that the homogenized class of models proposed herein will also be useful to express the effective behavior of more general network types, resulting from a future enhanced understanding of key features of the actin cytoskeleton.

Throughout, vectors are denoted by boldface lowercase latin characters, second-order tensors by boldface lowercase greek characters, and fourth-order tensors by boldface uppercase latin characters. All tensor components are referred to an orthonormal basis. All nonboldface characters are considered scalar quantities. When indicial notation is used, an index appearing twice in a term indicates sum over it in the range 1 to 3.

2 Model Formulation

We propose and detail here a class of models to approximate the homogenized linear elastic response of cross-linked polymer networks such as the actin cytoskeleton. We idealize network filaments as cylindrical rods of cross-sectional area A with extensional, bending, and torsional stiffness. The configuration of each filament can be characterized by its midpoint position and direction in space when unstressed. We consider three-dimensional, infinite networks that are periodic with period L in three orthogonal directions. The unit cell of the network is then a cubic box of side L. A cross-link between two filaments may be formed whenever the distance between the two is within a specified distance. We denote with ρ the volume fraction of filaments in the network, i.e., the quotient between the total volume occupied by all filaments and the volume of the unit cell L3. Additionally, let ω(n) be the angular probability density of the volume fraction, which indicates the angular distribution of the volume fraction ρ. The value of ω(n)ρ dS gives the volume fraction of filaments with orientation in a neighborhood dS of the unit sphere surrounding n. It satisfies that its integral over the unit sphere S2 is equal to one.

Motivated by homogenization results for linear periodic composite materials (see e.g., [33]), we expect the elastic energy density of the homogenized model to be well approximated by the lowest elastic energy attainable by the network under suitably imposed boundary conditions on the unit cell. For this study, the points of intersection between the network filaments and the unit cell faces are constrained to follow an affine deformation. A class of less stringent boundary conditions requiring only periodic but not necessarily affine displacements on the boundary of the unit cell are also of interest, but will not be considered here. If x denotes the position vector of one such point and ε is a second-order tensor, then an affine deformation Fε(x) maps x to x+ε.x. Since we are first interested in the linearized elastic behavior of initially stress-free networks, only affine deformations in which ε is small and symmetric need to be considered. In this case, the elastic energy density of the homogenized models has the form w(ε)=1/2Cijklεijεkl, where Cijkl are the Cartesian components of the homogenized or effective elastic moduli. The elastic moduli should satisfy that Cijkl=Cklij=Cjikl (see, e.g., [34]), also known as major and minor symmetries, respectively.

Even though homogenization of the mechanical properties for general networks is still an open problem, there are some cases in which the elastic moduli are easily obtained. Consider for example the case of an initially unstressed network in which every filament has only extensional stiffness K// (no bending or torsion). All cross-links in the network are assumed to be able to freely rotate with no energetic cost. It is easily seen that upon constraining the filament ends to follow an affine deformation Fε, an equilibrium configuration of the network is obtained by mapping all cross-links with Fε as well (in the linear elastic regime of concern here, such network deformation has the lowest energy among all those that use the same boundary conditions; it may not be unique though). Herein, we call this class of networks ‘affine networks’. As we shall detail later, the homogenized elastic moduli in this case are given by

Cijkl=KS2ρω(n)ninjnknldS.
(1)

As nicely discussed in [12], this result will not strictly hold when the filament ends are not constrained to follow an affine deformation. Previous numerical studies have also shown that the linear scaling of the moduli with ρ is not longer valid for more complex networks, such as when bending stiffness is accounted for. Instead, power laws of the form (ρ-ρref)α for some exponent α≥1 and reference density ρref have been found to give a good approximation to the dependence of the shear or Young modulus with the density for isotropic networks, over spans of several orders of magnitude [10, 15]. The value of ρref is typically related to the percolation threshold, e.g., [23, 24].

We are therefore motivated to posit the following form (ansatz) for the homogenized elastic moduli

Cijkl=S2(ρρref)αω(n)Kijkl(n)dS.
(2)

Here, Kijkl(n) is a fourth-order tensor valued function of the direction vector n, while ρref and α are scalars, all of them to be determined, and we consider only volume fractions ρ>ρref. In the particular case of affine networks in equation (1), we have ρref=0, α=1 and

Kijkl(n)=Kninjnknl.
(3)

The homogenized elastic moduli in equation (2) should have major and minor symmetries for any choice of angular density ω(n). It then follows that K(n) should also obey these same symmetries, and hence Kijkl=Kklij=Kjikl. This leaves at most 21 independent components of K(n) to be determined, all of them functions over the unit sphere. However, we show next that these components cannot be arbitrary functions of n, but that these are fully determined by five independent constants and a specific functional dependence on n. The elegant derivation below significantly reduces the complexity of the model and sheds light into the physical interpretation of K(n).

A physically reasonable requirement on the elastic moduli Cijkl is that if the entire network is rigidly rotated by a rotation Q then the elastic moduli should rotate accordingly, also known as the material frame indifference requirement in continuum mechanics, see e.g., [35]. More precisely, after rotating the network by Q, a new network is obtained with an angular probability density ω(n)= ω(Q-1n) and elastic moduli

Cijkl=S2(ρρref)αω(n)Kijkl(n)dS=S2(ρρref)αω(Q1n)Kijkl(n)dS=S2(ρρref)αω(n)Kijkl(Qn)dS.
(4)

Correspondingly, by material frame indifference the elastic moduli of the rotated network should be given by CIJKL=QIiQJjQKkQLlCijkl, or equivalently,

CIJKL=QIiQJjQKkQLlS2(ρρref)αω(n)Kijkl(n)dS.
(5)

Equating Eqs. 4 and 5, and utilizing the fact that the resulting identity should be valid for any filament angular probability density ω(n), we obtain

KIJKL(Qn)=QIiQJjQKkQLlKijkl(n),
(6)

which precisely defines the dependence of K on n. Notice that determining K in a single direction n0 allows its calculation in any other direction n through an appropriately chosen rotation Q. In fact, there are an infinite number of rotations that map n0 to n. If Q denotes one such rotation, then any of the others is obtained as QR, where R is any rotation about n0. This implies that for any such R, with Q equal to the identity,

KIJKL(n0)=RIiRJjRKkRLlKijkl(n0),
(7)

or alternatively, the tensor K(n0) must be transversely isotropic about n0 (see e.g., [36]). There are consequently only five independent parameters in K(n0). These values are easily expressed in the Voigt basis [36], the 6×6 matrix representation of K(n0). Assuming that the axis of symmetry n0 is parallel to e3, the unit vector parallel to the third coordinate axis, we get

K(e3)=[K1111(e3)K1122(e3)K1133(e3)000K1122(e3)K1111(e3)K1133(e3)000K1133(e3)K1133(e3)K3333(e3)00000012(K1111(e3)K1122(e3))000000K1313(e3)000000K1313(e3)].
(8)

Here K1111(e3), K1122(e3), K1133(e3), K3333(e3), and K1313(e3) are the five independent constants to be determined. Once these constants are known, the entire tensor valued function over the unit sphere K(n) can be calculated using equation (6). Together with ρref and α, this makes a total of seven constants in equation (2) needed to completely specify the model.

In order to gain insight into the physical origin of the five constants which make up K(n), we revisit the case of affine networks and compute the homogenized moduli in equation (1) and K in equation (3) for this class of networks next. The elastic energy of a filament of per unit filament volume is

Wfilament(λ)=12Kλ2,
(9)

where λ is the axial strain. In the linear elastic regime and under an affine deformation Fε the value of λ for a filament oriented in a direction n can be computed as nrεrsns. The elastic energy density of the homogenized model is the result of adding the elastic energies contributed by each filament in the network divided by the volume of the unit cell, i.e.,

w(ε)=S2ρω(n)Wfilament(niεijnj)dS,
(10)

and the homogenized elastic moduli then follow as

Cijkl=2w(ε)εijεkl=S2ρω(n)KijklninjnknldS.
(11)

For this case, K3333(e3)=K//, which explicitly shows the connection between K3333(e3) and the extensional stiffness of the filaments. In fact, K3333(e3) is the only one out of the five different constants in K(e3) that is different from zero. Recall that this result is a direct consequence of the fact that the cross-links are mapped affinely. Thus, the nonzero values of the other coefficients could serve as indicators of the degree to which the cross-links are mapped affinely when K(e3) is computed for more general networks of filaments with extensional stiffness only.

These last observations provide the basis for an interesting remark about the model. The material frame indifference argument showed that the proposed ansatz is obtained by arranging a transversely isotropic material with elastic moduli K in different directions. For affine networks, the cross-links do not play a role in the deformation, i.e., they can be removed and the network deforms in the same way. This is reflected in the elastic moduli K, since only the extensional stiffness in the direction of the symmetry axis, the direction of the filament, is different from zero. However, the model can admit more general moduli K when the other four constants are also different from zero. The presence of nonaffinely mapped cross-links takes advantage of this generality. In this case, in addition to the extensional stiffness provided by filaments in the direction of the axis of symmetry, K has a contribution of the bulk network. For example, by pulling in the direction of the axis of symmetry, a contraction in the transverse directions may be obtained, i.e., a nonzero Poisson ratio (figure 2). The model is then able to represent the fact that each filament transmits forces along its own direction and, through its connections to the rest of the network, in transversal directions as well. This observation suggests that the proposed model has the best chances to be accurate when each filament encounters a similar type of cross-linking with the rest of the network regardless of its orientation, for example, a similar number of cross-links. If this were not the case, the transverse mechanical properties in K would be orientation dependent, in direct contradiction of equation (6).

Figure 2
Schematic depicting a simple, non-periodic network in which pulling the vertical filament in the axial direction results in a contraction in the transverse direction. In this case, the filament transmits forces along its own direction and, through its ...

Simple examination reveals some of the shortcomings of the proposed model. For example, consider an initially unstressed network formed by filaments whose elastic energy arises solely due to bending (no torsion or stretching) and whose cross-links are rigid, i.e., the angle between any two filaments at a cross-link cannot change. The filaments are arranged in a prismatic microstructure so that filament directions are parallel to three orthonormal vectors e1, e2, and e3, while the distance between consecutive cross-links in each direction, are L1, L2 and L3, respectively. As we shall show next, the homogenized linear elastic moduli for this network depend on quadratic moments of ω(n) (i.e., products of the type ω(ei)ω(ej)) which account for the number of intersections between filaments in directions ei and ej. Clearly, this quadratic dependence on ω(n) is not represented in the model, so we do not expect the ansatz to be predictive for networks in which most of the elastic energy is used to bend the filaments

For the sake of simplicity, we shall compute only the shear modulus G=C1212. Notice then that the periodicity of the microstructure enables us to consider a minimal unit cell with dimensions L1, L2 and L3, which substantially simplifies the computation. We are interested in the energy of the network when an affine deformation Fε, with ε=γ(e1[multiply sign in circle]e2+e2[multiply sign in circle]e1)/2 for some small constant γ>0, is used to map the positions of the cross-links in the unit cell, and hence all cross-links in the network. It is straightforward to check that by allowing the filaments to bend so that the resultant torque at each cross-link is identically zero, an equilibrium configuration of the network is obtained. For the boundary conditions under consideration, only the filaments parallel to e1 and e2 bend, and they do it within the e1-e2 plane (figure 3). The bending energy per unit length for each one of these filaments is given by K[perpendicular]κ2/2, where κ denotes the local curvature and K[perpendicular] the bending stiffness. The linear elastic energy per unit volume is then easily computed to be

w(γ)=1L1L2L324γ2K(L1+L2).
(12)

The volume fraction of the network is obtained as ρ=A(L1+L2+L3)/(L1L2L3), while the fraction of the total filament volume occupied by filaments parallel to ei is denoted with ωi=Li/(L1+L2+L3). It then follows that

w(γ)=ω1ω2ω3ω1+ω224γ2KA2
(13)

and that

G=2w(γ)γ2=ρ2ω1ω2ω3ω1+ω248KA2.
(14)

The last equation clearly shows that the stiffness for this type of networks depends on the product between values of the angular probability distribution in different directions. Additionally, the dependence of G on the angular probability distribution is highly nonlinear. These two facts are in direct contradiction to the proposed model in equation (2), which explicitly contains a linear dependence on ω(n). Notice, however, that the model is able to capture the nonlinear stiffening of the network with the volume fraction.

Figure 3
Three-dimensional network with prismatic microstructure (left) and representative unit cells when undergoing shear (right). If the cross-links are rigid, and the network undergoes shear, the filaments in the plane of shear bend. In contrast, if the cross-links ...

3 Numerical methods

Generally, ρref,α, and the five independent components of K(n) cannot be analytically solved for. We outline here a procedure to numerically determine these parameters. We demonstrate the validity of this procedure and the proposed model by performing numerical experiments with finite element models of a specific type of cross-linked networks. By using finite element models, we can solve for the homogenized elastic moduli of each network exactly, and thus have a ‘gold standard’ to which we can directly compare elastic moduli calculated using equation (2).

3.1 Finite element modeling and computation of elastic moduli

Finite element models of three-dimensional, periodic, cross-linked networks were constructed by placing filaments of length 350nm inside cubic domains of length 400nm (figure 4). These length scales were selected since they may be relevant for image-based finite element models of the actin cytoskeleton (i.e., finite element models meshed directly from three-dimensional image data). For example, a typical pixel/voxel size in fluorescent microscopy is 0.2~0.4μm, and actin filaments are approximately in length in vivo [37, 38, 39]. Filament centroids were randomly sampled from a uniform distribution in the box, while distributions with different degrees of anisotropy were used for the orientations. Without loss of generality, we idealized actin filaments as elastic rods with cylindrical cross sections. Actin filaments were modeled as Euler-Bernoulli beams with linear elastic extensional and torsional stiffness, with a diameter d=8nm and Young modulus E=1.8GPa [30]. Although there are a large number of different cross-linking proteins in vivo with a wide variety of lengths and microscopic properties (many of which have yet to be elucidated), for simplicity, we assume filaments closer than a rod diameter to each other were cross-linked by rods having identical properties as the actin filaments, with lengths equal to the distance between the two filaments that were cross-linked. Of direct consequence is that unlike two-dimensional networks (e.g.,, [8, 15]), the rod diameter played a critical role in the network connectivity, since the probability of a fixed number of filaments to intersect asymptotically vanishes with the rod diameter [40]. We imposed periodic boundary conditions on the walls of the unit cell by constraining the points of intersection between filaments and the cubic domain boundary to follow a prescribed affine deformation, while the torques at the same points were constrained to be continuous across the two parts of the filament lying in neighboring unit cells. The total strain energy density for a given deformation ε imposed on the boundary was computed using the finite element method (Abaqus, Providence, RI). Once the strain energy densities were found, the homogenized network elastic moduli were calculated as

Cijkl=2wεijεkl(0)[w(Δδij+Δδkl)w(Δδij)][w(Δδkl)w(0)]Δ2,
(15)

where Δ>0 is a small number. By using equation (15), calculating all 21 independent elastic moduli of the fully anisotropic elastic tensor required 27 analyses per network.

Figure 4
Typical three dimensional, cross-linked actin network used in this study. The actin filaments and cross-links are modeled as elastic rods. Shown is the Von Mises stress distribution calculated using the finite element method when the network undergoes ...

3.2 Numerical determination of ρref, α, and the five independent components of K(n)

Seventy networks with varying volume fractions were generated, and their elastic moduli numerically computed. Physiological volume fractions of filamentous actin within different types of cells and cytoplasmic regions are still under discussion. However, the average volume fraction within bovine aortic endothelial cells has been reported to be on the order of 1% [9]. Since we expect the volume fraction to be higher than the average within actin-rich regions such as the cortex or actin bundles, we constructed networks with volume fractions ranging from ρ=0-0.10. Note that “dangling” filament ends (i.e., filament segments attached to a cross-link at one end but to nothing at the other) were accounted for when quantifying volume fractions (these segments do not contribute to the elastic energy of the network). We determined ρref by finding the minimum volume fraction at which the elastic tensor was non-zero among all tested networks. We found α by plotting C3333 as a function of volume fraction, and fitting this curve to a power law scaling of the form (ρ-ρref)α (figure 5). Next, we computed the angular distribution ω(n) for a single network (ρ=0.05), and determined the five independent components of K(n) using a direct search optimization to perform a least square fit, i.e., such that they minimized the norm of the difference between the network’s elastic tensor calculated using finite elements, and the network’s elastic tensor calculated using equation (2). Notably, performing this optimization for one network or many of them simultaneously did not appreciably change the value of K(n).

Figure 5
Log-log plot of elastic moduli C3333 and C1313 versus volume fraction ρ. The elastic moduli scale as ~ρ1.6 (black lines). Thus, a value of α=1.6 was used in equation (2).

3.3 Predicting the elastic moduli of networks of varying degrees of anisotropy and volume fraction

Once ρref, α, and K(n) were determined, we sought to test the accuracy of equation (2) in predicting the elastic moduli of networks with different volume fractions and varying degrees of anisotropy. To this end, we generated additional networks for which the filament orientations were randomly sampled from a family of angular probability distributions ωξ, 0≤ξ≤1, constructed as a normalized linear combination of two real spherical harmonics S00 and S10, i.e.,

ωξ(n)=(1ξ)S00(n)+ξS10(n)S2[(1ξ)S00(n)+ξS10(n)]dS.
(16)

Here S00 denotes the isotropic distribution function over the unit sphere, while S10 stands for a dipole distribution aligned in the e3 direction (if ϕ is the angle between e3 and n, then S10(ϕ)=(√3cos(ϕ))/(2√π)). Constructing networks with ξ=0 resulted in networks with filaments isotropically aligned, while constructing networks with ξ=1 resulted in networks with filaments primarily aligned in the e3 direction. In general, constructing networks with ξ=1 resulted in networks with an approximately threefold increase in stiffness in C3333 compared to C1111 and C2222. Note that the number of cross-links per network did not appreciably change with the degree of anisotropy tested here. We generated 170 networks with different volume fractions (0≤ξ≤0.10) and degrees of anisotropy (0≤ξ≤1), for which we computed their angular distributions ω(n). Using equation (2), and the values of ρref,α, and K(n) found previously, a predicted value for the elastic moduli of each network was computed. These values were then compared against the exact elastic moduli of each network computed with the finite element model. For each network the error between the predicted and exact values for the elastic moduli, Cmodel and Cexact, respectively, was found as ||Cmodel-Cexact||2/||Cexact||2. Note that in general, we did not observe larger error in any one component of the elastic tensor for isotropic and anisotropic networks. Finally, in order to visualize the fourth-order elastic tensor C, we generated polar plots of the quantity Cijklninjnknl as a function of direction n. Physically, Cijklninjnknl gives the normal component of the traction in direction n generated by a uniaxial unit strain in the same direction.

4 Results

The numerically determined values of ρref and α were 1.6 and 0.01, respectively. The five independent components of K(n), in the form of equation (8), were K1111(e3)=5.4MPa, K1122(e3)=9.7MPa, K1133(e3)=13.0MPa, K3333(e3)=3.8GPa and K1313(e3)=37.0MPa (figure 6). Using these values, equation (2) predicted the elastic tensors of networks generated with constant volume fraction (ρ=0.05) and varying degrees of anisotropy (0≤ξ≤1) to within 14.5±0.6% (mean±SE, n=100 networks1). As figures figures77 and and88 show, there seems to be only a minimal dependence of the error on anisotropy. Next, we used equation (2) to compute predicted values for the elastic moduli of both isotropic (ξ=0) and anisotropic (ξ=1) networks generated over a range of different volume fractions (0≤ρ≤0.10). The results are shown in figure 9. For both isotropic and anisotropic networks, the error was found to be relatively constant with respect to the volume fraction, except for volume fraction values approaching ρref, for which the errors increase rapidly. For example, for networks with volume fractions far from ρref, the error was relatively small (ξ=0: 12.6±0.9%, n=48; ξ=1: 12.1±0.9%, n=48), while for networks with volume fractions greater than but near ρref(ρ=ρref-0.03), the error was substantially larger (ξ=0: 52.8±5.8%, n=13; ξ=1: 57.1±6.7%, n=14). For networks with volume fractions smaller than ρref the elastic tensor was identically zero, reflecting the fact that ρref is a good estimate of the percolation threshold for these networks.

Figure 6
Polar plots of the components of KIJ(n), the 6×6 matrix representation of K. The plot in the Ith row and Jth column is the polar plot of KIJ(n). Positive values are represented in gray, negative values are represented in black. The magnitude of ...
Figure 7
Error in the predicted elastic moduli as a function of the degree of anisotropy ξ. A minimal dependence of the error on anisotropy is observed.
Figure 8
Polar plots of Cijklninjnknl as a function of direction n using elastic moduli predicted by the model (Cmodel, left column), and generated using elastic moduli calculated using the finite element method (Cexact, right column). The top row represents an ...
Figure 9
Error in the predicted elastic moduli as a function of volume fraction ρ for networks with different degrees of anisotropy (ξ=0,1). For both isotropic and anisotropic networks, the error was constant far from ρref~0.01, ...

5 Discussion

We proposed here a novel class of models for the homogenized linear elastic response of cross-linked polymer networks. The model requires determination of only seven independent parameters: five constants to construct a fourth-order tensor valued function K(n) that incorporates different aspects governing network mechanical behavior (e.g., filament elasticity, cross-link type, etc.), and two constants, ρref and α, to account for nonlinear scaling of network elastic moduli with filament density.

We demonstrated the applicability and validity of these models by numerically determining ρref, α, and the five independent constants used to construct K(n) for a particular class of network. We obtained a value of 1.6 for α, indicating that the elastic moduli in our three-dimensional networks scaled with filament density as ~(ρ-ρref)1.6. This power law dependence is similar to the scaling exponent of 1.8 found in two-dimensional rod networks with stiff cross-links near the percolation threshold [15]. In addition, we found ρref to be 0.01, close to values obtained in several studies that suggest that for rods the critical volume fraction for percolation takes the form ρref~d/2l~0.01 [41, 42, 43]. Out of the five different values that determine K(n), we found that K3333(e3) was two or three orders of magnitude greater than the other components. This indicates that networks in the class considered herein predominantly absorb deformation by stretching of the filaments, instead of from bending or torsion. The nonzero value of the other coefficients indicates that either some degree of nonaffinity exists in the deformed position of the cross-links, or simply reflects the presence of the not freely rotating cross-links.

When we calculated the eigenvalues of K(e3) over all symmetric second-order tensors, we found them to be all positive. This implies that the homogenized linear elastic moduli will be positive definite over all symmetric second-order tensors for any volume fraction ρ>ρref and angular probability density ω(n), a fact that renders the deformations of any of the resulting homogenized models elastically stable.

Once we numerically determined K(n), ρref and α, we used equation (2) to compute the predicted values for the elastic tensors of networks over a range of volume fractions (0≤ρ≤0.10), and with varying degrees of anisotropy (0≤ξ≤1). We discuss here several trends in the error in the predicted homogenized elastic moduli with respect to volume fraction and anisotropy. First, although the parameters in equation (2) were numerically derived from isotropic networks, the model generally predicted the elastic moduli of isotropic and anisotropic networks equally well. However, we expect that the error may be substantial for networks with higher degrees of anisotropy that those tested in this study. For example, for a theoretical network in which every filament is aligned in the e3 direction, the model predicts finite axial stiffness in the transverse directions, which is evidently not correct. Second, we found that the error was large for networks with volume fractions near ρref, regardless of the degree of anisotropy. This is indicative of the fact the elastic moduli near ρref may scale differently with volume fractions at these very low values. Finally, we found that the error was relatively constant for volume fractions far from ρref. It would not be unexpected for the error to increase substantially for networks with volume fractions larger than those tested in this study. At the volume fractions studied here the mechanical response of the network is determined mainly by the stretching mechanical behavior of the filaments. However, at higher volume fractions other phenomena may become important, as some two-dimensional studies indicate [10, 13]. Alternatively, simple inspection of the previous example consisting in a cubic microstructure with filaments that only have bending stiffness suggests that at higher volume fractions the bending response may become dominant, since in this case the stiffness scales quadratically with the volume fraction. We have not yet been able to experimentally verify this regime due to the large computational cost associated with high density networks in three dimensions.

Overall, once calibrated the model predicted the elastic moduli for these networks over a range of volume fractions and degrees of anisotropy. In fact, the numerical results are surprisingly good, since the variation of the moduli with the anisotropy parameter was well captured, despite the fact that the model parameters were determined with isotropic networks only. More generally, we expect the proposed model, with perhaps some minor extensions, to be useful in different regimes and for other types of networks not considered herein as well, if calibrated accordingly. This, of course, has to be explicitly studied for each network type, e.g., different cross-linker microscopic properties, prestress in the network and length distribution of actin filaments, among others.

The distinguishing feature of this class of models is that they can account for filament angular distribution and volume fraction, two microstructural parameters that can be estimated from light microscopy data. For example, although individual filaments cannot be resolved in fluorescent images (a typical pixel/voxel size in fluorescent microscopy is 0.2~0.4μm, whereas actin filaments are approximately 8nm in diameter), one can approximate the number and orientations of filaments within each image voxel. Specifically, the volume fraction of filaments within each voxel can be estimated from the voxel intensity, while their orientations can be estimated from the local image texture [44, 45] to give an approximate angular distribution of filaments [44]. Thus, the approach presented here may serve to construct cell structural models that account for the anisotropic and heterogeneous distribution of actin filaments, something that has yet to be explored and validated against appropriate experiments. Despite the facts that the actin cytoskeleton is a highly dynamic structure capable of remodeling, and that cells can undergo large deformations, the linearized mechanical response studied herein may provide important insight for the interpretation of a wide class of microrheological experiments.

6 Acknowledgments

We would like to acknowledge grant AR45989 from the National Institutes of Health, the National Science Foundation Graduate Fellowship, and the Veterans Affairs Palo Alto Bone and Joint Center for funding.

Footnotes

1SE: standard error, n: number of samples

AMS Subject Classification: 74D05; 74K10; 74S05; 74E10

References

1. Ingber DE. Mechanobiology and diseases of mechanotransduction. Annals of Medicine. 2003;35:1–14. [PubMed]
2. Pavalko FM, Chen NX, Turner CH, Burr DB, Atkinson S, Hsieh YF, Qiu JY, Duncan RL. Fluid shear-induced mechanical signaling in MC3T3-E1 osteoblasts requires cytoskeleton-integrin interactions. American Journal of Physiology-Cell Physiology. 1998;44:C1591–C1601. [PubMed]
3. Chen NX, Ryder KD, Pavalko FM, Turner CH, Burr DB, Qiu JY, Duncan RL. Ca2+ regulates fluid shear-induced cytoskeletal reorganization and gene expression in osteoblasts. American Journal of Physiology-Cell Physiology. 2000;278:C989–C997. [PubMed]
4. Cioranescu D, Paulin J, Telega JJ. Homogenization of reticulated structures. Applied Mechanics Reviews. 2001;54:B64–B65.
5. Dodson CTJ. Spatial variability and the theory of sampling in random fibrous networks. Journal of the Royal Statistical Society B. 1971;33:88–94.
6. Miles RE. Random polygons determined by random lines in a plane. Proceedings of the National Academy of Sciences of the United States of America. 1964;52:901–907. [PubMed]
7. Goudsmit S. Random distribution of lines in a plane. Reviews of Modern Physics. 1945;17:321–322.
8. Head DA, Levine AJ, MacKintosh FC. Distinct regimes of elastic response and deformation modes of cross-linked cytoskeletal and semiflexible polymer networks. Physical Review E. 2003;68:061907. [PubMed]
9. Satcher RL, Dewey CF. Theoretical estimates of mechanical properties of the endothelial cell cytoskeleton. Biophysical Journal. 1996;71:109–118. [PubMed]
10. Astrom JA, Makinen JP, Alava MJ, Timonen J. Elasticity of Poissonian fiber networks. Physical Review Letters E. 2000;61:5550–5556. [PubMed]
11. Mackintosh FC, Kas J, Janmey PA. Elasticity of semiflexible biopolymer networks. Physical Review Letters. 1995;75:4425–4428. [PubMed]
12. Heussinger C, Frey E. Floppy modes and nonaffine deformations in random fiber networks. Physical review letters. 2006;97:105501. [PubMed]
13. Head DA, Levine AJ, MacKintosh EC. Deformation of cross-linked semiflexible polymer networks. Physical Review Letters. 2003;91:061907. [PubMed]
14. Heussinger C, Frey E. Stiff polymers, foams, and fiber networks. Physical Review Letters. 2006;96:017802. [PubMed]
15. Wilhelm J, Frey E. Elasticity of stiff polymer networks. Physical Review Letters. 2003;91:108103. [PubMed]
16. Shin JH, Gardel ML, Mahadevan L, Matsudaira P, Weitz DA. Relating microstructure to rheology of a bundled and cross-linked f-actin network in vitro. Proceedings of the National Academy of Sciences. 2004;101:9636–9641. [PubMed]
17. Sastry AM, Cheng X, Wang CW. Mechanics of stochastic fibrous networks. Journal of Thermoplastic Composite Materials. 1998;11:288–296.
18. Onck PR, Koeman T, van Dillen T, van der Giessen E. Alternative Explanation of Stiffening in Cross-Linked Semiflexible Networks. Physical Review Letters. 2005;95:178102. [PubMed]
19. Heussinger C, Frey E. The role of architecture in the elastic response of semiflexible polymer and fiber networks. Physical Review E. 2007;75:011917. [PubMed]
20. Ostoja-Starzewski M, Stahl DC. Random fiber networks and special elastic orthotropy of paper. Journal of Elasticity. 2000;60:131–149.
21. Tremblay RR, Day AR, Tremblay AMS. Splay rigidity in the diluted central-force elastic network. Physical Review Letters. 1986;56:1425. [PubMed]
22. DiDonna BA, Lubensky TC. Nonaffine correlations in random elastic media. Physical Review E. 2005;72:066619. [PubMed]
23. Feng S, Sen PN. Percolation on Elastic Networks: New Exponent and Threshold. Physical Review Letters. 1984;52:216–219.
24. Arbabi S, Sahimi M. Elastic properties of three-dimensional percolation networks with stretching and bond-bending forces. Physical Review B. 1988;38:7173–7176. [PubMed]
25. Chandran PLB, Victor H. Affine versus non-affine fibril kinematics in collagen networks: theoretical studies of network behavior. Journal of Biomechanical Engineering-Transactions of the Asme. 2006;128:259–270. [PubMed]
26. Head DA, Levine AJ, MacKintosh FC. Mechanical response of semiflexible networks to localized perturbations. Physical Review E. 2005;72:061914. [PubMed]
27. Wong AJ, Pollard TD, Herman IM. Actin filament stress fibers in vascular endothelial cells in vivo. Science. 1983;219:867–869. [PubMed]
28. Gardel ML, Shin JH, MacKintosh FC, Mahadevan L, Matsudaira P, Weitz DA. Elastic behavior of cross-linked and bundled actin networks. Science. 2004;304:1301–1305. [PubMed]
29. Wagner B, Tharmann R, Haase I, Fischer M, Bausch AR. Cytoskeletal polymer networks: The molecular structure of cross-linkers determines macroscopic properties. Proceedings of the National Academy of Sciences. 2006;103:13974–13978. [PubMed]
30. Gittes F, Mickey B, Nettleton J, Howard J. Flexural rigidity of microtubules and actin filaments measured from thermal fluctuations in shape. Journal of Cell Biology. 1993;120:923–934. [PMC free article] [PubMed]
31. Pinsky PM, van der Heide D, Chernyak D. Computational modeling of mechanical anisotropy in the cornea and sclera. Journal of Cataract and Refractive Surgery. 2005;31:136–145. [PubMed]
32. Baesu E, Zheng M, Beljic D. A continuum model for two-dimensional fiber networks. Mechanical Properties of Bioinspired and Biological Materials. 2004;844:235–241.
33. Cioranescu D, Donato P. An introduction to homogenization. Oxford University Press; Oxford: 1999.
34. Holzapfel GA. Nonlinear solid mechanics: A continuum approach for engineering. Wiley; Chichester: 2000.
35. Malvern LE. Introduction to the mechanics of a continuous medium. Prentice-Hall Englewood Cliffs; New Jersey: 1969.
36. Lai WM, Rubin D, Krempl E. Introduction to continuum mechanics. Pergamon Press; New York: 1993.
37. Podolski JL, Steck TL. Length distribution of f-actin in Dictyostelium-iscoideum. Journal of Biological Chemistry. 1990;265:1312–1318. [PubMed]
38. Hartwig JH, Shevlin P. The architecture of actin filaments and the ultrastructural location of actin-binding protein in the periphery of lung macrophages. Journal of Cell Biology. 1986;103:1007–1020. [PMC free article] [PubMed]
39. Medalia O, Weber I, Frangakis AS, Nicastro D, Gerisch G, Baumeister W. Macromolecular architecture in eukaryotic cells visualized by cryoelectron tomography. Science. 2002;298:1209–1213. [PubMed]
40. Dodson CTJ, Sampson WW. Spatial statistics of stochastic fiber networks. Journal of Statistical Physics. 1999;96:447–458.
41. Chaterjee AP. A model for the elastic moduli of three-dimensional fiber networks and nanocomposites. Journal of Applied Physics. 2006;100:054302.
42. Foygel M, Morris RD, Anez D, French S, Sobolev VL. Theoretical and computational studies of carbon nanotube composites and suspensions: Electrical and thermal conductivity. Physical Review B. 2005;71:104201.
43. Balberg I, Anderson CH, Alexander S, Wagner N. Excluded volume and its relation to the onset of percolation. Physical Review B. 1984;30:3933–3943.
44. Verkhovsky AB, Chaga OY, Schaub S, Svitkina TM, Meister J-J, Borisy GG. Orientational order of the lamellipodial actin network as demonstrated in living motile cells. Molecular Biology of the Cell. 2003;14:4667–4675. [PMC free article] [PubMed]
45. Lichtenstein N, Geiger B, Kam Z. Quantitative analysis of cytoskeletal organization by digital fluorescent microscopy. Cytometry Part A. 2003;54A:8–18. [PubMed]