Nonlinear Mechanics of Athermal Branched Biopolymer Networks

Feb 22, 2016 - force interactions simply because the formation of a connected network ..... system size approaching infinity, the nonaffine fluctuatio...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCB

Nonlinear Mechanics of Athermal Branched Biopolymer Networks R. Rens,†,‡ M. Vahabi,† A. J. Licup,† F. C. MacKintosh,*,† and A. Sharma† †

Department of Physics and Astronomy, Vrije Universiteit, Amsterdam 1081 HV, The Netherlands Institute of Physics, University of Amsterdam, Amsterdam 1098 XH, The Netherlands



ABSTRACT: Naturally occurring biopolymers such as collagen and actin form branched fibrous networks. The average connectivity in branched networks is generally below the isostatic threshold at which central force interactions marginally stabilize the network. In the submarginal regime, for connectivity below this threshold, such networks are unstable toward small deformations unless stabilized by additional interactions such as bending. Here we perform a numerical study on the elastic behavior of such networks. We show that the nonlinear mechanics of branched networks is qualitatively similar to that of filamentous networks with freely hinged cross-links. In agreement with a recent theoretical study,1 we find that branched networks also exhibit nonlinear mechanics consistent with athermal critical phenomena controlled by strain. We obtain the critical exponents capturing the nonlinear elastic behavior near the critical point by performing scaling analysis of the stiffening curves. We find that the exponents evolve with the connectivity in the network. We show that the nonlinear mechanics of disordered networks, independent of the detailed microstructure, can be characterized by a strain-driven second-order phase transition, and that the primary quantitative differences among different architectures are in the critical exponents describing the transition.

1. INTRODUCTION Disordered fibrous networks are ubiquitous in biology. Fibers composed of naturally occurring biopolymers interconnect with each other, forming an elastic network. The network architecture depends on the particular biopolymer and polymerization conditions. One class of networks consists of individual fibers, interconnected by small cross-linking proteins. Another architecture, often seen in biology, is that of branched structures in which fibers split into two or more strands that can merge with other strands. When these two processes of splitting and merging occur repeatedly, a branched elastic network is formed. Such networks appear naturally within, as well as around, animal and plant cells.2 Actin, a major component of the intracellular cytoskeleton, is known to form branched networks near the cell membrane.3−5 In loadbearing tissues, the major component of the extracellular matrix is collagen, which forms branched networks.6−10 An important difference between intracellular filaments (actin and intermediate filaments) and collagen fibers is that the latter are rather thick and have a longer persistence length of order ∼1 cm.11−14 Models based on stretching out thermal fluctuations can account quantitatively for the stiffening response of intracellular networks.15−23 Unlike intracellular filaments, however, thermal fluctuations are unlikely to play an important role in the elasticity of collagen networks.9 Fiber networks with controlled connectivity based on 2D triangular and 3D FCC based lattices1,24−31 as well as Mikado structures32−38 have been used as models for collagen networks. These models can account for several features of the nonlinear mechanics of collagen networks; however, it still remains unclear whether taking branching explicitly into © XXXX American Chemical Society

account has any significant impact on the reported mechanical behavior. Because of the relatively large persistence length of collagen, it is justified to study the mechanics of a collagen network using a athermal model.13,14 A branched disordered network with only central force interactions is, however, unstable.39 Only when the connectivity or coordination number z exceeds the threshold of 2d, where d is the spatial dimensionality, are central force interactions sufficient to stabilize a network. Below this threshold of marginal stability, such networks are have a vanishing linear elastic modulus. Most of the naturally occurring networks are submarginal with respect to centralforce interactions simply because the formation of a connected network occurs either by cross-linking of two different fibers (z = 4) or by branching (z = 3) of fibers. Thus, networks with branching are submarginal in both 2D and 3D. Collagen networks, however, are stable to small deformations. In fact, the linear elastic modulus is known to show quadratic scaling with the concentration of collagen.1,24,40 It follows that any model of branched networks must take into account additional interactions in a submarginal network to account for a finite linear modulus. Although there are several fields such as normal stresses, active stresses,41 and thermal fluctuations,42 the most relevant field in the case of collagen is stabilization by bending interaction that arise from the Special Issue: William M. Gelbart Festschrift Received: January 9, 2016 Revised: February 18, 2016

A

DOI: 10.1021/acs.jpcb.6b00259 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

cost associated with rotating filaments relative to each other; however, in branched biopolymer networks such as collagen, the fibers themselves consist of multiple fibrils, which are covalently cross-linked internally. The point at which two fibers seem to connect is actually where the collection of fibrils splits up and combines.50 High-resolution images of collagen reported in ref 51 show clear evidence of this feature. We refer to the splitting as a branching point, which has 3-fold local coordination number (z = 3) by construction. Because the local connectivity cannot exceed three, the network is strictly submarginal.26,39 A schematic of a branching point is shown in Figure 1b. In our description of the structure of branched

persistence of the filaments. This is suggested by the observed quadratic scaling of modulus with concentration, corresponding to a linear scaling of modulus with the bending rigidity,1,24,43,44 as found in several computational studies on disordered networks.1,24,27,32−36,45 Here we present numerical and analytical study of branched networks. With numerical methods we obtain the mechanical response of the branched networks to an externally imposed shear-deformation. Our analytical approach is based on a scaling analysis recently applied in ref 1, where it was shown that the nonlinear mechanics of athermal fibrous networks can be understood as a strain-driven critical phenomenon. It was demonstrated that such critical behavior is a generic feature of submarginal networks. With an approximate but highly accurate scaling relation including the two critical exponents f and ϕ, ref 1 was able to predict the nonlinear stiffening curves of collagen over a wide range of concentrations; however, the primary focus of that study was fibrous networks with freely hinged cross-links and average connectivity strictly above three. In the present study, the average connectivity of all of the networks considered is below three: We consider an extreme limit in which branching is the only mechanism of network assembly. We demonstrate that the same scaling analysis as in ref 1, with different critical exponents f and ϕ, can fully capture the nonlinear mechanics of branched networks considered in this study. In this study, all simulation results correspond to branched networks in 2D. Recently, it has been shown46 that there is no fundamental difference between the nonlinear mechanics of fibrous networks in 2D and 3D, provided that the average connectivity is below 4, the threshold of marginal stability in 2D. In fact, by scaling the computationally obtained modulus with the line density, the nonlinear mechanics in 2D and 3D show good quantitative agreement. We believe, therefore, that our results are equally applicable to branched networks in 3D. The rest of the paper is organized as follows. In Section 2 we describe our computational model, which takes into account the local branched structure. We then describe the different geometric variations of branched networks considered in this study. In Section 3 we present our results on the elastic behavior of branched networks over the entire elastic regime. This section includes the performed scaling analysis of stiffness curves obtained from different network geometries. In Section 4 we discuss several aspects of criticality in branched networks. In Section 5 the role of disorder on the criticality is explored by comparing distorted with undistorted undiluted honeycomb lattices. In Section 6 we show that the critical exponents are not constant but depend on the average connectivity of a network. Finally, in Section 7 we discuss our findings and draw conclusions.

Figure 1. (a) Underlying structure of collagen is thought of as multiple long running fibrils internally cross-linked to form a single bundle that occasionally branches out. At a branching point individual fibrils move in different directions forming two separate bundles. The fibrils, which belong to two different bundles, are bent and hence store bending energy. (b) Each individual bundle is connected between two network nodes. At each branching point there are three bond pairs that can potentially contribute to the bending energy; however, we disregard the possibility of bundle splitting due to highly bent fibers. This is shown above as the pair made of 0−2 and 0−3. In the model this is taken into account by assigning an individual bending rigidity to each pair of bonds, with the bending rigidity of the smallest initial angle set to 0.

networks we ignore the possibility of two bundles cross-linked to each other, thus ensuring that the connectivity cannot exceed three. This allows us to calculate the bending energy stored in the network by considering only the branching points. We adopt the stretching part of the Hamiltonian from the rigid rod model15,23,52 because the stretching behavior of most biopolymers is not expected to be different from rod-like filaments. The stretching contribution to the Hamiltonian for an individual bond in the network is defined as /stretch =

μ 2

∑ i,j

1 (|X⃗i − X⃗j | − l0, ij)2 l0, ij

(1)

where i and j run over all connected nodes. The distance between two connected nodes can be expressed in terms of the difference in position vectors X⃗ i − X⃗ j. Here μ is the effective stretch modulus for multiple cross-linked fibers that act together as one filament and is taken to be the same for all elements. The rest length of a bond that connects two nodes i and j is denoted by l0,ij, which is defined as the end-to-end distance between the nodes at zero strain. As can be seen from the quadratic form, the effect of both compression and extension of bonds with respect to this rest length is symmetric as in a linear Hookean spring. Because the underlying structure of the network consists of long running branched fibers, the Hamiltonian for the bending part should capture how the fibers share bending energy at the branching points. The bending energy of the network is given by

2. MODEL Our focus in this study is branched networks. Many systems are known to consist of an internal branched structure.47−49 Here we use collagen to motivate the network model as follows. The underlying structure of collagen is different from other filaments found in biological systems such as the intracellular filaments, microtubules, and F-actin. Intracellular filaments are relatively straight and are able to cross-link by specific proteins and ultimately form a network structure. Physically, the effect of a cross-link between two filaments is a point-constraint forcing the two points of attachment at the cross-link to move together. The relative orientation is irrelevant; that is, there is no energy B

DOI: 10.1021/acs.jpcb.6b00259 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B /bend =

∑ ⟨ijk⟩

κijk (l0, ij + l0, jk)

order. A selected node gets its bonds removed in a random direction until a local connectivity of maximal three remains. This process is repeated for all nodes. An inherent feature of this protocol is that the connectivity can drop below three for a small fraction of the nodes. After obtaining a network with a connectivity very close to three, we can dilute such a network again by removing bonds in the system with a probability 1 − pbond. Note that the dilution is with respect to a disordered reference state and not with a full triangular network. These two geometries, although both originating in highly regular lattices, are structurally very different. In particular, the difference is seen in the initial angle distribution as well as the distribution of length between every consecutive branching. Any difference in the nonlinear mechanics of these two differences must therefore originate in the underlying structure. Also, the influence of the regularity of lattice-based networks on mechanics is studied. For this, we consider distorted versions of the two lattices previously described, as shown in Figure 2c,d. Lattices are distorted by simply displacing every node in the network over a distance of δ/2 in a random direction, where δ is a randomly chosen number between 0 and dmax. The maximum displacement, dmax, can be varied to control the amount of distortion in the network and can reach a maximum of l0. The main impact of distortion is on the initial angle and bond length distribution. In all cases, periodic boundary conditions are imposed in all directions. At every branching point we have the freedom to assign a different bending rigidity to each adjacent pair of bonds as shown in Figure 1. Besides considering the case where a finite bending rigidity is assigned to each of the three pairs of bonds, we consider other scenarios as well. Consider a collection of bundles which at one point branches into two separated bundles as shown in Figure 1a, then one can regard the pair of bonds 0−2 and 0−3 as not sharing any bending energy. In the model we include this by setting the associated bending rigidity κ203 of this pair of bonds as shown in Figure 1b to zero. Since in the model there is no notion of an underlying directionality, we include the branching by setting the bending rigidity corresponding to the smallest of the three initial angles to zero. Physically it is intuitive to think of this as the branching direction. To rule out the effect of this arbitrary choice, we have considered other variants of this scheme, for example, randomly setting one of the three bending rigidities to zero. From these variations, we find that the choice of zero bending rigidity bond-pair does not change the qualitative features of nonlinear mechanics of branched networks. In this article we present the data for network in which κ̃ for the smallest angle at a branching point is set to zero. In this study, we focus on branched networks in 2D. The shear response of the network is determined in the following way. For a given imposed strain, γ, the network is subjected to a simple shear deformation by displacing each node in the network according to

(θijk − θ0, ijk)2 (2)

where ⟨ijk⟩ refers to a pair of connected bonds in the network. For any pair there is a bending rigidity, κijk, which indicates the magnitude of resistance to bending deformation. At a branching point we take κijk to be κ, for the smallest angle we put the bending rigidity to 0 (as indicated in Figure 1b). From the Hamiltonian, it follows that the bending force is modeled as a linear response to the change in angle, θijk, with respect to the initial branching angle θ0,ijk. The two elastic parameters in our model are the stretching rigidity, μ, and bending rigidity, κ. In simulations, rather than specifying both parameters, we vary the reduced bending κ rigidity defined as κ̃ ≡ μl 2 . This dimensionless parameter is a 0

measure of the relative strength of bending to stretching and is in fact the relevant parameter for comparison with experiments as it scales linearly with the concentration in experiments.1,24,46 All numerical simulations are performed by setting μ and l0 to 1. We report κ̃ used in the simulations. A honeycomb lattice is well-suited to represent a branched network. An undiluted lattice has a coordination number three at every node in the lattice and is therefore submarginal. Previous reports on collagen networks have reported a coordination number close to three.53 We perform random dilution of a full honeycomb lattice with a probability 1 − pbond, where pbond is the probability for a bond to exist. Dilution introduces structural disorder and brings the average coordination number below three. In the experimental study of ref 51, it can be seen that collagen networks show structural similarity with diluted honeycomb lattices. A schematic of a diluted honeycomb lattice is shown in Figure 2a. To determine the significance of branching, we also consider a triangular lattice, as shown in Figure 2b, which has been diluted in such a way to have a local connectivity never exceeding three. In the protocol of obtaining such a network we start with a full triangular network (connectivity 6) and visit nodes in random

⎛1 γ ⎞ X⃗aff = ⎜ ⎟X⃗ init ⎝0 1⎠

(3)

Here X⃗ init denotes the initial position vector of a node in the network and X⃗ aff denotes the positions after the affine deformation. We then minimize the elastic energy stored in the network using the conjugate gradient algorithm.54 The final configuration of the network nodes X⃗ deviates from the initially imposed affine deformation. These local rearrangements, X⃗ nonaff

Figure 2. Schematic of a typical network obtained after dilution and distortion of a lattice-based network. (a) Honeycomb lattice after dilution with bond probability pbond = 0.85 and (c) a triangular lattice with pbond = 0.95. (b,d) Corresponding distorted lattices. The applied distortion is dmax = 0.8. After applying distortion the obtained configuration is set to be the relaxed state of the network. C

DOI: 10.1021/acs.jpcb.6b00259 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 3. Schematic of nonaffine displacements. (a) In black, the initial network configuration. (b) In black, the network after an affine shear deformation. The affinely displaced network is not in force balance and hence the nodes undergo further displacement to achieve force balance. In red, the final configuration of the network including local nonaffine displacements. Green arrows denote the nonaffine displacements.

Figure 4. Upper panels: Differential shear modulus K versus shear-strain, γ, for (a) an undistorted honeycomb lattice with a connectivity of ⟨z⟩ = 2.68 and (b) an undistorted triangular lattice with ⟨z⟩ = 2.87. The solutions to eq 5 are plotted with dotted lines. In panel a the solid curve indicates the onset of stiffening γ0. The dashed curve indicates the strain at which critical transition occurs. The inset shows scaling of the linear modulus, G0, with the bending rigidity, κ̃; the solid line indicates a unit slope. Lower panels: K rescaled for different values of κ according to the scaling relation given by eq 4. The black curves in panels c and d correspond to the data of panels a and b, respectively. The blue curves in these panels correspond to the distorted version of the lattices (c) dmax = 0.8 and (d) dmax = 0.5 and are arbitrarily shifted vertically from the data for undistorted networks.

= X⃗ − X⃗ aff, are referred to as the nonaffine displacements, as shown schematically in Figure 3. In Section 4, we show that with increasing strain, nonaffine displacements in the network increase dramatically. We show that such rapid increase is a signature of the highly nonlinear mechanics of the network. By applying Lees−Edwards boundary condition to the system a shear-deformation is imposed on the length scale of the system size.55 Once we obtain the energy stored in the network as a

function of γ we calculate the shear modulus as K = the shear stress as σ =

1 ∂/ , A ∂γ

1 ∂ 2/ A ∂γ 2

and

where / = /stretch + /bend and

A is the total area that the 2D network occupies. As previously mentioned, the two simulation parameters are the reduced bending rigidity, κ̃ = κ/μl20, and the bond probability, pbond. The total energy density in the network is jointly determined by κ̃ and pbond through the line density, ρ, defined as the total length D

DOI: 10.1021/acs.jpcb.6b00259 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B ⎛ κ̃ ⎞ ⎟ K ≈ |Δγ | f .±⎜ ⎝ |Δγ |ϕ ⎠

per unit area of the network. The computationally obtained modulus and stress are therefore measured in units of μρ1,24,46 throughout this paper.

(4)

Here Δγ = γ − γc is the strain interval measured relative to the critical point γc. This scaling is analogous to that for the conductivity of random resistor networks and fiber networks as a function of connectivity. The crossover function .±(x) captures the functional dependence with the positive (+) branch corresponding to Δγ > 0 and the negative (−) branch to Δγ < 0. For γ ≪ γc, .−(x) approaches x and thus K ≈ κ̃. On the contrary, for γ ≫ γc, .+(x) approaches a constant value and thus K becomes independent of κ̃. In Figure 4c this function is used to collapse the stiffening curves of Figure 4a. In Figure 4d, the collapse corresponds to the stiffening curves presented in Figure 4. The two sets of collapse in each panel of Figure 4c,d correspond to the undistorted (black) and distorted (blue) versions of the network. To obtain a good collapse of the data, the value for the exponents f and ϕ were selected. As can be seen, the collapse consists of three branches. The linear regime corresponds to the lower branch with K ≈ κ̃. The stretching regime, independent of κ̃ has zero slope and corresponds to the upper horizontal branch. These two branches merge into a single branch which captures the approach toward the critical point, that is, Δγ → 0. The slope of this branch, f/ϕ, determines the functional dependence of modulus as K ≈ κf/ϕμ1−f/ϕ at the critical point, similar to refs 26, 28, and 41. In ref 1, it was suggested that the critical exponents depend on the average connectivity of the network. The two networks considered here, honeycomb and triangular lattice, have quite similar average connectivity and indeed we see that their critical exponents do not differ much. However, as we will show later in this study, the exponents evolve with the average connectivity in the network. It follows from the above that the nonlinear mechanics of branched networks can be quantitatively captured in exactly the same fashion as for the filamentous networks. The crossover function .±(x) describes the shear modulus of the network over the entire range of strain. Within the framework of critical phenomenon, an equation can be derived for the crossover function. This was shown to accurately describe the corresponding crossover phenomenon in ferromagnetic phase transition.56 Because the strain-driven criticality in networks is analogous to ferromagnetic phase transition, a corresponding equation for the crossover function .±(x) can be derived.1 The equation is given as

3. NONLINEAR MECHANICS OF BRANCHED NETWORKS In Figure 4 we show the shear modulus K versus shear strain curves for (a) honeycomb and (b) triangular lattice with an average coordination number of ⟨z⟩ = 2.68 and 2.87, respectively. These stiffening curves are obtained from simulations of undistorted (dmax = 0.0) networks. We first describe the qualitative features of a stiffening curve. Each of the stiffening curves can be characterized by an initial linear regime in which the elastic modulus does not vary with the strain γ. The linear modulus G0 scales with κ̃,43,44 as shown in the inset of Figure 4a. In this regime, the dominant mode of local deformations within networks is bending. Beyond a certain strain, denoted as γ0, the modulus increases rapidly with strain before it converges to a constant value at large strains. The strain γ0 indicates the onset of nonlinear stiffening, and here we arbitrarily define it as the strain at which the nonlinear modulus has become approximately three times the linear modulus. As can be seen, γ0 is practically independent of the bending rigidity. The independence of γ0 on the bending rigidity is in agreement with the recent work of Licup et al.,24 in which it was demonstrated that such behavior can be understood by using simple geometric arguments as opposed to energetics and is observed in experiments on collagen as well.1,24 We note that the stiffening curves for distorted networks (honeycomb and triangular) exhibit the same qualitative features. (See Figure 4c,d.) As previously mentioned, γ0 marks the onset of the nonlinear regime apparent in the rapid increase in the modulus with increasing strain. The nonlinear regime can be divided into two subregimes. The first of which extends from γ = γ0 to γ = γc, as shown in Figure 4a, and the second from γ = γc and above. As will be shown later, the strain γc marks the onset of rigidity in a submarginal network with only central force interactions. Experimentally, the regime for large strains, γ ≫ γc, where the modulus approaches a constant value, is elusive primarily due to the failure of material;24 however, as can be seen in Figure 4, it is in the first subregime where most of the nonlinear stiffening occurs. In the second subregime, the differential modulus increases much slower, approaching a constant value for large strains. In fact, γc can be identified as the inflection point of the stiffening curves for κ̃ → 0, suggesting that γc marks the strain at which a network enters into a new regime (subregime). We can quantitatively capture the nonlinear mechanics of branched networks over the entire strain range in the following way. Reference 1. recently showed that the straindriven transition of a submarginal network (with only central force interactions) from floppy to rigid states is analogous to a second-order phase transition. The relevant order parameter is the modulus K, which increases continuously from zero with strain increasing beyond the critical strain γc. With the reduced bending rigidity κ̃ acting as an auxiliary field, a scaling law was proposed that captures the crossover from the bend-dominated linear regime to the stretch-dominated high-strain regime of filamentous networks. We perform the same scaling analysis on the stiffening curves obtained from simulations on branched networks. The scaling function motivated in ref 1 is given by

(ϕ− f ) K ⎛ K1/ f ⎞ κ̃ ⎟ ⎜±1 + ≈ |Δγ | ⎠ |Δγ | f ⎝ |Δγ |ϕ

(5)

where the ± solutions correspond to the two branches of the collapse described by the corresponding scaling functions G±. We use this to analytically calculate the modulus for any given strain. Critical exponents, f and ϕ, obtained from the collapse according to eq 4 and the critical strain, γc, are given as inputs to eq 5. In Figure 4a,b we have superimposed the predicted modulus according to eq 5 on the simulated data. The predictions of eq 5 are in excellent agreement with the simulations for both honeycomb and triangular lattice-based networks, as shown in Figure 4a,b. As previously shown the nonlinear mechanics of branched networks exhibit qualitative and quantitative similarities to the filamentous networks. As for filamentous networks,1 eq 5 can E

DOI: 10.1021/acs.jpcb.6b00259 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(v) the red bonds indicate stretching as the major component of the local deformation. Comparing the lattice at these three strains, we can see how for higher strains the bonds increase their alignment. Between these two geometrical phases a crossover occurs. As the strain approaches the critical strain, the network enters a mixed phase regime with stretching becoming increasingly important with increasing strain. In Figure 6a (iii), the network configuration is shown very close to the critical strain. At the critical strain, the stretching and bending energies become comparable to each other. This is highlighted in the upper panel of Figure 6b, in which the ratio of bending and stretching energy is plotted. We note that at the critical strain we do not necessarily have percolation of red bonds. Beyond the critical strain, stretching becomes the dominant mode of deformation in the network. As shown in the upper panel of Figure 6b, stretching energy is negligible in comparison with the bending energy for γ < γc. Both bending and stretching energies increase rapidly in the range γ0 ≤ γ ≤ γc. At γ = γc, bending and stretching energies become comparable to each other. This observation highlights the difference between stiffening strain and critical strain. The onset of stiffening γ0 does not mark the onset of stretchdominated behavior. Stretch-dominated behavior commences at γc > γ0. In the presence of bending interactions, γ0 marks the onset of the bend-to-stretch transition with the transition completing itself at γ = γc. In the range γ0 ≤ γ ≤ γc, the network deforms in a highly nonaffine manner. The strain γc coincides with the strain at which the nonaffine fluctuations in the network diverge. The nonaffine fluctuations are a consequence of large-scale rearrangements in the network. A quantitative measure of these rearrangements is the differential nonaffinity δΓ(γ) expressed in terms of the node position by32,57−59

accurately predict the modulus as a function of strain for branched networks. On the basis of these observations, using filamentous networks to understand nonlinear mechanics of collagen seems to be justified; however, it is interesting in its own right to further test criticality in branched networks.

4. SIGNATURES OF CRITICAL BEHAVIOR IN BRANCHED NETWORKS A submarginal network, with only central force interactions (κ̃ = 0), when subjected to simple shear, becomes rigid for strains |γ| > γc. The transition from floppy to rigid states, from the perspective of critical phenomenon, must be continuous with the modulus K vanishing at γ = γc. More precisely, K should exhibit a power-law dependence of the type K ≈ |Δγ|f where Δγ = γ − γc ≥ 0. We show in Figure 5 the modulus K as a function

Figure 5. Elastic modulus for a honeycomb lattice plotted against Δγ+ for three different κ̃ values (different colors) for three different system sizes (different symbols), where Δγ+ is the distance to the critical point γc for points above γc, so Δγ > 0. Solutions to eq 5 for each bending rigidity are plotted with lines. In the limiting case of κ̃ = 0, we find scaling of K ≈ |Δγ|f. This result is plotted with the solid black curve. Simulation data approach this solution for larger systems with low bending rigidity.

δ Γ(γ ) =

1 l02Ndγ 2

N

∑ (dX⃗i i

aff

→ ⎯ − dXi)2

(6)

dX⃗ aff i

where is the imposed affine displacement vector and dX⃗ i is the true displacement vector of node i in response to the differential strain dγ. δΓ(γ) is somewhat analogous to susceptibility in a ferromagnetic phase transition. For a diluted honeycomb lattice the δΓ(γ) is shown in Figure 6c. As can be seen, δΓ(γ) peaks at γ = γc. The height of the peak increases with decreasing κ̃. In fact, in the thermodynamic limit of the system size approaching infinity, the nonaffine fluctuations will diverge at γ = γc as κ̃ → 0. The increase in fluctuations at the critical strain with decreasing κ̃ is analogous to increasing susceptibility at the Curie temperature in a ferromagnetic phase transition with decreasing magnetic field.

of Δγ+ for several system sizes N and κ̃ values. The solutions to eq 5, valid in the limit of the infinite system size, for a set of κ̃ values are also shown in the Figure 5. The black solid line corresponds to κ̃ = 0. For any system size, when κ̃ > 0, the modulus is finite at γc and scales with κ̃ such that in the limit of κ̃ → 0 the modulus vanishes at γc. The power-law scaling of modulus becomes increasingly apparent for lower values of κ̃. Such scaling will continue all the way up to γc in the limit of κ̃ → 0 only in the thermodynamic limit of N → ∞; however, it is evident that for very small bending rigidity and extremely large system sizes the scaling K ≈ |Δγ|f can be observed ever closer to the critical strain. In the presence of bending interactions, the network is stable for γ < γc with its linear modulus scaling linearly with κ̃. In Figure 6, we show how a network configuration evolves with the applied strain for a bending rigidity of κ̃ = 10−6. We show how the dominant mode of deformation changes from bending to stretching with increasing strain. When the applied strain is low enough to be in the linear regime, the network is in a bending phase, as shown in Figure 6a (i) and (ii), with the invoked bending of bonds indicated in blue. At large strains (γ ≥ 100%), the network can be regarded as being in the stretching phase because in this limit almost all of the bonds are stretched out due to geometric alignment. In Figure 6a (iii) to

5. CRITICALITY IN A FULL HONEYCOMB LATTICE In the previous section, we showed that the diluted honeycomb lattice-based network exhibits a linear shear modulus governed by bending energy, followed by a rapid nonlinear stiffening beyond a threshold strain. This agrees with similar computational studies on other submarginal lattice-based networks and randomly connected Mikado networks.27,28,34−36 Lattice-based networks are appropriate to model the behavior of biological networks, provided that structural disorder is added to those lattices.1,46 Structural disorder can be introduced by diluting a full lattice, as shown in Section 2. Such bond dilution, however, introduces not only geometric disorder but also a reduction in the connectivity. The honeycomb lattice represents an F

DOI: 10.1021/acs.jpcb.6b00259 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 6. (a) Network undergoing a shear deformation. Here the bending rigidity is κ̃ = 10−6. For a range of shear-strain values, γ, ranging from 10 to 150%, the obtained network configurations are plotted. The bonds are colored according to the local stress, black in the absence of stress, blue for stress due to bending force, and red for stretching force. (i,ii) For low strains, the dominant deformation mode is bending, as shown in blue. (iii) At the critical strain, the network is in a mixed regime with both bending and stretching present. It can be seen how the bonds, belonging to a specific percolating cluster, transit to red, indicating the increasing stress close to the critical strain. Such a percolating cluster of red bonds will cease to exist on reducing κ̃ to sufficiently low values. (iv) Increasing the strain above the critical strain invokes many stretching modes throughout the network. (v) Finally, transiting to very large strains, the stretched bonds in the network align in the direction of the shear strain. In the lower panel of (b) the stiffness of the network is shown as a function of strain. The rapid increase in stiffness close to the indicated critical strain coincides with the observed geometrical transition in the network. In the upper panel of (b) the ratio of the two energy contributions is plotted to demonstrate the narrow range over which the transition between two phases occurs. Notice that the nonlinear stiffness up to the critical point should be attributed merely to bending interactions. (c) Differential nonaffinity δΓ(γ) as a function of deformation. The differential nonaffinity peaks at the critical strain γc.

interesting alternative structure because of its intrinsic submarginality in 2D, permitting the separate introduction of geometric disorder, independent of connectivity, as we demonstrate here. In a regular, full honeycomb lattice, because of the symmetry one expects that any nonaffine displacements within the network occur within a unit cell. It was shown in ref 60 that a regular honeycomb lattice is unstable to small deformations that lead to rotation of unit cells without stretching any bonds. Thus, such a network has a vanishing linear shear modulus. Interestingly, however, the regular honeycomb lattice becomes stable for any finite shear strain, making it intrinsically nonlinear, even at small strain. In the presence of bending interactions, the regular honeycomb lattice also becomes linearly stable. Moreover, as we show later, in the presence of finite spatial disorder, the full honeycomb lattice exhibits straindriven critical behavior at a finite strain threshold, analogous to the diluted honeycomb and generic submarginal networks.1 One can now specifically explore the role of spatial disorder, as opposed to dilution, in a full honeycomb lattice by distorting it while maintaining a constant connectivity of z = 3. In Figure 7 we show the elastic modulus for an undiluted honeycomb network. The blue data sets correspond to undistorted networks with varying bending stiffness, as discussed in Section 2. Because of the symmetry of the lattice, the displacement of nodes will be correlated for any given strain, and no unique strain value exists for which the correlation length diverges. We

Figure 7. Differential elastic modulus for an undiluted honeycomblattice with an average connectivity of three at each node. In blue, an undistorted fully regular lattice and in black a distorted lattice with dmax = 0.5. The regular lattice shows a vanishing modulus as κ̃ goes to zero. The dashed line indicates a slope 2 as expected scaling for this limit. For any finite κ̃ there exists a linear regime, the extent of which depends on the κ̃. As can be seen in the black data, disorder restores a finite strain-threshold γc at which the distorted network undergoes a continuous transition in rigidity. The inset shows the positive correlation between the critical strain and the amount of distortion.

can retain a 3-fold connectivity and now include structural disorder by distorting the nodes of the network. In Figure 7, in black, the elastic modulus is shown for such a distorted full G

DOI: 10.1021/acs.jpcb.6b00259 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

at ⟨z⟩ = 2.88. In Figure 9 we show the evolution of the critical exponents of both honeycomb and triangular lattice networks

honeycomb lattice. The disorder is added by displacing every node in a random direction over a random distance, as described in Section 2. As can be seen in Figure 7, when the applied distortion is nonzero, the network exhibits strainstiffening behavior analogous to the generic submarginal networks. Here, because of the structural disorder, nonaffine movements of nodes can occur independently, making it possible to relax the energy globally. Simulations were performed on networks with different levels of distortion, from which we obtain the critical strain (see inset of Figure 7). We find that the critical strain γc decreases with decreasing distortion and approaches zero as we get closer to the undistorted case. The full dependence of γc as a function of distortion will be investigated in a future study. We conclude that a submarginal network with any finite amount of disorder, either spatial by distortion or topological by dilution, will exhibit strain-driven critical behavior. The critical strain can be varied either by disorder in the network node positions or by varying the connectivity. Moreover, the critical strain shrinks continuously to zero as the disorder is decreased to zero, corresponding to the perfect lattice. This property could be used to design networks with a desired critical strain without changing the network connectivity.

Figure 9. (a) Values of f and ϕ/2 are plotted separately to show how the critical exponents are dependent on network connectivity ⟨z⟩. Major variation was found in exponent f rather than ϕ. (b) Ratio f/ϕ, the two critical exponents, plotted for different geometries as a function of the average number of connections z per node. A strong correlation between this ratio of exponents and the average connectivity ⟨z⟩ is found. Different network geometries show comparable exponents as long as the network connectivity is similar.

6. CONNECTIVITY DEPENDENCE OF THE CRITICAL EXPONENTS Here we study the evolution of the exponents f and ϕ obtained from eq 4 with the average connectivity in a network. To that end, we create both honeycomb and triangular networks with varying average connectivity in the range 2 < ⟨z⟩ < 3. The minimum connectivity close to 2 is a constraint due to connectivity percolation of the lattices we use. For each network at a given connectivity, we perform scaling analysis as in Section 3 on the stiffening curves to obtain the critical exponents f and ϕ. In Figure 8, we show collapse of the

with the average connectivity. We identify two important features. First, the exponents obtained from both honeycomb and triangular are practically indistinguishable, implying that the exponents are predominantly determined by the average connectivity and are independent of the detailed microstructure of the network. Second, the exponent ϕ shows a very weak dependence on the average connectivity. The exponent f varies from 0.45 to 0.82 over the entire range of connectivity considered here. Similar variation in critical exponents was reported for fiber networks in ref 1. The variation of the critical exponents with connectivity z that we observe is unexpected but not completely without precedent. As noted in ref 1, the variation in exponents occurs along a line of second-order transitions in the z−γ plane. This is similar to the behavior of certain 2D models, such as the Ashkin−Teller and 8-vertex models, which each exhibit continuously varying critical exponents along a critical line.61−63 Experimentally, such a variation in critical exponents with composition is also reported for certain quantum phase transitions.64,65 Interestingly, like quantum phase transitions, our critical behavior also corresponds to an athermal phase transition; however, in the presence of disorder, as is introduced by our dilution process, it is also possible that the exponents we observe are effective, corresponding to a crossover between critical exponents in the pure and disordered limits. The apparent variation in critical exponents reported for diluted Ising models has been attributed to such crossover behavior.66,67 Whether the varying exponents we observe represent a variation of critical exponents or are the result of crossover behavior is unclear. Interestingly, we observe an apparent evolution toward mean-field values of f = 1 and f/ϕ = 1/2 as we approach the regular honeycomb lattice.26 Nevertheless, the non-mean-field exponents found for diluted

Figure 8. For two different network connectivities the collapse is plotted with freely fitted critical exponents f and ϕ. The two systems are both undistorted honeycomb lattices, but have different average connectivity due to different amount of dilution of the network. The example curves illustrate how the ratio between the critical exponents f and ϕ depends on the connectivity.

stiffening curves obtained from undistorted honeycomb lattice for two different connectivities of ⟨z⟩ = 2.46 and 2.88. As can be seen, the critical exponents are different for the two connectivities. We note that the exponent ϕ for the two connectivities is practically the same. The value of the exponent f, however, increases significantly from 0.47 at ⟨z⟩ = 2.46 to 0.8 H

DOI: 10.1021/acs.jpcb.6b00259 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

critical exponents, f and ϕ. We found that the critical exponents, as expected, do not depend on the microstructure of the branched network. Both honeycomb and triangular lattice networks have the same critical exponents for a given average connectivity. The exponents, however, apparently evolve with the average connectivity. A similar evolution of exponents with connectivity has been reported in ref 1 for filamentous networks; however, the critical exponents were not found to show a significant difference between 2D and 3D systems. Interestingly, even the critical strain threshold did not show a dependence on dimensionality for networks sufficiently far below the 2D isostatic point. This suggests that, at least for transitions induced by finite strains on the order 10% or more, the network dimensionality becomes unimportant and the local connectivity or network topology dominates.

subisostatic networks are consistent with experiments on collagen networks.1

7. DISCUSSION AND OUTLOOK We studied the nonlinear mechanics of athermal branched networks. Branched networks are ubiquitous in biology. Collagen, the major component of extracellular matrix, forms branched network and can be considered to be relatively insensitive to thermal fluctuations under physiological conditions. Motivated by this, we considered two different microstructures: honeycomb and triangular, to represent athermal branched networks in 2D. Filamentous networks have been used to model the nonlinear mechanics of collagen networks;1,13,24,28,31 however, it remained unclear whether taking branching explicitly into account has any significant impact on the reported mechanical behavior. Reference 1 recently demonstrated that the nonlinear mechanics of fibrous networks can be understood within the framework of athermal critical phenomenon. The strain-driven transition of a submarginal network, with no bending interactions, from floppy to rigid states is a continuous second-order phase transition. Once the continuous nature of transition is established, finite bending rigidity can be included as an auxiliary field. Scaling analysis from the theory of critical phenomenon yields a scaling law that can quantitatively capture the entire nonlinear mechanics of the network for any bending rigidity and strain. Because critical behavior is rather insensitive to the microscopic details of interactions within network elements, we expect the elastic behavior of branched networks to be qualitatively similar to that of the filamentous networks. Indeed, our study on branched networks demonstrates the equivalence between the two systems. To compare computationally obtained modulus with the experiments, one must relate the two parameters, κ̃ and ρ, to the experimentally controlled variables. The most relevant experimental variable is the protein concentration. It has been shown that under the assumptions of an athermal network composed of beam-like linear elements, both ρ and κ̃ scale linearly with the volume fraction of the network.1,24,46 In submarginal networks a linear regime of constant stiffness can exist due to nonaffine rearrangements that couple to the bending rigidity. On application of sufficient strain, such networks undergo rapid strain stiffening. In a symmetric latticebased network, such as the undiluted honeycomb lattice, such nonaffine rearrangements are prohibited by the regularity and periodicity of the lattice. We showed that in this system indeed no linear regime is found in the absence of fiber bending rigidity. Furthermore, we showed that it is possible to restore a finite nonzero stiffening strain threshold by introducing spatial disorder into the regular lattice while maintaining a constant connectivity z = 3. In fact, the amount by which the lattice is distorted governs the magnitude of the critical strain. Nonlinear mechanics of an undiluted honeycomb lattice-based network will be explored in detail in a future study; however, the present study reveals the importance of spatial or geometric disorder in addition to the submarginal requirement. Except in simulations, because it is practically impossible to have a lattice-based network with perfect regularity and periodicity, even undiluted regular lattice-based networks are expected to show critical mechanical behavior. The scaling law, governing the crossover from benddominated to stretch-dominated stiffness, is a function of the distance to the critical strain, the bending rigidity, and the two



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was carried out on the Dutch national einfrastructure with the support of SURF Cooperative. This work is part of the research programme of the Foundation for Fundamental Research on Matter (FOM), which is financially supported by The Netherlands Organisation for Scientific Research (NWO).



REFERENCES

(1) Sharma, A.; Licup, A. J.; Jansen, K. A.; Rens, R.; Sheinman, M.; Koenderink, G. H.; MacKintosh, F. C. Strain-controlled Criticality Governs the Nonlinear Mechanics of Fibre Networks. Nat. Phys. 2016, DOI: 10.1038/nphys3628. (2) Alberts, B.; Johnson, A.; Lewis, J.; Raff, M.; Roberts, K.; Walter, P. Molecular Biology of the Cell; Taylor & Francis Group, 2007. (3) Fletcher, D. A.; Mullins, R. D. Cell Mechanics and the Cytoskeleton. Nature 2010, 463, 485−492. (4) Pollard, T. D.; Cooper, J. A. Actin, a Central Player in Cell Shape and Movement. Science 2009, 326, 1208−1212. (5) Mullins, R. D.; Heuser, J. A.; Pollard, T. D. The Interaction of Arp2/3 Complex with Actin: Nucleation, High Affinity Pointed End Capping, and Formation of Branching Networks of Filaments. Proc. Natl. Acad. Sci. U. S. A. 1998, 95, 6181−6186. (6) Comper, W. D. Extracellular Matrix; CRC Press, 1996; Vol. 1. (7) Hay, E. D. Cell Biology of Extracellular Matrix; Springer Science & Business Media, 1991. (8) Keene, D. R.; Engvall, E.; Glanville, R. W. Ultrastructure of type VI Collagen in Human Skin and Cartilage Suggests an Anchoring Function for this Filamentous Network. J. Cell Biol. 1988, 107, 1995− 2006. (9) Motte, S.; Kaufman, L. J. Strain Stiffening in Collagen I Networks. Biopolymers 2013, 99, 35−46. (10) Piechocka, I. K.; van Oosten, A. S.; Breuls, R. G.; Koenderink, G. H. Rheology of Heterotypic Collagen Networks. Biomacromolecules 2011, 12, 2797−2805. (11) Lindstrom, S. B.; Kulachenko, A.; Jawerth, L. M.; Vader, D. A. Finite-strain, Finite-size Mechanics of Rigidly Cross-linked Biopolymer Networks. Soft Matter 2013, 9, 7302−7313. (12) Sun, Y.-L.; Luo, Z.-P.; Fertala, A.; An, K.-N. Stretching Type II Collagen with Optical Tweezers. J. Biomech. 2004, 37, 1665−1669. (13) Stein, A. M.; Vader, D. A.; Weitz, D. A.; Sander, L. M. The Micromechanics of Three-dimensional Collagen-I Gels. Complexity 2011, 16, 22−28.

I

DOI: 10.1021/acs.jpcb.6b00259 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (14) Yang, L.; Van der Werf, K. O.; Fitié, C. F.; Bennink, M. L.; Dijkstra, P. J.; Feijen, J. Mechanical Properties of Native and Crosslinked Type I Collagen Fibrils. Biophys. J. 2008, 94, 2204−2211. (15) MacKintosh, F. C.; Käs, J.; Janmey, P. A. Elasticity of Semiflexible Biopolymer Networks. Phys. Rev. Lett. 1995, 75, 4425− 4428. (16) Morse, D. C. Viscoelasticity of Concentrated Isotropic Solutions of Semiflexible Polymers. 2. Linear Response. Macromolecules 1998, 31, 7044−7067. (17) Gardel, M. L.; Shin, J. H.; MacKintosh, F. C.; Mahadevan, L.; Matsudaira, P.; Weitz, D. A. Elastic Behavior of Cross-linked and Bundled Actin Networks. Science 2004, 304, 1301−1305. (18) Storm, C.; Pastore, J. J.; MacKintosh, F. C.; Lubensky, T. C.; Janmey, P. A. Nonlinear Elasticity in Biological Gels. Nature 2005, 435, 191−194. (19) Wagner, B.; Tharmann, R.; Haase, I.; Fischer, M.; Bausch, A. Cytoskeletal Polymer Networks: The Molecular Structure of Crosslinkers Determines Macroscopic Properties. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 13974−13978. (20) Tharmann, R.; Claessens, M.; Bausch, A. Viscoelasticity of Isotropically Cross-linked Actin Networks. Phys. Rev. Lett. 2007, 98, 088103(1−4). (21) Janmey, P. A.; McCormick, M. E.; Rammensee, S.; Leight, J. L.; Georges, P. C.; MacKintosh, F. C. Negative Normal Stress in Semiflexible Biopolymer Gels. Nat. Mater. 2007, 6, 48−51. (22) Lin, Y.-C.; Yao, N. Y.; Broedersz, C. P.; Herrmann, H.; MacKintosh, F. C.; Weitz, D. A. Origins of Elasticity in Intermediate Filament Networks. Phys. Rev. Lett. 2010, 104, 058101(1−4). (23) Broedersz, C. P.; MacKintosh, F. C. Modeling Semiflexible Polymer Networks. Rev. Mod. Phys. 2014, 86, 995−1036. (24) Licup, A. J.; Münster, S.; Sharma, A.; Sheinman, M.; Jawerth, L. M.; Fabry, B.; Weitz, D. A.; MacKintosh, F. C. Stress Controls the Mechanics of Collagen Networks. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 9573−9578. (25) Das, M.; MacKintosh, F.; Levine, A. J. Effective Medium Theory of Semiflexible Filamentous Networks. Phys. Rev. Lett. 2007, 99, 038101(1−4). (26) Broedersz, C. P.; Mao, X.; Lubensky, T. C.; MacKintosh, F. C. Criticality and Isostaticity in Fibre Networks. Nat. Phys. 2011, 7, 983− 988. (27) Broedersz, C. P.; Sheinman, M.; MacKintosh, F. C. Filamentlength-controlled Elasticity in 3D Fiber Networks. Phys. Rev. Lett. 2012, 108, 078102(1−5). (28) Feng, J.; Levine, H.; Mao, X.; Sander, L. M. Nonlinear Elasticity of Disordered Fiber Networks. arXiv preprint arXiv:1507.07519 2015. (29) Mao, X.; Stenull, O.; Lubensky, T. C. Effective-medium Theory of a Filamentous Triangular Lattice. Phys. Rev. E 2013, 87, 042601(1− 14). (30) Mao, X.; Stenull, O.; Lubensky, T. C. Elasticity of a Filamentous Kagome Lattice. Phys. Rev. E 2013, 87, 042602(1−14). (31) Feng, J.; Levine, H.; Mao, X.; Sander, L. M. Alignment and Nonlinear Elasticity in Biopolymer Gels. Phys. Rev. E 2015, 91, 042710(1−9). (32) Head, D. A.; Levine, A. J.; MacKintosh, F. C. Deformation of Cross-Linked Semiflexible Polymer Networks. Phys. Rev. Lett. 2003, 91, 108102(1−4). (33) Wilhelm, J.; Frey, E. Elasticity of Stiff Polymer Networks. Phys. Rev. Lett. 2003, 91, 108103(1−4). (34) Onck, P. R.; Koeman, T.; Van Dillen, T.; Van der Giessen, E. Alternative Explanation of Stiffening in Cross-linked Semiflexible Networks. Phys. Rev. Lett. 2005, 95, 178102(1−4). (35) Heussinger, C.; Schaefer, B.; Frey, E. Nonaffine Rubber Elasticity for Stiff Polymer Networks. Phys. Rev. E 2007, 76, 031906(1−12). (36) Conti, E.; MacKintosh, F. C. Cross-linked Networks of Stiff Filaments Exhibit Negative Normal Stress. Phys. Rev. Lett. 2009, 102, 088102(1−4).

(37) Sharma, A.; Sheinman, M.; Heidemann, K.; MacKintosh, F. C. Elastic Response of Filamentous Networks with Compliant Crosslinks. Phys. Rev. E 2013, 88, 052705(1−9). (38) Heidemann, K. M.; Sharma, A.; Rehfeldt, F.; Schmidt, C. F.; Wardetzky, M. Elasticity of 3D Networks with Rigid Filaments and Compliant Crosslinks. Soft Matter 2015, 11, 343−354. (39) Maxwell, J. C. L. On the Calculation of the Equilibrium and Stiffness of Frames. Philos. Mag. 1864, 27, 294−299. (40) Yang, Y.-l.; Leone, L. M.; Kaufman, L. J. Elastic Moduli of Collagen Gels Can Be Predicted From Two-dimensional Confocal Microscopy. Biophys. J. 2009, 97, 2051−2060. (41) Sheinman, M.; Broedersz, C. P.; MacKintosh, F. C. Actively Stressed Marginal Networks. Phys. Rev. Lett. 2012, 109, 238101(1−5). (42) Dennison, M.; Sheinman, M.; Storm, C.; MacKintosh, F. C. Fluctuation-stabilized Marginal Networks and Anomalous Entropic Elasticity. Phys. Rev. Lett. 2013, 111, 095503(1−5). (43) Kroy, K.; Frey, E. Force-extension Relation and Plateau Modulus for Wormlike Chains. Phys. Rev. Lett. 1996, 77, 306−309. (44) Satcher, R. L., Jr; Dewey, C. F., Jr Theoretical Estimates of Mechanical Properties of the Endothelial Cell Cytoskeleton. Biophys. J. 1996, 71, 109−118. (45) Broedersz, C.; MacKintosh, F. Molecular Motors Stiffen Nonaffine Semiflexible Polymer Networks. Soft Matter 2011, 7, 3186− 3191. (46) Licup, A.; Sharma, A.; MacKintosh, F. C. Elastic Regimes of Subisostatic Athermal Fiber Networks. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2016, 93, 012407(1−12). (47) Sahimi, M.; Arbabi, S. Mechanics of Disordered Solids. II. Percolation on Elastic Networks with Bond-bending Forces. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47, 703−712. (48) May, S.; Bohbot, Y.; Ben-Shaul, A. Molecular Theory of Bending Elasticity and Branching of Cylindrical Micelles. J. Phys. Chem. B 1997, 101, 8648−8657. (49) Lin, Y.; Wang, Y.; Zheng, J.; Yao, K.; Tan, H.; Wang, Y.; Tang, T.; Xu, D. Nanostructure and Linear Rheological Response of Comblike Copolymer PSVS-g-PE Melts: Influences of Branching Densities and Branching Chain Length. Macromolecules 2015, 48, 7640−7648. (50) Silver, F. H.; Freeman, J. W.; Seehra, G. P. Collagen Selfassembly and the Development of Tendon Mechanical Properties. J. Biomech. 2003, 36, 1529−1553. (51) Cavalcante, F. S.; Ito, S.; Brewer, K.; Sakai, H.; Alencar, A. M.; Almeida, M. P.; Andrade, J. S.; Majumdar, A.; Ingenito, E. P.; Suki, B. Mechanical Interactions Between Collagen and Proteoglycans: Implications for the Stability of Lung Tissue. J. Appl. Physiol. 2004, 98, 672−679. (52) Kratky, O.; Porod, G. Rö ntgenuntersuchung Gelö ster Fadenmoleküle. J. R. Neth. Chem. Soc. 1949, 68, 1106−1122. (53) Lindström, S. B.; Vader, D. A.; Kulachenko, A.; Weitz, D. A. Biopolymer Network Geometries: Characterization, Regeneration, and Elastic Properties. Phys. Rev. E 2010, 82, 051905(1−5). (54) Press, W. H. Numerical Recipes: The Art of Scientific Computing, 3rd ed.; Cambridge University Press, 2007; pp 515−520. (55) Lees, A. W.; Edwards, S. F. The Computer Study of Transport Processes Under Extreme Conditions. J. Phys. C: Solid State Phys. 1972, 5, 1921−1929. (56) Arrott, A.; Noakes, J. E. Approximate Equation of State For Nickel Near its Critical Temperature. Phys. Rev. Lett. 1967, 19, 786− 789. (57) DiDonna, B. A.; Lubensky, T. C. Nonaffine Correlations in Random Elastic Media. Phys. Rev. E 2005, 72, 066619(1−23. (58) Liu, J.; Koenderink, G. H.; Kasza, K. E.; MacKintosh, F. C.; Weitz, D. A. Visualizing the Strain Field in Semiflexible Polymer Networks: Strain Fluctuations and Nonlinear Rheology of F-Actin Gels. Phys. Rev. Lett. 2007, 98, 198304(1−4). (59) Hatami-Marbini, H.; Picu, R. C. Scaling of Nonaffine Deformation in Random Semiflexible Fiber Networks. Phys. Rev. E 2008, 77, 062103(1−4). (60) Alexander, S. Amorphous Solids: their Structure, Lattice Dynamics and Elasticity. Phys. Rep. 1998, 296, 65−236. J

DOI: 10.1021/acs.jpcb.6b00259 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (61) Ashkin, J.; Teller, E. Statistics of Two-Dimensional Lattices with Four Components. Phys. Rev. 1943, 64, 178−184. (62) Baxter, R. J. Eight-Vertex Model in Lattice Statistics. Phys. Rev. Lett. 1971, 26, 832−833. (63) Kadanoff, L. P.; Brown, A. C. Correlation Functions on the Critical Lines of the Baxter and Ashkin-Teller Models. Ann. Phys. 1979, 121, 318−342. (64) Butch, N. P.; Maple, M. B. Evolution of Critical Scaling Behavior near a Ferromagnetic Quantum Phase Transition. Phys. Rev. Lett. 2009, 103, 076404(1−4). (65) Fuchs, D.; Wissinger, M.; Schmalian, J.; Huang, C.-L.; Fromknecht, R.; Schneider, R.; Löhneysen, H. v. Critical Scaling Analysis of the Itinerant Ferromagnet Sr1−xCax RuO3. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 174405(1−7). (66) Kühn, R. Critical Behavior of the Randomly Spin Diluted 2D Ising Model: A Grand Ensemble Approach. Phys. Rev. Lett. 1994, 73, 2268−2271. (67) Calabrese, P.; Parruccini, P.; Pelissetto, A.; Vicari, E. Crossover Behavior in Three-dimensional Dilute Spin Systems. Phys. Rev. E 2004, 69, 036120(1−16).

K

DOI: 10.1021/acs.jpcb.6b00259 J. Phys. Chem. B XXXX, XXX, XXX−XXX