Simulated Cytoskeletal Collapse via Tau Degradation

We present a coarse-grained two dimensional mechanical model for the microtubule-tau bundles in neuronal axons in which we remove taus, as can happen in various neurodegenerative conditions such as Alzheimers disease, tauopathies, and chronic traumatic encephalopathy. Our simplified model includes (i) taus modeled as entropic springs between microtubules, (ii) removal of taus from the bundles due to phosphorylation, and (iii) a possible depletion force between microtubules due to these dissociated phosphorylated taus. We equilibrate upon tau removal using steepest descent relaxation. In the absence of the depletion force, the transverse rigidity to radial compression of the bundles falls to zero at about 60% tau occupancy, in agreement with standard percolation theory results. However, with the attractive depletion force, spring removal leads to a first order collapse of the bundles over a wide range of tau occupancies for physiologically realizable conditions. While our simplest calculations assume a constant concentration of microtubule intercalants to mediate the depletion force, including a dependence that is linear in the detached taus yields the same collapse. Applying percolation theory to removal of taus at microtubule tips, which are likely to be the protective sites against dynamic instability, we argue that the microtubule instability can only obtain at low tau occupancy, from 0.06–0.30 depending upon the tau coordination at the microtubule tips. Hence, the collapse we discover is likely to be more robust over a wide range of tau occupancies than the dynamic instability. We suggest in vitro tests of our predicted collapse.


Introduction
Neurodegenerative conditions such as Alzheimer's disease (AD), Parkinson's, and those resulting from chronic traumatic encephalopathy (CTE) damage arising in contact sports, represent a massive public health threat that annually impacts tens of millions worldwide. Finding routes to prevention or treatment prior to irreversible changes in the affected cells remain important goals. Part of the difficulty with this task is the complexity of the conditions themselves. For example, in Alzheimer's disease, while it is generally agreed, per the ''amyloid cascade hypothesis'' (see [1] and [2]) that the initial trigger of AD is the production and subsequent aggregation of A-beta protein (ABP) into oligomers (ABO), the way in which the ABOs trigger cell damage and eventual death remains unclear.
Neuro-fibrillary tangles (NFTs), consisting largely of hyperphosphorylated tau protein aggregates, have been observed in postmortem tissues of AD victims. These NFTs are correlated with regions of A-beta aggregates and neuronal damage or death [3] confirms that ABOs trigger processes which lead to modifications of the tau protein. NFTs and tau associated damage also arise in CTE [4]. Tau proteins play a critical role in the microtubule bundles (MTBs) of the axon cytoskeleton, by promoting the assembly of axonal microtubules and inhibiting dynamic instability (DI) [5], and by crosslinking the microtubules in the bundles. The integrity of the MTBs is critical for neuronal function: they provide mechanical stability to axons for decades of human life, and they serve as the highways for fast and slow axonal transport of neurotransmitters and organelles to and from the synaptic regions of the neuron. Clearly, at the cellular level, a significant marker of irreversible damage in conditions affecting tau proteins would be loss of integrity of the MTBs, and thus a point of no return for intervention.
Studies of cultured neurons exposed to high concentrations of ABOs have confirmed that the ABOs can be internalized, and once internalized these trigger production of proteinases which can cut taus, and kinases which can phosphorylate taus [6]. Both are important, since full length taus are needed to assure the proper spacing of microtubules, and hyperphosphorylated taus (HPTs) may tend to dissociate from the negatively charged microtubule surfaces.
If tau detachment from microtubule binding sites exceeds 50%, this should lead to collapse of MTBs due to the dynamic instability of individual microtubules. In experiments on cultured hippocampal neurons exposed to massive concentrations of ABOs (1000-10000 times physiological ABP concentrations in vivo) MT loss is indeed observed [7], with the lateral density decreasing by a factor of 25 compared to unexposed cells. However, these high ABO concentrations may obscure more subtle changes in the MTBs. Specifically, since the taus provide transverse mechanical stability, MTB transverse rigidity may be lost when a majority of the taus are still intact: random removal of connecting springs in a lattice of otherwise noninteracting cylinders leads to a loss of transverse mechanical stiffness at the rigidity percolation threshold near 60% of spring occupancy as we discuss below. Moreover, as the HPTs no longer favorably associate with the surface of microtubules, they can lower the free energy by gaining translational entropy through reduced excluded volume between microtubules, yielding a fluctuation mediated attraction between microtubules. This will, in effect, add an inward radial pressure to the MTBs as taus are removed.
In this paper, we develop a two dimensional mechanical model for the MTBs and consider the percolation transition as we remove taus between the microtubules (Fig. 1). Our method confirms that in the absence of the depletion force, the transverse rigidity is lost when the concentration of tau binding is reduced to nearly 60% of the original value. We denote this ratio of bound tau concentration to the maximum value as the tau occupancy p: Adding the depletion force via a simple excluded volume model leads to a first order collapse that depends upon the ratio of the tau spring force to the concentration of HPT. This collapse arises because the steady inward pressure of the nearest neighbor MT mediated depletion force contracts the microtubule bundles to a point where next-nearest-neighbor microtubules can experience attraction. For physiologically attainable values of tau spring force and HPT concentrations, the collapse transition can easily arise above 40% bound tau occupancy, while we argue that the dynamic instability is likely to occur below this value. Hence, with or without the depletion force we believe there will be irreversible mechanical damage to the MTBs at high bound tau concentrations, well before the dynamical instability is relevant. In formulating treatment strategies for tau related diseases, these conclusions could be very important.

Overview of Model
Consistent with experiment, the MTB is modeled initially as a simple hexagonal lattice of one micron long rigid microtubules of radius R~12:5nm (blue disks in Figure 1) linked by tau springs (red lines). There are approximately l = 50 tau springs per micron in a microtubule pair with spring constant k eff , and the mean microtubule center to center distance is taken to be 45 nm [8]. The springs are modeled harmonically in this initial simulation where compression is small on a per spring basis. We have simulated five MTB sizes of varying initial radii R mtb (1) at full tau occupancy p~1: R mtb (1)~100, 150, 200, 250, 300 nm with, respectively, 19, 37, 73, 117, and 163 microtubules. The radius is measured as the maximum center-to-center MT separation from the central MT to a boundary one for the hexagonal shaped bundles. Because each 50 nm increase in radius adds one hexagonal shell of microtubules, we describe bundle sizes using ''n-shell'' where n = 2, 3, 4, 5, 6 for R mtb (1). We neglect surface tension at the boundary due to the combined action of the axon membrane and actin cytoskeleton, as direct measurement shows this to be small for axons [9]. We employ the steepest descent relaxation algorithm (equivalent to highly damped molecular dynamics), where microtubules are iteratively moved a distance proportional to the net force experienced until the positions are converged to within a fractional tolerance factor [10]. We present details of the relaxation algorithm in the Supplemental Information. The ratio of initial resting spring length ' 0 to lattice spacing a is denoted g, which we vary in the case of no depletion force. We note that we are not including neurofilaments in our model, which are not critical for stability of the MTBs (though they are for axons) [11].
The taus are intrinsically disordered proteins with a single microtubule binding domain. Hence, they must be at least in dimer form to link two microtubules [12]. The N-terminus of the tau is predominantly negative, while the middle (M) region preceding the microtubule binding domain is predominantly positive. As argued in Ref. [12], the tau can dimerize in this region through salt bridges between the N/M regions, which take up about 200 residues, giving a total peptide length to this domain of L NM~2 00a R &80 nm where a R is the length of a residue, about 0.4 nm. The simplest model for the spring constant, in the low extension/compression regime, is of a paired set (one per monomer) of (n s {1) entropic springs of length ' s~LNM =(n s {1) between n s salt bridges where we assume the unstructured entropic spring peptides have a short persistence length j p^1 nm. The result is (c.f., this result in Phillips et al.) [13] We thus estimate k eff &0.05 pN/nm. Proteins with well defined secondary structure of course have much higher effective spring constants under compression, potentially reaching more than 10 pN/nm, as can be inferred from works on green flourescent protein [14] and fibrinogen [15] for example. Figure 2 shows the excluded volume model for the depletion force, derived for cylindrical geometry following the spherical example of Phillips et al. [16]. Once taus are phosphorylated, they no longer associate with the microtubules within the annulus between R, Rzr. To leading order in the volume difference, the effective interaction potential V D (Z)=L per unit length L between a pair of microtubules with center-to-center distance Z and Z=(2(Rzr))~f is given by Here P O is the osmotic pressure of the nonbinding hyperphosphorylated tau dimers given by where k B is Boltzmann's constant, T the temperature, and t 2p Â Ã the concentration of hyperphosphorylated tau dimers. The corresponding depletion force is taken as the negative gradient of V D with respect to Z: This gives a depletion force per unit length F D (Z)=L of In our simulations, of course the tau dimers are not spherical objects, and their precise radius of gyration is not known, although we anticipate the dimers to be at least 20 nm long (the wall-to-wall distance between microtubules) and 5-10 nm thick. For simplicity we have taken r~R = 12.5 nm. We find, in the region of collapse(Z&35nm), that the magnitude of the depletion force per unit length divided by the osmotic pressure w D (Z) is weakly dependent upon r and given approximately by w D varies only from 256 nm to 285 nm as r varies from 10-15 nm. Hence, the main determinant of the depletion force is the magnitude of the osmotic pressure. We note that taking the radius of gyration of the tau dimers to be 10-15 nm makes them larger than the observed (calculated) value of 6.5(6.9) nm for tau monomers [17]. However, these monomers are in a compressed conformation in which the positive and negative charges comingle. The relatively extended conformation necessary for the dimer must be larger.
The upper limit on the osmotic pressure associated with hyperphosphorylated tau dimers is clearly the bound concentration of dimers, t 2b ½ , in the nondegraded axon. Using the volume between microtubules, we estimate t 2b ½ = 320 micromolar, which gives an osmotic pressure of 800 Pa. We note that the density of bound taus in the MT bundles is significantly larger than the^2-4 micromolar estimated concentrations of free tau monomers, as has been observed before [18].
It is convenient to characterize our model bundles with depletion force in terms of dimensionless ratio r of the spring force constant per unit length compared to the osmotic pressure, which is given by For l = 50/micron, k eff = 0.05 pN/nm, and P O = 800 Pa, r~3:125.

Test Case: simulation of taxol stabilized microtubules with PEO intercalants
We note that the depletion force has been observed in in vitro experiments with taxol stabilized microtubule bundles intercalated with PEO polymers (and no taus) [19]. In this case the only additional interaction between the micotubules is via screened Coulomb coupling, which is estimated by the force per unit microtubule surface area P ES with the Debye length L D = 1.47 nm. We add this force to our depletion force to test our steepest descents equilibration as shown in Fig. 3. Over a factor of 100 variation in P O , the agreement between theory and the data in Fig. 11 of Needleman et al. [19] is good. Note that given the large rest length of the taus, the interaction P ES is irrelevant prior to any collapse and so we neglect it in our model for neuronal microtubules.

Model MT repulsion for axonal MTBs
In our simulations with taus for axonal MTBs, we mimic the combined screened electrostatic repulsion/steric repulsion between microtubules by a Morse potential of the form where we take D e~1 0 {4 J and a = 1 nm 21 . While this potential does provide a small minimum near the microtubule surfaces, its main effect is to inhibit collapse beyond the region where microtubules are in contact, which is sufficient for our purposes.

Test case: Rigidity Percolation with no Depletion Force
In the depletion force free case, we have computed the second order elastic constant C (Fig. 4) which measures the stiffness of the microtubules against transverse compression for a variety of different ratios g of the initial microtubule-to-microtubule separation to spring resting length. In all cases, the scaled curve for C(p,g)=C(1,1) is identical, since unlike the work of Ref. [20] we do not hold our springs in a fixed frame under tension but rather allow them to relax, so the result is essentially equivalent to the case where bond length equals rest length in Ref. [20]. This gives a transverse rigidity collapse (zero transverse elastic constant) at a critical tau occupancy p c &0:63 in agreement with Ref. [20]. Note that we do not see a sharp threshold due to finite size effects; we have confirmed that the rounding of the transition is reduced as a function of increasing system size. We stress here that by vanishing elastic constant, we mean that the MTB will initially lose rigidity against radial compression, until the MTs can come into direct steric contact, at which point the elastic constant will jump to a value suitable for close packed MTs. Hence, even without inclusion of the depletion force, we anticipate a significant mechanical degradation (transverse rigidity loss) at a high p value which exceeds the range (pƒ0:3) argued below to be relevant for the dynamic instability mechanism.

MTB Collapse with Depletion Force
With the depletion force present, we find a first order volume collapse as p is reduced, as shown in Fig. 5, for one value of the dimensionless parameter r. This figure shows a typical plot of the normalized transverse radius of gyration of the microtubule bundle R MTB (p)=R MTB (1) decreasing slowly with decreasing p until the critical percolation threshold p c is reached, at which point R MTB drops precipitously to the hexagonal closed packed steric limit R MTB (p c )~2nR~25n nm, where as defined above, n is the number of hexagonal shells. The parameter r controls the position of the collapse, as shown in Fig. 6(a), where we show curves for pwp c such that the left endpoint is at p c : This collapse is independent of the size of our starting bundle showing that it is not an artifact of small bundles in our simulation as shown in Fig. 6(b) where we plot the results for a single r value several bundle sizes. Each curve in Fig. 6 is averaged over five runs and is truncated at p~p c ; the horizontal bar located at the collapse region represents (by its width) the range of p c values obtained. Fig. 7 shows that holding r fixed while co-varying P O ,k eff yields the same collapse curves for pwp c . Note that below the collapse the compressional Simulated Cytoskeletal Collapse via Tau Degradation PLOS ONE | www.plosone.org rigidity will increase, as it is no longer dominated by the taus but by the MTs themselves (which become close packed below the collapse). We illustrate the collapse in the movie (MovieS1.mp4) provided in the supplemental material, which depicts the n~2 bundle as taus are removed.
In the majority of of our simulations we have taken the osmotic pressure constant at all p values, but of course if the depletion force is due to the phosphorylated tau dimers, we expect P O to be proportional to 1{p: The lavender curve in Fig. 6(a) shows the result for which the osmotic pressure is taken to be proportional to 1{p, reaching the same value as for the yellow curve at the collapse threshold p~p c : Clearly to within the range of p c values obtained over our different random number seeds, we obtain the same p c for constant and varying P O models. The lavender curve varies more substantially in radius, because at p?1, there is a different equilibrated radius suitable for zero depletion force.

Heuristic Mean Field Theory (MFT) of MTB Collapse
The origin of the collapse becomes clear if we consider a heuristic mean field theory (MFT) for the microtubule bundles. The assumptions which go into this are as follows: (i) as we remove taus we replace k eff with pk eff . (ii) The average symmetry around a given microtubule remains hexagonal, and the nearest neighbor and next nearest neighbor coordination numbers are 6. (iii) There is a hard wall when the microtubule center-to-center separation reaches 2R: (iv) Initially, the next nearest neighbor microtubules are at distances R NNN §2(Rzr), but as the system contracts, they reach within the depletion force active zone.
Because of the latter assumption, the potential energy acquires an extra ''kick'' when the next-nearest neighbor MTs enter the depletion zone radius. A simple model energy per unit length E MFT which captures this is where Z is the mean MT center-to-center separation and L is the MT length. Here E M FT is the total energy of the MTB within the mean field theory. The third term accounts for the next nearest neighbor interactions. We supplement this by an infinite wall at Zƒ2R.  [19]. The red circles are our data for osmotic pressures given in Ref. [19], the blue triangles are the data in that reference. doi:10.1371/journal.pone.0104965.g003 The key idea of the MFT is illustrated in Fig. 8. There, we show several curves for n MFT~EMFT =(2P O (Rzr) 2 ) vs. f~Z=2(Rzr) as p is varied. For pwp c , e.g., p~1, there is a minimum at an equilibrated value Z&1:6(Rzr) due to the balance of depletion attraction and spring repulsion. As shown in Fig 7 for 2(Rzr) = 50 nm and r = 3.125, the enforced hard wall minimum at Z~2R becomes degenerate with the large Z minimum as p is reduced to the critical value p c = 0.815. Below the critical value for p, the hard wall minimum is lower in energy. Hence, the transition at p c is a first order jump. In Fig. 9 we show a Figure 4. Transverse rigidity transition in the absence of the depletion force. As described in the text, we compute the 2nd order elastic constant for radial compression as a function of the occupied spring fraction p and (b) scale to the p~1,g~1:0 limit. In contrast to Ref. [20], all g values yield the same results when scaled, because their lattice boundary is fixed while ours relax to equilibrium. In all cases, the transverse rigidity percolation threshold is p c &0:6, in good agreement with the g~1:0 case of Ref. [20] doi:10.1371/journal.pone.0104965.g004 The driving feature for the collapse is the bootstrapping effect of diminished bundle radius driving increased depletion attraction which in turn pulls the next nearest neighbors in to overwhelm the resisting tau entropic springs.
Another important feature shared by the MFT and the full simulation is that a fixed r value driven by a fixed P O will produce an MTB collapse at the same p c as does a varying r(p) with P O (p)~P O (p c )(1{p)=(1{p c ): The latter model for the osmotic pressure is appropriate for the assumption that the removed taus themselves produce the depletion force. Competition with Dynamic Instability: Range of possible collapse dominance It has been presumed in previous work that the most likely route of degradation for the MTBs under tau degradation is via dynamic instability, in which the MTs stochastically vanish due to dominance of depolymerization when there are insufficient stabilizing taus. By augmenting our model with the following assumptions, we can produce a heuristic estimate for the p threshold below which dynamic instability should be important: Depolymerization initiates at the ends of MTs, so the critical stabilizing molecules binding to MTs should do so at the ends. (ii) Given ample evidence for two binding sites of tau on MTs [21], including one in the MT lumen [22], we assume that taus at the ends provide insurance against dynamic instability [23]. (iii) Evidence from studies with the small molecule eribulin suggest that binding of one molecule at the end is sufficient to inhibit dynamic instability [24].
Thus we assume that it is necessary to remove endpoint taus to generate dynamic instability. An estimate of where this will occur depends upon (i) the tau coordination at the ends, and (ii) the probability of site percolation on our lattice (if we remove a continuous backbone of microtubules through dynamic instability we lose lattice integrity). The site percolation threshold for two dimensional triangular lattices is 0.5 [25]. A conservative upper bound on the threshold for dynamic instability is thus to assume an end coordination of 1, and the probability of removing 2 end taus from a microtubule is thus (1{p) 2 : 50% of end taus would then be removed when (1{p DI ) 2~1 =2, giving p DIz~( 1{1= ffiffi ffi 2 p )~0:29. If we assume instead that the coordination is 6 at the ends, we obtain (1{p DI{ ) 12~1 =2 as the criterion, giving the occupancy of taus to generate DI at p DI {~( 1{2 {1=12 )~0:056.
In fact, our model is oversimplified by assuming that all taus bind equally to the MTs. The presence of two binding sites with different affinities (lumen binding taus have a higher affinity than surface taus) suggest that the stabilizing end taus are harder to remove with phosphorylation, so that the above estimates for the onset of DI are likely to be reduced further.
This leads us to conclude with confidence that the collapse scenario here is likely to occur over the dynamic instability for p §0:29, and is certainly very probable for p §0:06. With the more conservative criterion (removal of one tau from each end) then we find for r~15nm that the collapse is viable over DI for k eff ƒ0:14 pN/nm, while for r~10nm the collapse is viable for k eff ƒ0:067 pN/nm.

Discussion
Having demonstrated that a MTB collapse is likely to occur for a range of reasonable physical parameters under conditions that phosphorylated taus are removed from microtubules, it is reasonable to ask whether this can be observed or may have implications correlating with existing observations. First, as mentioned earlier, there is strong evidence of MT dynamic instability in cultured neurons exposed to massive doses of ABO over a period of hours, at 1-5 micromolar concentration [7] compared to estimated AB concentrations of .25-.5 nanomolar in vivo [26]. An exposure of mature cultured neurons over a period of days to more reasonable subnanomolar concentrations of AB dimers (but not monomers) does not induce microtubule damage at the resolution of the experiments but does cause beading of hyperphosporylated taus proximate to the microtubules and neuritic damage [27]. We conjecture that the high concentration ABO experiments induce massive phosphorylation of the taus, while the lower concentration exposures may not cause enough tau removal to induce our proposed MTB transverse collapse.
We propose two direct experimental tests to observe this novel collapse. First, we suggest using taxol stablized microtubule bundles in vitro per Ref. [19] with taus binding along the MT external surfaces, then exposing the MTBs to phosphorylating kinases. The taxols will inhibit dynamic instability, so this will be a clean test of whether the transverse collapse scenario proposed here is viable. Second, one can measure the response of the axon to ABO exposure in cultured neurons with microbead compression of the neurons. This has been done for initial ABO exposure on cultured axon free N2a cells and on HT22 mouse hippocampal cell nuclei where it has been shown that the mechanical properties of the membrane change are consistent either with increased osmotic pressure due to calcium influx or a stiffening of the membrane under 30 minutes of exposure to micromolar concentrations of ABO [28]. We would advocate monitoring of the axon directly over time with microbead compression to look for evidence of collapse in the force response curves. If the depletion force is irrelevant, then at some point the rigidity loss will lead to a change in the measured elastic constant. If there is a collapse then there will be a change in the monitored force in contact with the bead. Even at high exposures of ABO, if the collapse scenario is correct there should be an intermediate time scale in which the collapse is observable prior to the onset of MT dynamic instability. Studies of axonal compression under ABO exposure in late stages of cell life have not, to our knowledge, been carried out.
One observation of patients with AD is an atrophy of white (axonal) matter in affected regions of the brain. Functional MRI measurements of water diffusion anisotropy for patients with mild cognitive impairment gives evidence for axonal damage [29]. Additionally, it has been determined that changes in white matter volume of AD patients compared to elderly control patients can  2 ) vs. scaled mean microtubule separation f~Z=(2(Rzr)) for 2(Rzr)~50nm and r = 3.125. As p is reduced, the potential at f~0:5, which is the separation at microtubule contact, reaches a lower value than the minimum which evolves from f~0:82 at p~1. doi:10.1371/journal.pone.0104965.g008 reach values of 16.4% in parahippocampal tissues and 21.3% in entorhinal tissues [30]. While cell death in late stages of AD can certainly decrease tissue volume, the cause of this atrophy in living patients has not been positively ascribed to cell death and could arise from the collapse scenario presented here, although we caution that the spatial resolution is certainly insufficient for positive attribution. Axonal dystrophy is a hallmark of later stage Alzhimer's disease [31].
The collapse, if verified, has implications for treatment strategies for Alzheimer's. First, it is more likely to occur than MT dynamic instability, and thus represents the first irreversible signpost in late stage Alzheimer's. Second, treatments which aim to stabilize microtubules, such as those which are taxol based [32], will not inhibit this collapse. Rather, assuming as seems likely that phosphorylation is the source of tau removal from MTBs, approaches which inhibit activated kinases would appear more promising [33].
We plan future studies on combined mechanics/kinetics simulations of the MTBs under conditions of phosphorylation to attempt to determine the time scale for the collapse (which is not addressed here). We also will examine fully three dimensional models for the MTBs which include MT bending to provide a more realistic set of data for experimental tests on cultured cells.

Percolation and Equilibration Algorithm
The percolation simulation algorithm is developed from an algorithm used by Tang and Thorpe in their simulation of rigidity percolation of an elastic triangular net [20]. Tang and Thorpe's simulation analyzed the percolation of a triangular net under a constant tension, as if stretched onto a frame. In this case, each node is treated as a lattice point attached to six stretched Hooke springs. In our case, however, each node is attached to six compressed Hooke springs, and subject to a distance-dependent attractive depletion force. Therefore, Tang and Thorpes algorithm requires some modification in order to be applicable to our case.
Our simulation begins by populating the center points of the nodes and the tau springs between them. As they are populated, an index is assigned to each node and each spring connection site. Then, a random number generator chooses a (random) spring index and removes one spring. This loss of a spring perturbs the balanced forces in the system and destroys the mechanical equilibrium. The simulation enters a loop, which runs over all node position indices. Beginning with i~0, the simulation finds the net force vector acting on node i by summing the force vectors from the spring forces, depletion forces, and steric forces, viz. F F i~F F spring,i zF F D,i zF F steric,i : We then update the positionr r i (t) to the new fictive time tzdt by adding a distance proportional to the net forcẽ where a is an adjustable parameter used to control the step size. We use a~10 {6 nm/pN. This fictive molecular dynamics assumes the springs are overdamped or critically damped so that the nodes do not oscillate as they relocate. Once the last node is moved, we check to see if the system has equilibrated to within a tolerance of 10 {6 nm at each position during the previous loop; if not, we rerun the force relaxation loop until convergence at this tolerance is achieved. We average our results over 10-15 different starting random number seeds.
To estimate the second order elastic constant C at a given p value, we uniformly stretch or compress the area by a small dimensionless strain E~Da=a where a is the average bond length, reequilibrate at the new strain value, and take C from where E is the total potential energy. We take E sufficiently small that we are in the quadratic in strain region.