Alignment and Nonlinear Elasticity in Biopolymer Gels
Alignment and Nonlinear Elasticity in Biopolymer Gels Jingchen Feng,1 Xiaoming Mao,2 and Leonard M. Sander3 1
arXiv:1402.2998v1 [cond-mat.soft] 12 Feb 2014
Bioengineering Department and Center for Theoretical Biological Physics, Rice University, Houston TX, 77251-1892 2 Department of Physics, University of Michigan, Ann Arbor MI 48109-1040 3 Physics & Complex Systems, University of Michigan, Ann Arbor MI 48109-1040 (Dated: February 14, 2014) We present a Landau type theory for the non-linear elasticity of biopolymer gels with a part of the order parameter describing induced nematic order. Our point of view is that all of the non-linear elastic behavior of these materials can be attributed to the onset of fiber alignment with induced strain. We suggest an application to contact guidance of cell motility in tissue. We compare our theory to simulation of a disordered lattice model for biopolymers. We treat homogeneous deformations such as shear and extension, and also treat the case of a localized perturbation which is intended to be a simple model for a contacting cell in a medium. PACS numbers: 62.20.D-, 87.10.Pq, 87.17.Jj
I.
INTRODUCTION
Bio-polymer gels are complicated materials consisting of arrays of cross-linked polymer chains. Almost all of these materials show non-linear elasticity, usually strainstiffening. Typically, the shear modulus increases by an order of magnitude under applied strain [1]. These materials are important in biology. For example, in tissue there is almost always a gel, the extra cellular matrix (ECM) which gives tissue its structure. The biopolymer collagen-I is the most common constituent of the ECM. When cells migrate within tissue they crawl by attaching to the ECM. This is important in a variety of circumstances: wound healing and cancer invasion are important examples [2]. When cells move in this fashion they deform the ECM quite a lot by pulling on it. However deformation also affects cell motility. In particular, if the ECM is aligned, cells tend to move along the aligned direction – this is called contact guidance [3, 4]. Since cells deform and align the ECM, and their motility is affected by this deformation, mechanically mediated cell-cell interactions are to be expected [2, 4]. The purpose of this paper is to formulate a theory for the non-linear elasticity in biopolymers to give a framework for studying the mechanical interaction between moving cells. The phenomenon of contact guidance is known from experiment, but no real understanding of its microscopic nature is available. An obvious step towards such understanding is to quantify the alignment of the fibers in the ECM. To do this, we take over a concept from the theory of liquid crystals, the nematic order parameter. Of course, collagen-I does not have a nematic instability, but it can be aligned by stress to give a non-zero value of the nematic tensor, Q, describing the local directions of the polymers, 1 Q(r) ≡ hˆ ν νˆ − Ii. d
(1.1)
Here νˆ is the unit vector pointing along the orientation
of the polymers within the volume element at r, and the average h· · · i is over all fibers within the volume element, weighted by the fiber length. In what follows we will use a scalar, q, 0 ≤ q ≤ 1, to characterize the strength of the alignment. We define q = [d/(d − 1)]λm , where λm is the largest eigenvalue of Q. For example, in two dimensions, q = hcos(2θ)i, where θ is the direction of the polymer chain. Now we can speculate about contact guidance. Suppose we have a cell moving randomly so that it has an effective diffusion coefficient, D◦ if there is no alignment. Then, based on symmetry, the simplest way to account for alignment would be: D = D◦ (1 + αQ),
(1.2)
where α is a coefficient which could be measured. This idea motivates us to think about the the nonlinear elasticity of collagen-I using Q as a central ingredient. Purely from the point of view of elasticity, this idea is reasonable. One explanation [5, 6] for strainstiffening in biopolymers like collagen-I is that for small strains the elasticity is dominated by the (small) bending modulus – different parts of the disordered material turn with respect to one another so that the deformation is non-affine. For large strains the material is aligned, and must stretch; the large stretching modulus determines the elastic response; the deformation is affine as in ordinary elasticity theory. It is plausible that the transition occurs when q achieves an appreciable value. We describe biopolymers composed of linear elements whose collective response is non-linear: the non-linearity arises from a “hidden” variable, the alignment. In what follows we pursue this notion by formulating a Landau-deGennes theory for biopolymers along the lines suggested by [7] for nematic elastomers. We test the framework by fitting our theory to a disordered lattice model [8, 9] to determine the coefficients. Then we insert a “cell”, i.e., a localized source of deformation, to begin to address the questions alluded to above.
2 II.
MODEL
We model the nonlinear elastic energy of bio-polymer gels with the following effective free energy: Z h λ (Tr v)2 + µ ˜Tr v2 + g(Tr v)3 F = 2 i − tTr(v · Q) + V (Q) dd r, (2.1) where λ, µ ˜ are the Lam´e coefficients and v is the (left) Cauchy-Green strain tensor. The strain tensor depends on the displacement field u(r): vij ≡ (∂i uj + ∂j ui + ∂l ui ∂l uj ) /2,
(2.2)
It transforms as a tensor in the deformed state. The tilde on µ ˜ means that it is a “bare” value. The observed shear modulus, µ, is gotten by renormalizing µ ˜ by coupling to Q, as we will see. The potential for the nematic order, V (Q), is the energy cost to align the network against the constraints. The non-linear part of the energy cost due to volume change is given by g. The presence of this term departs from our view that non-linearity should be attributed to alignment. Note that the term vanishes for shear deformations. We will show below that a proper account of disorder (which is certainly present in biopolymers in the regime that interests us) can give a more fundamental view of non-linearity in volume deformation. For the moment, we retain g as an effective parameter. The term −tTr(v · Q) is the leading order coupling between strain and alignment allowed by symmetry. The order parameter, Q, which describes the alignment induced by strain, transforms as a tensor in the deformed state so that it has the same symmetry as v. Their contraction is a scalar that can enter the free energy. Because Q is traceless, in leading order the alignment responds to pure shear and not to hydrostatic deformations. For small deformations, linear elasticity, we keep the leading order term in the expansion of V , Vl (Q) =
A TrQ2 . 2
(2.3)
˜ is Minimizing F for fixed v gives Q = (t/A)˜ v , where v the traceless part of v. Thus strain induced alignment is a linear response for small deformations. Since this relation is determined by symmetry, we expect it to hold for all biopolymer gels. The alignment tensor can be eliminated from F , recovering linear elasticity: Z h i λ (2.4) (Trv)2 + µTrv2 dd r, Fl = 2 in which the effective shear modulus: µ=µ ˜−
t2 , 2A
is reduced from the bare value.
(2.5)
FIG. 1. A small area of the network before (left) and after (right) extension along the horizontal axis; dots denote crosslinks. The polymer fibers strongly align along the direction of the extension, and bend according to the constraints at the crosslinks (an example of such a bent fiber is on the left). The bending in this local area is exhausted when a large fraction of the fibers is parallel to the direction of extension, such as the left section of the bent fiber, and the other one as well. Any further extension necessarily involves stretching the fibers, which costs more energy.
At larger deformations, nonlinear terms in V start to dominate and Q becomes smaller than its linear response value. Therefore the effective shear modulus increases from its value in the linear regime, giving strain stiffening. In the extreme case of Q reaching a saturated value Qmax independent of v, we get µ ˜ ≃ µ. The nonlinearity in V (Q) can be understood intuitively. As we mentioned, at small strain the deformation in a dilute bio-polymer network consists mostly of bending of fibers. At greater strain, the deformation crosses over to stretching of fibers, because the easy bending modes are “exhausted,” and the increase of Q saturates, as shown in Fig. 1. In contrast, for denser networks, even for small deformations the system is stretching dominated, and the crossover disappears. We expect the nonlinearity to be strong for dilute networks and weak for dense ones. This is verified in simulations, see below.
III.
TRIANGULAR LATTICE MODEL
To calibrate the free energy of Eq. (2.1) we apply it to a simplified lattice model for bio-polymer gels based on a triangular lattice of lattice constant a, in which each bond is present with a probability p [8, 9]. Straight lines in this lattice, which have average length (1 − p)−1 , are identified as fibers with stretching stiffness k and bending stiffness κ. The lattice sites are freely rotating crosslinks. The Hamiltonian is: E=
k X κ X 2 gij (|Rij | − a)2 + gij gjk ∆θijk , (3.1) 2a 2a hiji
hijki
where gij = 1 for bonds that are present, and 0 for removed ones. The first term is the stretching energy; |Rij | is the distance between sites i and j in the deformed state. The second term is bending; hijki labels three consecutive sites along a straight line in the reference state, and ∆θijk the change of angle determined by them. The linear elasticity of this model is largely controlled by the central force isostatic point at pc ≃ 2/3. For
3
V (Q) =
A C D Tr Q2 + Tr Q4 + Tr Q6 2 4! 6!
(3.2)
The odd terms vanish by symmetry in two dimensions. Since we are using an expansion of this type, we should not expect the scheme to work deep in the non-linear regime; we expect quantitative results for relatively small q, say q ≤ 0.3. At this level there are 6 parameters in the theory, {λ, µ ˜, g, t, A, C, D}. Consider the linear regime, which is characterized by three independent slopes for stress-strain curves in shear and in hydrostatic extension, and the nematic-strain curve in shear; the corresponding slopes are the shear modulus, µ = µ ˜ − t2 /(2A), the bulk modulus, K = µ + λ, and t/A. By symmetry, there cannot be any average induced nematic order in the hydrostatic extension case, and this is what we find. The rest of the parameters are gotten in the nonlinear regime. We fit them using simulation for hydrostatic ex-
10
4
Stress
103
κ =10−3 ,Shearing Simulation: p=0.80 Fitting with Theory: p=0.80 Simulation: p=0.60 Fitting with Theory: p=0.60
0.4 0.3
1
100
10
4
10-2
10-1
Strain
0.0 0.0
100
κ =10−3 ,Hydrostatic deformation
0.2
0.4
Strain
0.6
0.8
1.0
500
Simulation: p=0.80 Fitting with Theory: p=0.80 Simulation: p=0.60 Fitting with Theory: p=0.60
400
Nonlinear term C
105
0.2 0.1
10-1 10-2 -3 10
Simulation: p=0.80 Fitting with Theory: p=0.80 Simulation: p=0.60 Fitting with Theory: p=0.60
0.5
102
10
κ =10−3 ,Shearing
0.6
Alignment order q
105
103
Stress
p > pc the deformations are mostly affine; below pc disorder induced by the removal of bonds leads to nonaffine response. This is because the bending stiffness of the fibers is much smaller than the stretching stiffness: κ/(ka2 ) ≪ 1. When κ = 0, the system reduces to a central-force lattice, which exhibits a rigidity percolation transition at pc . The network with weak bending can be viewed as central-force network with a relevant perturbation of bending stiffness, and as a result the elasticity of such networks can be thought of as a crossover at the central-force isostatic point. The crossover has been studied in [9]. We studied this lattice numerically in the nonlinear regime using 64x64 lattices. We applied three types of homogeneous deformations: pure shear, hydrostatic extension, and simple extension and measured the stress and the nematic order as functions of strain, γ, at various values of p and κ/(ka2 ). For pure shear, some results are shown in Fig. 2. We have observed that, in agreement with our theory, the nematic order tensor Q does have the form η(γ)˜ v where η(γ) is a scalar. Thus the orientation of the nematic order is determined and we can use q = hcos 2θi, to characterize the strength of the alignment. The alignment is a nontrivial function of the strain as we see in the Figure. The slope is t/A for small strain from the discussion above. These results are consistent with the picture outlined above. Shear-stiffening is strong for p ≪ pc and vanishes for p > pc . The characteristic strain, γ ∗ , where the strain-stiffening takes place is large at small p and vanishes near pc . These observations are consistent with the notion of the exhaustion of bending modes as the network enters the stretching dominated regime. They are related to stiffening in jammed packings, [10], where scaling γ ∗ ∼ |∆z| occurs. Here ∆z is the extra coordination above isostaticity, equivalent to p − pc in our model. To make the comparison with our lattice simulation quantitative, we use Eq. (2.1) in two dimensions, and expand the potential V , up to sixth order in Tr Q:
300
102
200
101
10
10
100
0
-1
10-3
10-2
Strain
10-1
100
0 55
60
65
70
75
80
Probability p
(c)
85
90
95
100
(d)
FIG. 2. Top left: Stress-strain curves, log scale, for shear in lattices with p = 0.6, 0.8, with κ/(ka2 ) = 10−4 showing the qualitative change as p passes through the isostatic point, pc =2/3. For the dilute lattice, p < pc , there is alignmentinduced shear-stiffening. Top right: the strength of the nematic order as a function of strain. The theoretical curves are fits to Eq. (2.1) as discussed in the text. The data points are from simulation and dashed lines are the fitted results. Bottom left: Fits to the elastic energy from Eq. (2.1) for hydrostatic deformation as a function of strain. The stress is the derivative of this curve. The data points are from simulation and dashed lines are the fitted results. Bottom right: the fitted parameter C as a function of p.
tension for g and pure shear for the rest, separately via least-square fitting. The comparison between simulation and the fit are shown in Figure 2. As mentioned above, strain-stiffening depends strongly on p. This is captured by the parameter C, which controls where the potential becomes strongly nonlinear and increases the shear modulus. The fitted value of C as a function of p is shown in Fig. 2, bottom right. Quite remarkably, and in accord with our interpretation, there is a sharp peak at pc , showing that γ ∗ vanishes as p → p− c . We used the fitted parameters to calculate the stressstrain and nematic-strain curves for simple extension. In Fig. 3 we compare to simulation. The agreement is fairly good, except for q at large alignments; c.f. the discussion below. IV.
CELLS IN ECM
With our results for homogeneous deformations in hand, we turn to a schematic representation of a cell in the ECM. We can make a very simplified picture by cutting a circular hole in the material at some point, and applying a negative pressure to the exposed surface to represent the contractile force of the cell on the material [11]. This is done in our Landau theory by applying
4 κ =10−4 ,Simple Extension Simulation: p=0.80 Prediction from Theory: p=0.80 Simulation: p=0.60 Prediction from Theory: p=0.60
103 10
2
101
0.8 0.7
Alignment order q
0.6
Energy
100
0.5
10-1 10-2
0.4
10-3
0.3
10-4 10-5
0.2 0.1
10-6 10
κ =10−4 ,Simple Extension Simulation: p=0.80 Prediction from Theory: p=0.80 Simulation: p=0.60 Prediction from Theory: p=0.60
-7
10-3
10-2
10-1
Strain
0.0 0.0
100
0.2
0.4
(a)
Strain
0.6
0.8
1.0
(b)
FIG. 3. Comparison of the simulation results (data points) and theory predictions (dashed lines) for the simple extension case at p = 0.6, 0.8. The theory uses parameters fitted from the pure shear and the hydrostatic extension cases.
displacement u
0.12 0.10
Prediction Far field solution Simulation
0.20 Alignment order q
Prediction Far field solution Simulation
0.10
0.80
0.05
0.60
0.02
0.40 4
6
8
10
Distance from center node r
12 14
6
8
10
Distance from center node r
(a)
12
14
(b)
FIG. 4. Comparison of the simulation results (data points) and theory predictions (solid lines) for a model for a cell in ECM. The theory uses the parameters derived for the homogeneous deformation case.
boundary conditions. For the lattice model, we apply forces to nodes of the lattice at radius of several lattice spacings and relax the lattice. To use Eq. (2.1) in this case, we note that we have two variables which depend on the distance, r, from the center of the hole: they are the radial component of the displacement, u(r), and the strength of the radial alignment, q(r). If we minimize F with respect to u, q we find two coupled equations: 2u′ 2t 2u − 2 = r r λ + 4µ 2t ′ u u + = V ′ (q). 3 r
u′′ +
q′ q + 3 r
Alignment localization of this type is captured by our theory through the nonlinear potential V . For small r, q is large and the nonlinear potential leads to a weaker dependence of q(r) on r; that is, alignment is more or less constant near the cell [2]. We used the parameters from the homogeneous deformation case to calculate u(r), q(r) by numerically solving Eq. (4.1); see Fig. 4. We also studied the same problem using simulation by choosing a node and pulling in on nodes a few lattice spacings away. The two kinds of results agree, as shown in Fig. 4. V.
SUMMARY AND EXTENSIONS
In summary, we have given a Landau-type theory for elastic non-linearity in biopolymer gels. The central focus is to introduce the alignment tensor as a part of the order parameter. This corresponds to the notion that for many biopolymers, particularly collagen-I, the non-linear behavior is due to switching of deformation from bending to stretching as the fibers become parallel. It also should, eventually, allow an attempt to treat multi-cell systems using Eq. (1.2) as a starting point. Of course, this calculation is only the beginning of a treatment for cells. We should probably represent the cell as a force dipole rather than a spherical object: moving cells are polarized. And we need to put more than one cell in the system to see the interactions. We can suggest a more satisfying alternative to the intrinsic non-linearity in the system represented by g. We prefer it because our point of view is that all of the complex behavior results from changes in alignment. So far, Eq. (2.1) has no coupling to volume deformation – i.e. only the traceless part of v induces alignment. However, even for hydrostatic expansion, there will be local alignment and weak parts of the disordered system will rotate. We propose to treat this by introducing a random tensor field, h(r) with hhi = 0 which accounts for disorder; see also [13]. In the lattice model it is related to the missing bonds. Then we can add a term to the integrand in Eq. (2.1) of the form −Tr vTr (h · Q). In the same way as above, in the linear regime this renormalizes the first Lam´e coefficient (and hence the bulk modulus, K): λ → λ − Trhh2 i/A.
(4.1)
For large r the deformation is small, we expect the system to be linear. If we use V ≈ Aq 2 /2 we recover the familiar result for the deformation in two dimensions [12] u(r) ∼ 1/r. The second equation of (4.1) shows that q ∼ 1/r2 . As r decreases the deformations get larger, and we enter the nonlinear elasticity regime. As discussed in Ref. [2] there will be a critical radius within which nonlinear effects are important, and alignment is large.
As in the previous case the renormalization decreases in the non-linear regime. The leading order effect is the same as the term in g that we added by hand. ACKNOWLEDGMENTS
We would like to thank F. MacKintosh, C. Broedersz, M. Das, H. Levine, and E. Ben-Jacob for useful conversations. J-C Feng is supported by the National Science Foundation Center for Theoretical Biological Physics (Grant PHY-1308264).
5
[1] C. Storm, J. Pastore, F. MacKintosh, T. Lubensky, and P. Janmey, Nature 435, 191 (2005). [2] L. M. Sander, Journal of biomechanical engineering 135, 71006 (2013). [3] P. A. Agudelo-Garcia, J. K. De Jesus, S. P. Williams, M. O. Nowicki, E. A. Chiocca, S. Liyanarachchi, P.-K. Li, J. J. Lannutti, J. K. Johnson, and S. E. Lawler, Neoplasia (New York, NY) 13, 831 (2011). [4] D. Vader, A. Kabla, D. Weitz, and L. Mahadevan, PloS one 4 (2009). [5] P. Onck, T. Koeman, T. van Dillen, and E. van der Giessen, Physical Review Letters 95, 178102 (2005). [6] A. M. Stein, D. A. Vader, D. A. Weitz, and L. M. Sander, Complexity 16, 22 (2011). [7] T. Lubensky, R. Mukhopadhyay, L. Radzihovsky, and
X. Xing, Physical Review E 66, 011702 (2002). [8] M. Das, F. C. MacKintosh, and A. J. Levine, Phys. Rev. Lett. 99, 038101 (2007). [9] C. P. Broedersz, X. Mao, T. C. Lubensky, and F. C. MacKintosh, Nature Physics 7, 983 (2011). [10] M. Wyart, H. Liang, A. Kabla, and L. Mahadevan, Phys. Rev. Lett. 101, 215501 (2008). [11] M. Buenemann, H. Levine, W.-J. Rappel, and L. M. Sander, Biophysj 99, 50 (2010). [12] L. D. Landau, E. M. Lifshitz, A. M. Kosevich, and L. P. Pitaevskii, Theory of elasticity (Pergamon Press, Oxford; New York, 1986). [13] M. Sheinman, C. P. Broedersz, and F. C. MacKintosh, Phys. Rev. E 85, 021801 (2012).
arXiv:1402.2998v1 [cond-mat.soft] 12 Feb 2014
Bioengineering Department and Center for Theoretical Biological Physics, Rice University, Houston TX, 77251-1892 2 Department of Physics, University of Michigan, Ann Arbor MI 48109-1040 3 Physics & Complex Systems, University of Michigan, Ann Arbor MI 48109-1040 (Dated: February 14, 2014) We present a Landau type theory for the non-linear elasticity of biopolymer gels with a part of the order parameter describing induced nematic order. Our point of view is that all of the non-linear elastic behavior of these materials can be attributed to the onset of fiber alignment with induced strain. We suggest an application to contact guidance of cell motility in tissue. We compare our theory to simulation of a disordered lattice model for biopolymers. We treat homogeneous deformations such as shear and extension, and also treat the case of a localized perturbation which is intended to be a simple model for a contacting cell in a medium. PACS numbers: 62.20.D-, 87.10.Pq, 87.17.Jj
I.
INTRODUCTION
Bio-polymer gels are complicated materials consisting of arrays of cross-linked polymer chains. Almost all of these materials show non-linear elasticity, usually strainstiffening. Typically, the shear modulus increases by an order of magnitude under applied strain [1]. These materials are important in biology. For example, in tissue there is almost always a gel, the extra cellular matrix (ECM) which gives tissue its structure. The biopolymer collagen-I is the most common constituent of the ECM. When cells migrate within tissue they crawl by attaching to the ECM. This is important in a variety of circumstances: wound healing and cancer invasion are important examples [2]. When cells move in this fashion they deform the ECM quite a lot by pulling on it. However deformation also affects cell motility. In particular, if the ECM is aligned, cells tend to move along the aligned direction – this is called contact guidance [3, 4]. Since cells deform and align the ECM, and their motility is affected by this deformation, mechanically mediated cell-cell interactions are to be expected [2, 4]. The purpose of this paper is to formulate a theory for the non-linear elasticity in biopolymers to give a framework for studying the mechanical interaction between moving cells. The phenomenon of contact guidance is known from experiment, but no real understanding of its microscopic nature is available. An obvious step towards such understanding is to quantify the alignment of the fibers in the ECM. To do this, we take over a concept from the theory of liquid crystals, the nematic order parameter. Of course, collagen-I does not have a nematic instability, but it can be aligned by stress to give a non-zero value of the nematic tensor, Q, describing the local directions of the polymers, 1 Q(r) ≡ hˆ ν νˆ − Ii. d
(1.1)
Here νˆ is the unit vector pointing along the orientation
of the polymers within the volume element at r, and the average h· · · i is over all fibers within the volume element, weighted by the fiber length. In what follows we will use a scalar, q, 0 ≤ q ≤ 1, to characterize the strength of the alignment. We define q = [d/(d − 1)]λm , where λm is the largest eigenvalue of Q. For example, in two dimensions, q = hcos(2θ)i, where θ is the direction of the polymer chain. Now we can speculate about contact guidance. Suppose we have a cell moving randomly so that it has an effective diffusion coefficient, D◦ if there is no alignment. Then, based on symmetry, the simplest way to account for alignment would be: D = D◦ (1 + αQ),
(1.2)
where α is a coefficient which could be measured. This idea motivates us to think about the the nonlinear elasticity of collagen-I using Q as a central ingredient. Purely from the point of view of elasticity, this idea is reasonable. One explanation [5, 6] for strainstiffening in biopolymers like collagen-I is that for small strains the elasticity is dominated by the (small) bending modulus – different parts of the disordered material turn with respect to one another so that the deformation is non-affine. For large strains the material is aligned, and must stretch; the large stretching modulus determines the elastic response; the deformation is affine as in ordinary elasticity theory. It is plausible that the transition occurs when q achieves an appreciable value. We describe biopolymers composed of linear elements whose collective response is non-linear: the non-linearity arises from a “hidden” variable, the alignment. In what follows we pursue this notion by formulating a Landau-deGennes theory for biopolymers along the lines suggested by [7] for nematic elastomers. We test the framework by fitting our theory to a disordered lattice model [8, 9] to determine the coefficients. Then we insert a “cell”, i.e., a localized source of deformation, to begin to address the questions alluded to above.
2 II.
MODEL
We model the nonlinear elastic energy of bio-polymer gels with the following effective free energy: Z h λ (Tr v)2 + µ ˜Tr v2 + g(Tr v)3 F = 2 i − tTr(v · Q) + V (Q) dd r, (2.1) where λ, µ ˜ are the Lam´e coefficients and v is the (left) Cauchy-Green strain tensor. The strain tensor depends on the displacement field u(r): vij ≡ (∂i uj + ∂j ui + ∂l ui ∂l uj ) /2,
(2.2)
It transforms as a tensor in the deformed state. The tilde on µ ˜ means that it is a “bare” value. The observed shear modulus, µ, is gotten by renormalizing µ ˜ by coupling to Q, as we will see. The potential for the nematic order, V (Q), is the energy cost to align the network against the constraints. The non-linear part of the energy cost due to volume change is given by g. The presence of this term departs from our view that non-linearity should be attributed to alignment. Note that the term vanishes for shear deformations. We will show below that a proper account of disorder (which is certainly present in biopolymers in the regime that interests us) can give a more fundamental view of non-linearity in volume deformation. For the moment, we retain g as an effective parameter. The term −tTr(v · Q) is the leading order coupling between strain and alignment allowed by symmetry. The order parameter, Q, which describes the alignment induced by strain, transforms as a tensor in the deformed state so that it has the same symmetry as v. Their contraction is a scalar that can enter the free energy. Because Q is traceless, in leading order the alignment responds to pure shear and not to hydrostatic deformations. For small deformations, linear elasticity, we keep the leading order term in the expansion of V , Vl (Q) =
A TrQ2 . 2
(2.3)
˜ is Minimizing F for fixed v gives Q = (t/A)˜ v , where v the traceless part of v. Thus strain induced alignment is a linear response for small deformations. Since this relation is determined by symmetry, we expect it to hold for all biopolymer gels. The alignment tensor can be eliminated from F , recovering linear elasticity: Z h i λ (2.4) (Trv)2 + µTrv2 dd r, Fl = 2 in which the effective shear modulus: µ=µ ˜−
t2 , 2A
is reduced from the bare value.
(2.5)
FIG. 1. A small area of the network before (left) and after (right) extension along the horizontal axis; dots denote crosslinks. The polymer fibers strongly align along the direction of the extension, and bend according to the constraints at the crosslinks (an example of such a bent fiber is on the left). The bending in this local area is exhausted when a large fraction of the fibers is parallel to the direction of extension, such as the left section of the bent fiber, and the other one as well. Any further extension necessarily involves stretching the fibers, which costs more energy.
At larger deformations, nonlinear terms in V start to dominate and Q becomes smaller than its linear response value. Therefore the effective shear modulus increases from its value in the linear regime, giving strain stiffening. In the extreme case of Q reaching a saturated value Qmax independent of v, we get µ ˜ ≃ µ. The nonlinearity in V (Q) can be understood intuitively. As we mentioned, at small strain the deformation in a dilute bio-polymer network consists mostly of bending of fibers. At greater strain, the deformation crosses over to stretching of fibers, because the easy bending modes are “exhausted,” and the increase of Q saturates, as shown in Fig. 1. In contrast, for denser networks, even for small deformations the system is stretching dominated, and the crossover disappears. We expect the nonlinearity to be strong for dilute networks and weak for dense ones. This is verified in simulations, see below.
III.
TRIANGULAR LATTICE MODEL
To calibrate the free energy of Eq. (2.1) we apply it to a simplified lattice model for bio-polymer gels based on a triangular lattice of lattice constant a, in which each bond is present with a probability p [8, 9]. Straight lines in this lattice, which have average length (1 − p)−1 , are identified as fibers with stretching stiffness k and bending stiffness κ. The lattice sites are freely rotating crosslinks. The Hamiltonian is: E=
k X κ X 2 gij (|Rij | − a)2 + gij gjk ∆θijk , (3.1) 2a 2a hiji
hijki
where gij = 1 for bonds that are present, and 0 for removed ones. The first term is the stretching energy; |Rij | is the distance between sites i and j in the deformed state. The second term is bending; hijki labels three consecutive sites along a straight line in the reference state, and ∆θijk the change of angle determined by them. The linear elasticity of this model is largely controlled by the central force isostatic point at pc ≃ 2/3. For
3
V (Q) =
A C D Tr Q2 + Tr Q4 + Tr Q6 2 4! 6!
(3.2)
The odd terms vanish by symmetry in two dimensions. Since we are using an expansion of this type, we should not expect the scheme to work deep in the non-linear regime; we expect quantitative results for relatively small q, say q ≤ 0.3. At this level there are 6 parameters in the theory, {λ, µ ˜, g, t, A, C, D}. Consider the linear regime, which is characterized by three independent slopes for stress-strain curves in shear and in hydrostatic extension, and the nematic-strain curve in shear; the corresponding slopes are the shear modulus, µ = µ ˜ − t2 /(2A), the bulk modulus, K = µ + λ, and t/A. By symmetry, there cannot be any average induced nematic order in the hydrostatic extension case, and this is what we find. The rest of the parameters are gotten in the nonlinear regime. We fit them using simulation for hydrostatic ex-
10
4
Stress
103
κ =10−3 ,Shearing Simulation: p=0.80 Fitting with Theory: p=0.80 Simulation: p=0.60 Fitting with Theory: p=0.60
0.4 0.3
1
100
10
4
10-2
10-1
Strain
0.0 0.0
100
κ =10−3 ,Hydrostatic deformation
0.2
0.4
Strain
0.6
0.8
1.0
500
Simulation: p=0.80 Fitting with Theory: p=0.80 Simulation: p=0.60 Fitting with Theory: p=0.60
400
Nonlinear term C
105
0.2 0.1
10-1 10-2 -3 10
Simulation: p=0.80 Fitting with Theory: p=0.80 Simulation: p=0.60 Fitting with Theory: p=0.60
0.5
102
10
κ =10−3 ,Shearing
0.6
Alignment order q
105
103
Stress
p > pc the deformations are mostly affine; below pc disorder induced by the removal of bonds leads to nonaffine response. This is because the bending stiffness of the fibers is much smaller than the stretching stiffness: κ/(ka2 ) ≪ 1. When κ = 0, the system reduces to a central-force lattice, which exhibits a rigidity percolation transition at pc . The network with weak bending can be viewed as central-force network with a relevant perturbation of bending stiffness, and as a result the elasticity of such networks can be thought of as a crossover at the central-force isostatic point. The crossover has been studied in [9]. We studied this lattice numerically in the nonlinear regime using 64x64 lattices. We applied three types of homogeneous deformations: pure shear, hydrostatic extension, and simple extension and measured the stress and the nematic order as functions of strain, γ, at various values of p and κ/(ka2 ). For pure shear, some results are shown in Fig. 2. We have observed that, in agreement with our theory, the nematic order tensor Q does have the form η(γ)˜ v where η(γ) is a scalar. Thus the orientation of the nematic order is determined and we can use q = hcos 2θi, to characterize the strength of the alignment. The alignment is a nontrivial function of the strain as we see in the Figure. The slope is t/A for small strain from the discussion above. These results are consistent with the picture outlined above. Shear-stiffening is strong for p ≪ pc and vanishes for p > pc . The characteristic strain, γ ∗ , where the strain-stiffening takes place is large at small p and vanishes near pc . These observations are consistent with the notion of the exhaustion of bending modes as the network enters the stretching dominated regime. They are related to stiffening in jammed packings, [10], where scaling γ ∗ ∼ |∆z| occurs. Here ∆z is the extra coordination above isostaticity, equivalent to p − pc in our model. To make the comparison with our lattice simulation quantitative, we use Eq. (2.1) in two dimensions, and expand the potential V , up to sixth order in Tr Q:
300
102
200
101
10
10
100
0
-1
10-3
10-2
Strain
10-1
100
0 55
60
65
70
75
80
Probability p
(c)
85
90
95
100
(d)
FIG. 2. Top left: Stress-strain curves, log scale, for shear in lattices with p = 0.6, 0.8, with κ/(ka2 ) = 10−4 showing the qualitative change as p passes through the isostatic point, pc =2/3. For the dilute lattice, p < pc , there is alignmentinduced shear-stiffening. Top right: the strength of the nematic order as a function of strain. The theoretical curves are fits to Eq. (2.1) as discussed in the text. The data points are from simulation and dashed lines are the fitted results. Bottom left: Fits to the elastic energy from Eq. (2.1) for hydrostatic deformation as a function of strain. The stress is the derivative of this curve. The data points are from simulation and dashed lines are the fitted results. Bottom right: the fitted parameter C as a function of p.
tension for g and pure shear for the rest, separately via least-square fitting. The comparison between simulation and the fit are shown in Figure 2. As mentioned above, strain-stiffening depends strongly on p. This is captured by the parameter C, which controls where the potential becomes strongly nonlinear and increases the shear modulus. The fitted value of C as a function of p is shown in Fig. 2, bottom right. Quite remarkably, and in accord with our interpretation, there is a sharp peak at pc , showing that γ ∗ vanishes as p → p− c . We used the fitted parameters to calculate the stressstrain and nematic-strain curves for simple extension. In Fig. 3 we compare to simulation. The agreement is fairly good, except for q at large alignments; c.f. the discussion below. IV.
CELLS IN ECM
With our results for homogeneous deformations in hand, we turn to a schematic representation of a cell in the ECM. We can make a very simplified picture by cutting a circular hole in the material at some point, and applying a negative pressure to the exposed surface to represent the contractile force of the cell on the material [11]. This is done in our Landau theory by applying
4 κ =10−4 ,Simple Extension Simulation: p=0.80 Prediction from Theory: p=0.80 Simulation: p=0.60 Prediction from Theory: p=0.60
103 10
2
101
0.8 0.7
Alignment order q
0.6
Energy
100
0.5
10-1 10-2
0.4
10-3
0.3
10-4 10-5
0.2 0.1
10-6 10
κ =10−4 ,Simple Extension Simulation: p=0.80 Prediction from Theory: p=0.80 Simulation: p=0.60 Prediction from Theory: p=0.60
-7
10-3
10-2
10-1
Strain
0.0 0.0
100
0.2
0.4
(a)
Strain
0.6
0.8
1.0
(b)
FIG. 3. Comparison of the simulation results (data points) and theory predictions (dashed lines) for the simple extension case at p = 0.6, 0.8. The theory uses parameters fitted from the pure shear and the hydrostatic extension cases.
displacement u
0.12 0.10
Prediction Far field solution Simulation
0.20 Alignment order q
Prediction Far field solution Simulation
0.10
0.80
0.05
0.60
0.02
0.40 4
6
8
10
Distance from center node r
12 14
6
8
10
Distance from center node r
(a)
12
14
(b)
FIG. 4. Comparison of the simulation results (data points) and theory predictions (solid lines) for a model for a cell in ECM. The theory uses the parameters derived for the homogeneous deformation case.
boundary conditions. For the lattice model, we apply forces to nodes of the lattice at radius of several lattice spacings and relax the lattice. To use Eq. (2.1) in this case, we note that we have two variables which depend on the distance, r, from the center of the hole: they are the radial component of the displacement, u(r), and the strength of the radial alignment, q(r). If we minimize F with respect to u, q we find two coupled equations: 2u′ 2t 2u − 2 = r r λ + 4µ 2t ′ u u + = V ′ (q). 3 r
u′′ +
q′ q + 3 r
Alignment localization of this type is captured by our theory through the nonlinear potential V . For small r, q is large and the nonlinear potential leads to a weaker dependence of q(r) on r; that is, alignment is more or less constant near the cell [2]. We used the parameters from the homogeneous deformation case to calculate u(r), q(r) by numerically solving Eq. (4.1); see Fig. 4. We also studied the same problem using simulation by choosing a node and pulling in on nodes a few lattice spacings away. The two kinds of results agree, as shown in Fig. 4. V.
SUMMARY AND EXTENSIONS
In summary, we have given a Landau-type theory for elastic non-linearity in biopolymer gels. The central focus is to introduce the alignment tensor as a part of the order parameter. This corresponds to the notion that for many biopolymers, particularly collagen-I, the non-linear behavior is due to switching of deformation from bending to stretching as the fibers become parallel. It also should, eventually, allow an attempt to treat multi-cell systems using Eq. (1.2) as a starting point. Of course, this calculation is only the beginning of a treatment for cells. We should probably represent the cell as a force dipole rather than a spherical object: moving cells are polarized. And we need to put more than one cell in the system to see the interactions. We can suggest a more satisfying alternative to the intrinsic non-linearity in the system represented by g. We prefer it because our point of view is that all of the complex behavior results from changes in alignment. So far, Eq. (2.1) has no coupling to volume deformation – i.e. only the traceless part of v induces alignment. However, even for hydrostatic expansion, there will be local alignment and weak parts of the disordered system will rotate. We propose to treat this by introducing a random tensor field, h(r) with hhi = 0 which accounts for disorder; see also [13]. In the lattice model it is related to the missing bonds. Then we can add a term to the integrand in Eq. (2.1) of the form −Tr vTr (h · Q). In the same way as above, in the linear regime this renormalizes the first Lam´e coefficient (and hence the bulk modulus, K): λ → λ − Trhh2 i/A.
(4.1)
For large r the deformation is small, we expect the system to be linear. If we use V ≈ Aq 2 /2 we recover the familiar result for the deformation in two dimensions [12] u(r) ∼ 1/r. The second equation of (4.1) shows that q ∼ 1/r2 . As r decreases the deformations get larger, and we enter the nonlinear elasticity regime. As discussed in Ref. [2] there will be a critical radius within which nonlinear effects are important, and alignment is large.
As in the previous case the renormalization decreases in the non-linear regime. The leading order effect is the same as the term in g that we added by hand. ACKNOWLEDGMENTS
We would like to thank F. MacKintosh, C. Broedersz, M. Das, H. Levine, and E. Ben-Jacob for useful conversations. J-C Feng is supported by the National Science Foundation Center for Theoretical Biological Physics (Grant PHY-1308264).
5
[1] C. Storm, J. Pastore, F. MacKintosh, T. Lubensky, and P. Janmey, Nature 435, 191 (2005). [2] L. M. Sander, Journal of biomechanical engineering 135, 71006 (2013). [3] P. A. Agudelo-Garcia, J. K. De Jesus, S. P. Williams, M. O. Nowicki, E. A. Chiocca, S. Liyanarachchi, P.-K. Li, J. J. Lannutti, J. K. Johnson, and S. E. Lawler, Neoplasia (New York, NY) 13, 831 (2011). [4] D. Vader, A. Kabla, D. Weitz, and L. Mahadevan, PloS one 4 (2009). [5] P. Onck, T. Koeman, T. van Dillen, and E. van der Giessen, Physical Review Letters 95, 178102 (2005). [6] A. M. Stein, D. A. Vader, D. A. Weitz, and L. M. Sander, Complexity 16, 22 (2011). [7] T. Lubensky, R. Mukhopadhyay, L. Radzihovsky, and
X. Xing, Physical Review E 66, 011702 (2002). [8] M. Das, F. C. MacKintosh, and A. J. Levine, Phys. Rev. Lett. 99, 038101 (2007). [9] C. P. Broedersz, X. Mao, T. C. Lubensky, and F. C. MacKintosh, Nature Physics 7, 983 (2011). [10] M. Wyart, H. Liang, A. Kabla, and L. Mahadevan, Phys. Rev. Lett. 101, 215501 (2008). [11] M. Buenemann, H. Levine, W.-J. Rappel, and L. M. Sander, Biophysj 99, 50 (2010). [12] L. D. Landau, E. M. Lifshitz, A. M. Kosevich, and L. P. Pitaevskii, Theory of elasticity (Pergamon Press, Oxford; New York, 1986). [13] M. Sheinman, C. P. Broedersz, and F. C. MacKintosh, Phys. Rev. E 85, 021801 (2012).
