P H Y S I C A L R E V I E W L E T T E R S week ending5 SEPTEMBER 2003VOLUME 91, NUMBER 10

Elasticity of Stiff Polymer Networks

Jan Wilhelm and Erwin FreyHahn-Meitner-Institut, Abteilung Theorie, Glienicker Strasse 100, D-14109 Berlin, Germany

Fachbereich Physik, Freie Universitat Berlin, Arnimallee 14, D-14195 Berlin, Germany(Received 31 March 2003; published 5 September 2003)

108103-1

We study the elasticity of a two-dimensional random network of rigid rods (‘‘Mikado model’’). Theessential features incorporated into the model are the anisotropic elasticity of the rods and the randomgeometry of the network. We show that there are three distinct scaling regimes, characterized by twodistinct length scales on the elastic backbone. In addition to a critical rigidity percolation region and ahomogeneously elastic regime we find a novel intermediate scaling regime, where the elasticity isdominated by bending deformations.

DOI: 10.1103/PhysRevLett.91.108103 PACS numbers: 87.16.Ka, 62.20.Dc, 82.35.Pq

depend on the number of cross-links and the density of hinges’’). Similar to thermally fluctuating semiflexible

The elasticity of cells is governed by the cytoskeleton,a partially cross-linked network of relatively stiff fila-ments forming a several 100 nm thick shell called theactin cortex [1]. While the statistical properties of singlecytoskeletal filaments are by now relatively well under-stood [2,3], theoretical concepts for the elasticity of stiffpolymer networks are still evolving. One major openquestion is how stresses and strains are transmitted insuch networks. In synthetic gels formed by rather flexiblechain molecules the response to macroscopic externalforces is — on the level of single filaments — isotropicand entropic in origin. It is generally believed that macro-scopic stresses are transmitted in such a way that localdeformations within the network stay affine, i.e., that theend-to-end distance of individual filaments follows themacroscopic shear deformation [4]. In contrast, the build-ing blocks of the actin cortex are semiflexible polymers,whose hallmark is an extremely long persistence length‘p comparable to the total contour length ‘. As a con-sequence, the response of such stiff polymers to externalforces shows a pronounced anisotropy [5]. Consider a semi-flexible polymer with one end clamped at a fixed orienta-tion. When forces are applied at the other end transverseto the tangent vector at the clamped end, the responsemay be characterized by a transverse spring coefficientk?�‘� � 3�=‘3 proportional to the bending modulus �.Whereas this response is of purely mechanical origin, thelinear response to longitudinal forces is due to the pres-ence of thermally excited undulations which make theaverage end-to-end distance of the polymer shorter thanits contour length. The corresponding effective springcoefficient kk�‘� � 6�2=�kBT‘4� is proportional to �2=Tindicating the breakdown of linear response for very stifffilaments [2]. In a typical network one expects the dis-tance between cross-links ‘c to be much smaller than thepersistence length and filament length. Hence we havekk�‘c�=k?�‘c� � 2‘p=‘c � 1, i.e., the elastic response ofthe filaments is indeed highly anisotropic.

The anisotropic elastic properties of individual fila-ments suggests that the macroscopic elasticity of networks

0031-9007=03=91(10)=108103(4)$20.00

filaments, but also on the geometry and architecture ofthe network. For some very regular networks such as atriangular lattice the longitudinal spring coefficient kkdominates the macroscopic moduli [6] since the networkcannot be strained without a change of the end-to-enddistance of individual polymers. In other regular networkarchitectures, the softer bending modes would be domi-nant [7]. Naturally, this will lead to a very differentprediction for the elastic modulus of the network. It isnot at all obvious what type of network geometry (elon-gation dominated versus bending dominated) is relevantin less ideal structures with a significant amount of dis-order as found in typical cytoskeletal networks.

As a first step towards understanding the elasticity ofstiff polymer networks we consider a two-dimensionalmodel defined as follows (see Fig. 1). We generate therandom network by placing N linelike objects of equallength ‘ on a plane with area A � L2 such that bothposition and orientation of the filaments are uniformlyrandomly distributed. Periodic boundary conditions inboth directions are used. Upon increasing the line density� � N‘=A there is a critical threshold �c for geometricpercolation [8]. Numerical simulations [9] show that thecorrelation length � ��� �c�

�� of the incipient infinitepercolation cluster scales with a critical exponent � �4=3, identical to the value obtained for random site per-colation on a lattice [10]. Transport of scalar quantitieslike the conductivity is also in the same universality classas lattice models [11]. In order to study the transport ofnonscalar quantities such as shear stress we need tospecify how forces are transmitted between the buildingblocks of the network. In our Mikado model the buildingblocks are homogeneous elastic rods characterized by aYoung modulus E and a circular cross section of radius r.Wherever two rods intersect, they are connected by across-link with zero extensibility. In the cytoskeletonone finds a variety of linker proteins with a range ofmechanical properties [12]. Here we restrict ourselves tocross-links that either fix the relative orientation of therods (‘‘stiff cross-links’’) or allow free rotation (‘‘free

2003 The American Physical Society 108103-1

FIG. 1 (color). Typical networks atlow and high density. Dangling bonds,not contributing to the elasticity, havebeen cut off. The stress distribution isshown in false colors; the load on a fila-ment increases from blue to red. The leftpicture is for ��10, system size L�10,and an aspect ratio � � 0:0001. 99:99%of the strain energy is stored in bendingmodes. In contrast, the right pictureshows a network for ��50, L�2, and� � 0:01, where only 5% of the strainenergy is in bending modes; the remain-der is stored in compression modes. Forthe choice of units see the main text.

P H Y S I C A L R E V I E W L E T T E R S week ending5 SEPTEMBER 2003VOLUME 91, NUMBER 10

polymers, the elastic response of a stick segment betweentwo neighboring cross-links is characterized by lengthdependent force constants for compression or elongation,kcomp�‘c� � �r2E=‘c, and bending, kbend�‘c� � k?�‘c� ��3=4��r4E=‘3c. The distance between two cross-links ‘cshows a Poissonian distribution, where the average dis-tance of cross-links along a filament scales as the inverseof the line density, �‘‘c � �=� [8]. While this is a purelymechanical model that does not exhibit the temperaturedependent longitudinal linear response force constant kkof thermally fluctuating semiflexible polymers givenabove, it still captures the essential feature that the com-pressional stiffness is much larger than the bending stiff-ness kcomp�‘c�=kbend�‘c� � �4=3�‘2c=r2 � 1 for typicaldensities of the network. It does not account for stericeffects due to thermal fluctuations of the filaments, whichgive rise to the plateau modulus in solutions [13].

Consider the ground state energy of the network as afunction of the deviations of the positions of all cross-links and the rod orientations at the cross-links fromtheir initial values. For small deformations of the net-work, this function can be approximated by a quadraticform that vanishes for vanishing deviations, as — by con-struction — the undeformed network is not prestressed.The matrix representing the quadratic form can be com-puted from the network geometry and the elastic con-stants. Note that this linearized model cannot account forthe effect of buckling instabilities appearing for finitedeformations. To analyze the elastic properties of themodel network, a shear deformation respecting the peri-odic boundary conditions is enforced by demanding thatcorresponding points on the left and right boundary ofthe simulation cell undergo equal displacements while thedisplacements of corresponding points on the upper andlower boundary of the cell must agree vertically but differhorizontally by a distance � � �L, where � is the shearstrain. The orientation of the rods at corresponding pointson the boundary are required to be equal. The remainingdegrees of freedom are then allowed to relax, i.e., theharmonic approximation to the energy of the network isminimized in the presence of the constraints. The deriva-

108103-2

tive of the resulting energy of the deformed state withrespect to the strain � is proportional to the shear modu-lus. This reduces the determination of the modulus of agiven network to the solution of a linear equation.However, for interesting parameters (thin rods), the prob-lem is numerically highly unstable as we are searchingfor the lowest point of a complicated high-dimensionalvalley with extremely steep slopes but hardly varyingbase altitude. The results presented below were obtainedusing the commercially available finite element solverNastran by MSC Software.

In the following discussion we choose the rod length ‘as unit of length and �=‘3 as unit for the elastic modulus.Then the independent parameters are density �, systemsize L and aspect ratio � � r=‘ of the rods. Note that thelatter is a measure of the relative magnitude of compres-sional to bending stiffness.

We start with an analysis of the elasticity close to thepercolation threshold. For stiff cross-links we find that thepercolation threshold is the same for rigidity as for con-nectivity percolation, �c � 5:71. For free hinges a higherline density �c � 6:7 is needed for the network to be-come rigid. This agrees well with recent results, �c �6:68, for stiff fiber networks [14], where the cross-linksare fixed in space but the angles between the fibers canvary. In both cases, we find that the shear modulus Gvanishes like G� ��� �c�

� as the line density ap-proaches the critical value �c. For our numerical analysiswith finite systems we expect the shear modulus to obeythe following finite size scaling law

G � L��=� h�L= � ; (1)

where the scaling function behaves as h�x� � x�=�

and h�x� � 1 for large and small values of x � L= ,respectively. Figure 2 shows that the data collapse worksvery well for densities ranging from values close to �c upto � 20. For the data shown, L ranges from 2 to 30. Forlarger densities, systematic deviations are clearly visible,indicating that the length scale r � �‘ becomes relevanthere. We will discuss the behavior in the high densityregime in more detail below. We get the best data collapse

108103-2

P H Y S I C A L R E V I E W L E T T E R S week ending5 SEPTEMBER 2003VOLUME 91, NUMBER 10

in the critical region if we choose the values 2:4 0:2 and2:3 0:2 for the critical exponent �=� in the case of stiffcross-links and free hinges, respectively. Since the dif-ference between the exponents is within the statisticalerror, we can make no definite conclusion whether net-works with free hinges and stiff cross-links belong todifferent universality classes. The rigidity exponent � 3:15 0:2 is significantly lower than in other classes ofcontinuum percolation models, such as the ‘‘Swiss-cheesemodel,’’ where � 5 [15,16]. It is also lower than thevalue � 4 for lattice models with bond-bending forces[10,17]. Hence it seems likely that the Mikado modelconstitutes a new universality class for rigidity percola-tion. Similar results have been found in Ref. [18].

We now come back to the systematic deviations fromthe scaling law, Eq. (1), at densities above � 20. To un-derstand these better, let us have a closer look at the shearmodulus as a function of r � �‘ for densities not tooclose to the percolation threshold. In this regime the shearmodulus becomes independent of system size for moder-ately large systems; for the following results we havechosen systems satisfying L= � 200. Figure 3 showsthe shear modulus as a function of � for a series of den-sities; we have communicated a preliminary version ofthese data in Ref. [19]. Note that kbend�‘� is effectivelykept constant since we are measuring all elastic constantsin units of �=‘3. There are two strikingly different re-gimes. For high densities and/or thick rods (� * 0:1),where compressional stiffness is lower or comparable tothe bending stiffness (lower right part of Fig. 3), the shearmodulus scales linearly with the filament compressionalmodulus and the number of filaments per unit area, G�

-2 0 2 4log10(Lδρν

)

-5

0

5

10

15

log 10

h

5.566.56.76.86.977.17.27.5810152030406080100

ν=4/3 µ/ν=2.3 ρc=6.7

G=L−µ/ν

h(Lδρν)

FIG. 2 (color). Double logarithmic plot of the scaling func-tion h�x� for the shear modulus of the ‘‘Mikado model’’ withfree hinges as a function of x � Lj��j� with �� � �� �c fora series of densities � indicated in the graph. Note that for finitesystems the shear modulus is also nonzero below �c (lowerbranch in the plot). The corresponding plot for stiff hinges issimilar, but exhibits a different �c � 5:71 (see text).

108103-3

��� �c���2. Such a linear regime has also been found in

a series of studies on random fiber networks [20]. It is bynow well established that the elastic modulus in thisregime can be described quantitatively in terms of effec-tive medium models [21]. Hence, in the high line densityregime the network behaves as a homogeneously elasticmedium, dominated by the compressional modulus of theindividual filaments. As a consequence, local deforma-tions follow a macroscopic shear in an affine way. Thishas to be contrasted with the elastic behavior for slenderrods with low aspect ratios (� 10�5 for the higherdensities), where bending becomes the softer mode. Wefind an extended plateau region, which broadens signifi-cantly with lowering the line density, where the shearmodulus becomes completely independent of kcomp�‘� ���2kbend�‘� [19]. This strongly suggests that in this re-gime the macroscopic elasticity of the network is domi-nated by the bending stiffness of the filaments. Thisconclusion is corroborated by the observation that almostall of the energy stored in the deformed network is ac-counted for by transverse deformation of the rods (com-pare Fig. 1). Another remarkable feature of the plateauregime is the strong dependence of the shear modulus online density. We find G� ��� �c�

�0with a rather large

exponent �0 6:7. Figure 3 suggests the crossover scalingansatz

G � ��� �c��0g ���� �c�

�0� � 0��0=�0 ~gg��= 0� ; (2)

where we have defined a new length scale 0 ���� �c�

��0. For this ansatz to reduce to the modulus

expected in the affine region, the scaling function g�x�needs to scale as g�x� � x�2 for x � 1 and the exponentsneed to obey the scaling relation �0 � 2�0 � 1. In theplateau regime, g�x� is expected to be constant. As shownin Fig. 4, we obtain an excellent scaling collapse for overalmost 8 orders of magnitude in the scaling variable x ��= 0 using �0 � 2:83 or equivalently �0 � 6:67 and thecritical line density �c 5:71 associated with connectiv-ity percolation. Additionally, the scaling function g�x�

10-8

10-6

10-4

10-2

100

α

100

105

1010

G

ρ = 15ρ = 20ρ = 30ρ = 40ρ = 60ρ = 80ρ = 100ρ = 150

FIG. 3 (color). Double logarithmic plot of the shear modulusG as a function of � for fixed kbend�‘�. Data shown for freehinges.

108103-3

-4 -2 0 2 4 6log10(α δρν

)

-12

-10

-8

-6

-4

-2lo

g 10g

ρ = 15ρ = 20ρ = 30ρ = 40ρ = 60ρ = 80ρ = 100ρ = 150

ν=2.83 µ=6.67 ρc=5.71G=δρµ

g(α δρν)

FIG. 4 (color). Scaling plot of the shear modulus for freehinges for a series of densities above � � 15 indicated in thegraph (same data as in Fig. 3). Data collapse to the crossoverscaling form, Eq. (2), is obtained with �0 � 2:83. Note thathere and in all other figures the unit of length is ‘ and the unitof the shear modulus G is �=‘3.

P H Y S I C A L R E V I E W L E T T E R S week ending5 SEPTEMBER 2003VOLUME 91, NUMBER 10

displays the expected behavior. Meeting both of theserequirements is highly nontrivial, and gives strong evi-dence for the anomalous scaling law in Eq. (2).

The existence of such a broad scaling regime far fromthe percolation threshold is a surprising and intriguingfeature of stiff polymer networks. Its physical origin isdistinct from the critical scaling regime, and governed bya new length scale 0. While the geometrical significanceof 0 is yet unclear, one may speculate that the anomalousscaling behavior is a subtle consequence of the interplaybetween quenched randomness of the network structureand long-range correlation effects induced by the stiff-ness of the filaments. An immediate consequence of thescaling form, Eq. (2), is the existence of a crossover linedensity �cross scaling as ‘�cross � ��1=�0

, where we havereintroduced units of length ‘. This implies that increas-ing filament length at constant line density drives thesystem towards the affine regime, in accord with Ref. [18].

While these results for an idealized two-dimensionalmodel are certainly not straightforwardly applicable tothree-dimensional cytoskeletal networks, one may stilltry to get an idea of the scales involved. We expect thatnetwork densities can be compared roughly by using theaverage distance ‘c between intersections as a measure: Acytoskeletal network might have ‘c 0:1�m with fila-ment lengths of 2�m and filament radii of 4 nm. Thesevalues correspond to a two-dimensional line density of�20 and an aspect ratio of �0:002, which wouldplace a typical actin network in the bending dominatedintermediate regime at T � 0. For T > 0, a faithful treat-ment of thermal fluctuations of the network would re-quire the inclusion of the polymer configuration asadditional degrees of freedom. However, in a first approxi-mation, the fluctuations of the polymers between cross-links could be represented by replacing the longitudinal

108103-4

force constant kcomp used in the present analysis by thethermal linear response coefficient kk discussed in theintroduction augmented by numerical prefactors to reflectthe local boundary conditions. As the dependence of kk onthe random distance between cross-links is quite differentfrom that of kcomp, this may well lead to different scalingbehavior. Additional contributions to the modulus areexpected when thermal fluctuation of the cross-link po-sitions are taken into account [22]. Understanding thefull complexity of cytoskeletal networks certainly meritsfurther theoretical and experimental work. Future inves-tigations may among many other questions want to ad-dress three-dimensional systems, polydispersity, thermalfluctuations, or even the kinetics of the cross-linkingmolecules.

We thank M. Alava and K. Kroy for discussions, andP. Benetatos for a critical reading of the manuscript.

[1] B. Alberts et al., Molecular Biology of the Cell (GarlandPublishers, New York, 1994), 3rd ed.

[2] J. Wilhelm and E. Frey, Phys. Rev. Lett. 77, 2581 (1996).[3] L. LeGoff, O. Hallatschek, E. Frey, and F. Amblard,

Phys. Rev. Lett. 89, 258101 (2002).[4] M. Doi and S. F. Edwards, The Theory of Polymer

Dynamics (Clarendon Press, Oxford, 1986).[5] K. Kroy and E. Frey, Phys. Rev. Lett. 77, 306 (1996).[6] F. MacKintosh, J. Kas, and P. Janmey, Phys. Rev. Lett. 75,

4425 (1995).[7] R. Satcher and C. Dewey, Biophys. J. 71, 109 (1996).[8] G. Pike and C. Seager, Phys. Rev. B 10, 1421 (1974).[9] Y. Leroyer and E. Pommiers, Phys. Rev. B 50, 2795

(1994).[10] D. Stauffer and A. Aharony, Introduction to Percolation

Theory (Taylor & Francis, London, 1994), 2nd ed.[11] I. Balberg, N. Binenbaum, and N. Wagner, Phys. Rev.

Lett. 52, 1465 (1984).[12] L. Limozin and E. Sackmann, Phys. Rev. Lett. 89, 168103

(2002).[13] B. Hinner et al., Phys. Rev. Lett. 81, 2614 (1998).[14] M. Latva-Kokko and J. Timonen, Phys. Rev. E 64, 066117

(2001).[15] S. Feng, B. I. Halperin, and P. N. Sen, Phys. Rev. B 35,

197 (1987).[16] L. Benguigui, Phys. Rev. B 34, 8176 (1986).[17] S. Arbabi and M. Sahimi, Phys. Rev. B 47, 703 (1993).[18] D. A. Head, A. J. Levine, and F. C. MacKintosh, Phys.

Rev. Lett. (to be published).[19] E. Frey, K. Kroy, J. Wilhelm, and E. Sackmann, in

Dynamical Networks in Physics and Biology, edited byG. Forgacs and D. Beysens (Springer-Verlag, Berlin,1998).

[20] V. Raisanen, M. Alava, K. Niskanen, and R. Nieminen,J. Mater. Res. 12, 2725 (1997).

[21] J. Astrom et al., Phys. Rev. E 61, 5550 (2000).[22] M. Plischke and B. Joos, Phys. Rev. Lett. 80, 4907 (1998).

108103-4