Statistical Mechanics of a Cat's Cradle

It is believed that, much like a cat's cradle, the cytoskeleton can be thought of as a network of strings under tension. We show that both regular and random bond-disordered networks having bonds that buckle upon compression exhibit a variety of phase transitions as a function of temperature and extension. The results of self-consistent phonon calculations for the regular networks agree very well with computer simulations at finite temperature. The analytic theory also yields a rigidity onset (mechanical percolation) and the fraction of extended bonds for random networks. There is very good agreement with the simulations by Delaney et al. (Europhys. Lett. 2005). The mean field theory reveals a nontranslationally invariant phase with self-generated heterogeneity of tautness, representing ``antiferroelasticity''.

It is thought that the cytoskeleton consists of a polymeric network which is maintained under tension [1,2,3], much like the child's game made of strings known as a "cat's cradle". Unlike a macroscopic "cat's cradle", the cytoskeleton is subject to fluctuating forces from both thermal motions and the nonequilibrium noise due to motor proteins and constant polymerization and depolymerization events [4].
Upon freezing molecular fluids that resist compression undergo transitions to rigid solids. Here we show that, a caricature of the cytoskeleton, a model of interconnecting ropes that easily buckle also undergoes a variety of phase transitions. We study networks of ropes that bear tension when they are stretched beyond a critical threshold, but otherwise buckle. While the model's interactions contrast with the usual compression-bearing mechanical interaction between atoms or molecules, the "inverted" nonlinearity of ropes still leads to an onset of rigidity like that in crystals or glasses, but when the system is expanded rather than compressed. Depending on where we set the baseline of environmental pressure or tension, we may also say the system buckles upon compression.
The cytoskeleton actually has both tension-bearing and compressional elements. The present theory, can be extended to include both elements. Elsewhere we have discussed the nonequilibrium aspects [5] that are also needed to capture the essential features of the tensegrity model [2,3] of cell mechanics [6], but here we limit ourselves to the effects of equilibrium fluctuations.
Delaney, Weaire, and Hutzler pioneered the study of such buckling lattices with neighboring sites that interact with an "inside-out" potential u(r; l) = Θ(r−l)· k 2 ·(r−l) 2 , where Θ() is the Heaviside unit step function [7]. The energy stored in a bond vanishes when the length of the bond r is less than a critical value l, but increases quadratically beyond the cutoff l. A two (three) dimensional network of bonds is thermodynamically characterized by its tension (negative pressure), as a function of its area (volume) −p 2 = ∂F/∂s (−p 3 = ∂F/∂v).
In the low temperature limit, the entropy term can be ignored and elementary statics suffices to describe regular networks. The square network, for example, has a tension as a function of the area expansion given by −p sq 2 (s/s 0 ) = k(1 − s 0 /s) Θ(s/s 0 − 1). Rigidity only occurs for areas exceeding s 0 . Tensions for other lattices can also be computed p tri 2 = 2 √ 3p sq 2 and p hex 2 = p tri 2 /3 for the triangular and the hexagonal network respectively. For the 3D cubic network one finds at T = 0, −p sc yielding a rigidity threshold at v 0 . We have set s 0 (or v 0 ) as the area (or the volume) of the corresponding network constructed from all bonds at the unit length. The tension of 3D networks has a maximum as a function of volume, because at large volume the energy density of the system decreases as v −1/3 . At high temperature or equivalently for soft springs, even without disorder, static models are insufficient and one always needs a statistical mechanical description. If the network is disordered, i.e., having a broad distribution of onset length l, finding the tension is a nontrivial statistical problem even in the limit of low temperature [7].
To analytically treat statistical cat's cradles, we adopt the self-consistent phonon technique, initially developed for anharmonic lattice solids with compression bearing interactions [8,9]. In this approach one replaces the pairwise interaction in the equilibrium density matrix exp[−β ij u ij ( r i − r j )] by a nearly equivalent effective single-body term as Decoupling the manybody effects can be carried out by assuming a self-consistent one-body potentialũ has a quadratic form The self-consistent phonon method has not only been applied to crystals but also to noncrystalline solids [10], network glasses [11], and to problems of nonequilibrium structural assembly [12]. The self-consistent phonon method can be made the beginning of a systematic expansion [8], but already has been shown to work very well for systems made up of compressive elements.
We study both 2D and 3D systems. Long wavelength phonons lead to a diverging mean square displacement (MSD) with increasing system size for 2D systems, in a strict sense. However this divergence is weak (logarith-mic) and there is still a crossover reflecting rather rapid decrease in fluctuations of the mean square nearest neighbor distance upon extension. In fact, the self-consistent phonon method, when used for the 2D system predicts the correct onset of rigidity for the random 2D network at T = 0 as we shall see.
In the self-consistent phonon theory the partition function z(ξ) = [ d r exp(−ǫ(ξ) − α(ξ)r 2 )] N is evaluated as a function of expansion ratio ξ = s/s 0 or v/v 0 . Both the on-site free energy ǫ and the phonon frequency α depend on the ratio of the mean neighbor distance R to the critical extension l, R/l = D √ ξ. We can use the chain rule to calculate tension, e.g., −p 2 = ∂f ∂s = − ∂ ln z(ξ(R)) β∂s(R) . In three dimensions the Mayer f -bond v( R; l, α) for the interaction u(r; l) when averaged over the Gaussian fluctuation of a neighboring particle around its mean position R having its own spring constant α gives This can be expressed using the error function erf(). The 2D averaged Mayer f -bond can be expressed as a numerical integral of modified Bessel function I 0 .
An interesting property that can be calculated within the self-consistent phonon theory is the fraction of the bonds bearing force (that are taut), i.e., having a length larger than the intrinsic cutoff l. These are functions of α and R: The self-consistent phonon frequency α for several different regular networks is shown in the left panels of Fig. 1. Unlike the situation for the systems made of hardspheres [5], there is always some nonzero self-consistent solution α. Entropy always leads to some rigidity. There is, however, a rapid crossover from low α to high α accompanying the expansion of the system. Generally either high coordination number z or larger microscopic rigidity-temperature ratio b (= kβ/2) lead to a larger α. The fraction of "loose" bonds, (1 − q) predicted from the phonon theory is plotted in the right panels of Fig. 1. This fraction decreases as the network expands.
We have shown there is a distinct bifurcation of α that also occurs when the volume ratio ξ becomes small enough for 3D networks with high rigidity-temperature ratio b, This bifurcation appears in all three types of 3D networks studied here. The fcc networks most clearly display this bifurcation, followed by the bcc, and then the sc networks at smaller ξ and higher b. As shown in Fig. 2, usually there is only one stable solution of the self-consistent equation α tagged = f (α neighbor ). However, under certain conditions (small volume and large rigiditytemperature ratio, or pictorially, when the slope at the intercept is not in the range of between -1 and 1), the self-consistent solution becomes unstable, and another pair of solutions α l and α h satisfying α l = f (α h ) and α h = f (α l ) emerges as the stable solution. This stable pair solution indicates the possibility of having a nontranslationally invariant solution, in which a low α particle is surrounded by high α neighbors and vice versa. The resulting spatial heterogeneity is self-generated. In the bcc and sc networks, this solution can be easily interpreted physically: there is a symmetry-broken phase of bcc or sc networks that can be thought as an "antiferromagnetic" (AF) one with the underlying lattice being split to two subnetworks (sc → fcc or bcc → sc): one of the two equal subnetworks bears more tension than the other. We call this the "antiferroelastic" phase. In the fcc case, since there are four particles in the unit cell, one high α particle cannot be completely surrounded by low α particles. Nevertheless one can have two possible ordered AF ground states similar to the crystal structures of AuCu and AuCu 3 . To test this possibility, we set up a more elaborate self-consistent solution having four values of α = (α 0 , α x , α y , α z ) at (000), (011), (101), (110) respectively, α 0 = f (α x , α y , α z ), α x = f (α 0 , α z , α y ), α y = f (α z , α 0 , α x ), α z = f (α y , α x , α 0 ). By explicitly labeling all four particles in a unit cell (as done similarly in the AF Ising fcc system [13]), we searched but failed to find these explicit symmetry-broken solutions. Thus we believe this frustration of fcc network would lead to a random phase which is analogous to a self-generated glass. This possibility may potentially be resolved using a replica calculation [14,15].
We checked our self-consistent phonon theory using explicit computer simulations of the network with thermal noise. We show the results for the bcc case in Fig. 3. The simulations were performed with stochastic dynamics on n 3 unit cells with periodic boundary condition for various rigidity-temperature ratios b and nonlinearity parameters l for various values of n. We see the self-consistent phonon theory agrees quite well with the simulation. For the simulation of size 5 3 , the errors for the MSD of the particle position between theory and simulation are only about 1%. We extrapolated to the thermodynamic limit based on a series of simulation of different sizes of networks. A Similar good agreement is achieved for 2D networks when the MSD of relative neighbor position from simulation is compared with theory.
To check the prediction of spontaneous symmetry breaking of neighboring α values requires a large ratio of l/R ∼ 10 necessitating simulation using a very large system at very low temperature. We estimate the necessary system size to be at least 10 6 . We were unable to carry out such a large scale simulation.
We plot the tension vs. extension in Fig. 4 (a) and (b). Above a critical tension, there are two solutions for p(s/s 0 ) = p. The solution with ∂p/∂s < 0 is thermodynamically stable for a constant pressure ensemble.
The system always appears to be stable in two dimensions. For 3D networks, As shown in Fig 4(b), there is single (unstable) solution at high temperature for constant p 3 ; but there are three (one stable and two unstable) solutions at low temperature. The values of the critical rigidity-temperature ratio for this transition are b * = 1.04, 0.76, and 0.48 for sc, bcc, and fcc networks respectively. Thus the network is intrinsically unstable for the constant p ensemble for high temperature or low spring constant, i.e. b < b * , regardless of supplied tension −p. For b > b * , there is a critical tension −p * above which the system is no longer stable. Generally a larger b implies a larger stable region.
In the large s/s 0 or v/v 0 limit, the Dulong-Petit law provides pressure that can be based on dimension counting alone, giving p 2 − p o = N k B T /(2S) and p 3 − p o = N k B T /(2V ). Here p o = p(T = 0) is the pressure from the mechanical static calculations at zero temperature. The factor 1 2 arises because only the potential energy is counted here, the kinetic energy counts the other half. Thus for large expansion ratios, p 2 → const. and p 3 → 0, which are consistent with the tension curves shown in Fig. 4

(a) and (b).
So far we have focussed on the effects of thermal disorder alone. We now discuss the effects of bond disorder. For illustration, we concentrate on the 2D random network, especially in the βk → ∞ limit. Following Delaney et al. [7], we studied a "flat" random distribution of string lengths. More general distributions present no new challenges to our method.
As shown in Fig. 4(c) for the triangular network, phonon frequencies α rise quickly with expansion at zero temperature when b = ∞. The smallest area s * /s 0 satisfying α −1 (s/s 0 ) = 0 is the critical value for the onset of macroscopic rigidity (mechanical percolation). Since the spring constant k and temperature are always bundled together in the treatment, once the system is taut, α → ∞ for 1/b = 0. For finite temperature, α is always finite. Th rigidity onset s * /s 0 calculated here quantitatively agrees with the simulation results in Ref. [7] as a function of width of distribution of string length w. As expected, a larger w means a lower α generally. We also display in Fig. 4 the fraction of vertices with given number of taut bonds P n (d) and the mean fraction of taut bonds q (e). These properties have also been obtained by the simulations of Delaney et al. We can only make a qualitative comparison for these two properties since our calculations were performed with finite b, corresponding to T = 0 while the simulations were strictly done at T = 0. Our theory does not allow us to calculate the loose bond fraction at T = 0 due to the bundling of β and k. An approximation was made for the calculation of p(n), the fraction of vertices with given number n bonds taut based on the mean fraction of taut bonds q, i.e., p(n) = Z! n!(Z−n)! q n (1 − q) Z−n . Here Z is the coordination number.
As shown in Fig. 4(f), for bond-disordered 3D networks, again there is a rigidity onset, which is similar to those of 2D networks. Among the cases studied only the . The effects of the random bond distribution w on the phonon frequency α are shown in (c) with the inversion 1/α shown in the inset of (c) to better illustrate the critical area for onset of rigidity. The portion of vertices with n bonds tauted Pn and the mean taut ratio q are shown in (d) and (e) respectively. We show the inversion 1/α as functions of volume for various 3D random networks in (f). w = 1 fcc network showed the antiferroelastic structural inhomogeneity bifurcation. The self-consistent phonon method suggests random bond-disorder destroys antifer-roelasticity. This work was supported in part by NSF-CTBP.