Nanoindentation of 35 Virus Capsids in a Molecular Model: Relating Mechanical Properties to Structure

A coarse-grained model is used to study the mechanical response of 35 virus capsids of symmetries T = 1, T = 2, T = 3, pseudo T = 3, T = 4, and T = 7. The model is based on the native structure of the proteins that constitute the capsids and is described in terms of the C atoms associated with each amino acid. The number of these atoms ranges between 8 460 (for SPMV – satellite panicum mosaic virus) and 135 780 (for NBV – nudaureli virus). Nanoindentation by a broad AFM tip is modeled as compression between two planes: either both flat or one flat and one curved. Plots of the compressive force versus plate separation show a variety of behaviors, but in each case there is an elastic region which extends to a characteristic force . Crossing results in a drop in the force and irreversible damage. Across the 35 capsids studied, both and the elastic stiffness are observed to vary by a factor of 20. The changes in mechanical properties do not correlate simply with virus size or symmetry. There is a strong connection to the mean coordination number , defined as the mean number of interactions to neighboring amino acids. The Young's modulus for thin shell capsids rises roughly quadratically with , where 6 is the minimum coordination for elastic stability in three dimensions.


Introduction
Simple globular viruses protect their strands of RNA or DNA with remarkable self-assembled proteinic shells known as capsids. The chemical and thermal stability of capsids has been studied for decades, but their mechanical properties are only beginning to be determined using nanoindentation [1,2,3]. There is little basic understanding of how mechanical strength varies between viruses, how it affects function, and how it is related to virus structure. In this paper we use a coarse-grained structure-based model to explore the types of mechanical response that capsids may exhibit and their relation to capsid geometry and protein bonds.
The capsids of simple globular viruses are typically of icosahedral symmetry and they are assembled from one or several kinds of proteins. The proteins cluster into n-meric capsomeres that form morphological units of capsids (so the proteins act as subunits). The structure of capsids has been explained by Caspar and Klug [4] as resulting from a regular triangulation of a sphere and is thus governed by the triangulation number T such that the number of subunits is equal to 60T. A pseudo-T = 3 virus has the symmetry of a T = 3 virus, but in which either the number of subunits in a capsomere is larger than in the standard classification, or the subunits are not sequentially identical. The size of capsids tends to grow with T but the actual size also depends on the size of the subunits. One question will be whether mechanical properties depend on T or on the nature of connectivity in the network of interactions between amino acids.
The nanoindentation technique typically involves anchoring a capsid on a substrate and then pushing on it with the tip of an atomic force microscope [1][2][3]. The tip typically has a larger radius than the virus and can be made nearly flat over the virus. The force, F , initially grows linearly as the plate/tip separation decreases. The slope corresponds to an effective spring constant, k. When a characteristic force, F c , is reached, the capsid undergoes a sudden partial collapse. Beyond this point the deformation process becomes irreversible and trajectory-dependent. The values of F c and k depend on whether the capsid is full or empty, i.e. if the genetic content of the virus is removed. They also depend, to a lesser extent, on the rate of squeezing. Among the five capsids listed above, k varies between 0.09 N/m for HBV and 0.57 N/m for MVM when the capsids are empty. The characteristic forces vary between 0.6 nN (CCMV) and 1.2 nN (MVM).
Here, we extend our previous theoretical studies [16] of nanoindentation of CCMV and CPMV (cowpea mosaic virus) to 33 other capsids for which the native structure is known. The capsids studied, together with their names and PDB codes, are listed in Tables 1 and 2. Also shown are the number of amino acids N, the mean radius R R, the capsid width dR, defined as twice the rms variation in radius, and the mechanical properties. Note that the values of F c and k vary by more than a factor of 20.
The model we use has been developed and tested in the context of single protein manipulation [18][19][20][21][22]. The model is coarsegrained, structure-based, and it comes with an implicit solvent. These features allow for studies of much larger proteinic objects than more detailed atomic models. We study capsids with up to 135 780 amino acids (NBV) and with radii up to 284 Å (HK97). Despite these simplifications, our model is molecular in nature which makes it distinct from the elastic shell model considered by Gibbons and Klug [23,24]. While shell models with appropriate constitutive laws can reproduce the initial portion of experimental force curves, they give smoothly varying capsid deformations.
They do not capture the intrinsically heterogeneous nature of the non-covalent bonds that dominate the mechanical response of capsids or allow for breaking of these bonds. Coarse-grained simulations of CCMV show deformation is localized at the boundaries of capsomeres and that F c is associated with breaking of inter-protein bonds [16]. Moreover, the rate and temperature dependence of F c show that bond breaking is thermally activated.
Another striking prediction of the coarse-grained model is a large difference in the mechanical properties of CCMV and CPMV even though they have similar N, R R and dR. F c is about three times bigger for CPMV than for CCMV whereas k is bigger by an order of magnitude [16]. This would require very different elastic properties to be used in an elastic shell model. We found that the difference correlated with differences in the average coordination number in the native state SzT, including connections along the backbone and to nearby amino acids. The average coordination number is only 6.4 for CCMV and rises to 7.4 for CPMV. Here, by examining a wider variety of capsids, we are able to show that SzT is indeed an important factor in determining the maximum stiffness of capsids of a given R R and dR. The results for k can be scaled to give an effective Young's modulus E using the Figure 1. Nanoindentation of model for CCMV. The figure contains four panels, each consisting of two subpanels. The lower subpanel shows the variation of force with separation and the upper shows c -the fraction of unbroken native contacts for one simulation at the indicated parameters. The top two panels show the temperature dependence of F (s) for flat plates, with numbers indicating the value of k B T=e. The left shows three trajectories for k B T=e~0:3 from Ref. [16] and one for zero temperature that is uniformly higher. The right repeats one F (s) trajectory for k B T=e~0:3 and shows two trajectories for each other temperature. The value of c is plotted for one case at the high and low temperatures. The two lower panels show the role of curvature in one of the indenting surfaces. Here, k B T=e = 0.3 and numbers indicate the radius of the tip R in nm. The lower left panel compares single trajectories obtained for three values of R to the flat results above. The lower right panel shows three trajectories for R = 30 nm. doi:10.1371/journal.pone.0063640.g001 formula for thin shells. This local property shows a quadratic trend with SzT{6, where 6 is the minimum coordination number for mechanical stability in three dimensions.
In the following section we describe the model used for the simulations. We then present results for CCMV to illustrate the role of temperature and a finite radius, R, of the tip. A finite R has little effect on k but reduces F c . Next, results for five systems are compared with existing experimental data. The model captures quantitative trends and shows capsids can exhibit a variety of qualitatively different force-distance curves. The next section compares results for all 35 virus capsids and the final section presents a summary and conclusions.

Methods
We start by presenting the pertinent aspects of the model. The molecular dynamics [25] model we use is exactly the same as in ref. [16]. In particular, native contacts between C a atoms are identified based on the overlap of the (slightly enlarged) van der Waals spheres associated with the heavy atoms of the amino acids [26]. Only capsids whose native state is available [27,5] are considered. It should be noted that in most cases the positions of some amino acids are not determined. These are not part of the fixed capsid structure and do not scatter coherently. These segments are believed to dangle inside the capsid in a disordered configuration and may be important for encapsulating RNA or DNA [28,29]. Since they do not stay in fixed positions, we assume that they do not contribute to the mechanical stability of the capsid.
The interaction for the native contact between atoms i and j at distance r ij is described by the Lennard-Jones potential , where s ij is determined for each pair ij so that the potential minimum coincides with the experimentally determined native distance. A native contact is considered broken if r ij exceeds 1.5 s ij . Interactions between atoms that are not part of a native contact are purely repulsive and are given by a Lennard-Jones potential with length s 0 that is truncated at the position of the energy minimum 2 1=6 s 0~4 Å . Covalent couplings along the protein backbones are described with a harmonic potential with spring constant 50 e=Å 2 . This value is high enough that covalent bonds are effectively rigid during capsid deformation but small enough that the equations of motion can be integrated without reducing the time step significantly.
Simulations are performed at a constant temperature using a Langevin thermostat [30] that mimics the effect of the surrounding solvent. The Langevin damping is large enough that the system is  Table 1. Characteristics of the T = 1, T = 2, and T = 3 virus capsids that are studied in this paper. The first three columns show the acronym used, the PDB structure code and the common name. The next four columns give geometrical parameters from the PDB structure: the number of C a atoms describing the model capsid, their average coordination number vzw and average radius R R and the thickness of the capsid shell dR (defined as twice the rms variation in the radial direction). The final three columns give mechanical properties from simulations at k B T~0:3e: the initial spring constant k, the force at the onset of irreversibility F c , and the effective elastic modulus E obtained using Equation in the overdamped limit where inertia can be ignored. Based on mapping simulations to experimental measurements of dynamic quantities such as diffusion, the time unit t in our simulation corresponds to *1ns [31,32]. Except as noted, simulations are at temperature k B T~0:3 e, which was used in previous studies of capsids [16] and many mechanical studies of proteins [22,18,19]. At this temperature, most of the proteins have good folding properties within the coarse-grained model [18]. In comparing to experiments we will use our most recent estimate of the binding energy parameter e*110 pN/Å , which is based on comparison to forces from experiments on protein stretching [22]. This energy corresponds to about 800 K, which would imply room temperature is closer to 0.35 e=k B than 0.3 e=k B . This difference is comparable to the uncertainty in the estimate of e and we will show that mechanical properties generally vary slowly with temperature.
To model nanoindentation the capsid is placed between two repulsive plates. The interaction with flat plates is described by a repulsive potential that scales as z {10 0 [33], where z 0 is the distance between the plate and a C a atom. A curved plate is described by the truncated and shifted Lennard-Jones potential of R is 30 nm. Whenever we refer to the ''curved'' case of nanoindentation, we mean a situation in which one plate is flat and the other is curved. Flat plates are oriented normal to the z{ axis and capsids are oriented so that this axis coincides with the zaxis in the structure file [27,5]. This is the 2-fold icosahedral axis. We have considered two additional sqeezing directions for CCMV and CPMV in ref. [16]. There was some change in forceseparation curves with capsid orientation, but the values of k and F c did not change significantly.
In the initial state, the plates are far enough apart that they do not interact with the capsid. The plates are then brought together by increasing the speed of both plates symmetrically to a value v p over 2000 t. In most cases considered here, v p~0 :0025 Å =t which corresponds to a combined speed of 0.005 Å /t. This was found to be slow enough to produce quasistatic results in studies of protein stretching [20,21]. While this velocity (*500 mm/s) is higher than experimental velocities (0.1 to 10 mm/s), it is slow enough to allow stress to equilibrate across the capsid. This is monitored by checking that the forces the capsid exerts on the two plates are equal and opposite when averaged over a small range of separations. The magnitudes of these forces are averaged to obtain the total compressive force, F , that would be measured by an AFM. This force is studied as a function of the separation s between the plates. For the case of a curved surface the separation is measured between the closest points.

Results
Mechanostability of CCMV -dependence on the temperature and on the curvature of the plate We now illustrate the type and magnitude of changes produced by temperature and R by considering the CCMV capsid. Figure 1 shows examples of nanoindentation trajectories for this capsid under various conditions. Note that experimental data are typically plotted vs. the displacement of the tip towards the substrate, which becomes more positive as s decreases. In addition, raw AFM data needs to be corrected for the compliance of the cantilever, while there is no compliance in our simulations.
In the top left panel there are three trajectories F (s) obtained at k B T~0:3 e that are reproduced from ref. [16]. They, and especially the central highlighted one, provide a reference for all of the results presented here. The three trajectories nearly coincide at the beginning of the sqeezing process. In particular, there is no force for sw284 Å . The force then rises linearly with slope k to a characteristic force F c as s decreases to about 175 Å . This corresponds to compression by about 35% from the force onset which is comparable to what is observed in experiments on CCMV [10,11]. The value of F c is of order 5 e/Å , i.e. about 550 pN, which is also close to the experimental value of 600 pN [10,11]. In the linear regime, the indentation process depends only weakly on the squeezing speed and is reversible: reversing the velocity of the plates results in approximate retracing of the F (s) curves. When the separation is reduced beyond the linear regime, there is a sudden drop in the force and then a steep rise at small separations due to steric repulsion. Once the force has dropped, the deformation is not reversible. If the plates are retracted, the force falls rapidly to zero [16].
The onset of irreversibility in our simulations of CCMV at k B T~0:3e is associated with rupture of native contacts [16], rather than a buckling transition like that observed in elastic shell models [24,10,23]. The top subpanel of the top left panel in Figure 1 shows the fraction, c, of unbroken native contacts as a function of s. In the initial stage of compression, the number of unbroken contacts remains essentially equal to that in the native structure, although thermal fluctuations are strong enough for a few contacts to fluctuate between broken and unbroken states at k B T~0:3 e. At the end of the linear regime, nearly 50% of the contacts rupture as the force drops rapidly. Most of these bonds connect different proteins within the capsid, rather than different parts of the same protein [16]. These bonds do not reform on the time scale of our simulations, but may reform in experiments on the time scale of many minutes [3].
Several aspects of the results show that the bond rupture at F c is thermally activated. The top right panel of Figure 1 shows the temperature dependence of F (s). For all k B Tw0, the value of F c depends on trajectory. While the stiffness decreases only weakly as k B T increases, F c drops significantly and is almost completely suppressed at k B Tw0:4 e for CCMV. These results are consistent with thermal fluctuations being able to activate the transition at lower F as k B T increases and to produce run-to-run variations in F (s) trajectories. Our previous studies showed that reducing the speed of compression, and thus reducing the time for activation, raises F c [16].
For all nonzero temperatures shown, we find a sharp drop in c as the force drops below F c . However at k B T~0 there is a drop in F that is not associated with bond breaking and occurs at very small separations. Since there are no thermal fluctuations, the drop occurs at a well-defined separation. This case appears to reflect a buckling instability of the capsid like that predicted for elastic shell models [24,10,23].
We now consider the influence of tip curvature on capsid response. The lower left panel of Figure 1

Nanoindentation for capsids considered experimentally
As mentioned in the Introduction, nanoindentation measurements have also been made for MVM, NV, HBV, and HK97. Figure 2 shows F (s) traces for coarse-grained models of these capsids at k B T=e~0:3 and 0 with flat and curved walls. All lines show a rapid upward curvature in the first 5-10 Å after contact. In this range the repulsive potentials from the confining walls are just beginning to overlap with the outer atoms and it is not included in fitting the spring constant k. At smaller separations the traces illustrate the four types of behavior found for the full range of capsids studied later. These are a single peak at F c (MVM), a peak with a shoulder (NV), a series of gradually rising peaks before F c (HBV), and a peak followed by a long plateau (HK97). Note that curved walls do not change the type of force trace. As seen in the previous section, there is just a slight reduction in k and F c because the compressive stress is focussed on a smaller region. The magnitude of this shift is largest and most consistent between different runs for HK97 and HBV.
The traces for MVM and NV are relatively simple. At k B T~0:3e there is a single main peak followed by a sharp drop where about half the native contacts break. The main difference is that some traces for NV show a shoulder before F c and others    Table 1 break at the same separation (*330Å ). At k B T=e~0, the force peak is replaced by an inflection point and no bonds break. This indicates that the capsids are always mechanically stable against small perturbations, but can fail through thermally activated bond breaking at high enough temperatures and long enough times.
A series of peaks is superimposed on a steadily increasing force in the case of HBV. This complicates the determination of F c . For k B T=e~0:3, the majority of bond breaking occurs after the final peak at about 15 e=Å . A sufficient number break near s~180Å and F(s)*6 e=Å , to hinder reversibility and we take this as F c . In analyzing experimental data, the first sharp peak at *225 Å and F (s)*5 e= Å might be reported as F c even if retraction would have produced a relatively reproducible force trace. Thus there are greater systematic uncertainties in F c for capsids in Tables 1 and 2 that exhibit this type of force trace. The force trace at zero temperature shows peaks at very similar s with no bond breaking. This indicates that the capsid loses mechanical instability even in the absence of thermally activated bond breaking and the nature of these instabilities will be the focus of future studies.
Finally, in the case of HK97, there is a weakly articulated force peak followed by a plateau. The plateau ends with a drop for k B T=e~0:3 and flat walls, and extends to the onset of steric repulsion for curved walls and at zero temperature. Significant bond breaking does not occur till sv400 A which is far along the plateau. In addition, the first peak occurs near s~450 Å for all temperatures. Once again, this is indicative of a mechanical instability like the buckling instability seen for thin shells. Note that k is particularly difficult to determine for HK97. The initial slope for forces up to about 2e=Å is used in the table. At lower separations the flat wall results rise more steeply, while the tip results do not. Results for HK97 may be more sensitive to tip size than other capsids because it is the largest and has a radius that is more than double most other capsids. As the force grows, the tip may push through the capsid like a needle. A more extreme case was considered in Ref. [34] where the effective tip radius was only a few Å and it passed completely through the capsid shell. Figure 3 compares theoretical results for F c and k with experimental findings. Note the clear linear correlation between theoretical and experimental values. Indeed, using the value of e determined from matching forces from the model to protein stretching experiments [22] gives forces that are quantitatively similar to capsid experiments (dotted lines). Note that a variety of force curves are obtained for MVM. The highest and most pronounced peak in Fig. 2c of ref. [9] is at 1.2 nN, which is consistent with the trend shown in the lower panel of Figure 3. Some capsids showed no instablities up to this load, while others showed precursor peaks. These variations and observations of capsid geometry. indicate that capsid orientation and thermal activation of bond breaking are important and that failure is localized between proteins making up different trimers. Failure of interprotein bonds is consistent with our earlier results on CCMV [16].
Experimental values of k are better fit by increasing the correspondng theoretical values by a factor of 2 (solid line). As explained in Ref. [16], the only length scale included in the coarsegrained model is the separation between C a atoms in native contacts and this determines both the structure and the rate of change of forces with separation. The actual rate of change of  forces will be determined by the shorter distances separating individual atoms of the amino acids associated with each C a . This will increase k relative to the coarse-grained prediction as seen in Fig. 3. We conclude that the coarse-grained model can be used to predict trends in mechanostability for empty virus capsids and provide quantitative estimates for F c using e=Å = 110 pN while k should be doubled for this normalization.

Nanoindentation processes in other virus capsids
We now present the F (s) curves for the remaining capsids listed in Tables 1 and 2. Their selection was based entirely on availability of structural information and they can be viewed as an essentially random sample. Figures 4 and 5 present results for capsids belonging to the T = 1, T = 2, and T = 4 structural classes. Figure 6 corresponds to the T = 3 class, Figures 7 and 8 to the pseudo T = 3 class, and Figure 9 to the T = 7 class. Most of the results have been obtained with flat walls and are shown by solid lines.
For some capsids, dotted lines show results for a tip with radius of curvature R~30 nm. As in the previous section, the curved tip does not change the type of force curve and produces a small decrease in mechanical strength. In most cases there is little change in k. The exception is HK97 (Fig. 2) where k drops about 15%. Changes in F c tend to be larger than those in k, but the largest change is a 25% drop for POLIO (polio virus -type I Mahoney strain). Given the large number of viruses considered, we expect that these results provide reasonable bounds for the magnitude of changes in mechanical properties that would be produced by changing the radius of an AFM tip from nearly flat to a typical value of 30 nm.
In most cases, the F (s) curves show a well defined single force peak that is not very sensitive to the choice of trajectory. Several show a change in slope before the peak and in a few cases this developes into a weak peak: STNV (T = 1), SPMV (T = 1), BmDNV (T = 1), and SBMV (T = 3). More and larger peaks are seen for PhiX (T = 1), PIC (T = 2), TYMV (T = 3), and TRSV (pseudo T = 3). A few force curves show a weak peak followed by a  Figure 7 for the six remaining pseudo T = 3 capsids listed in Table 2. doi:10.1371/journal.pone.0063640.g008 . The different types of behavior are not uniquely associated with specific values of T, but weak plateaus are more common for T = 3 capsids and all but one pseudo T = 3 capsid exhibits a single sharp peak. There is also no clear correlation with capsid size. Capsids with weak peaks and plateaus tend to have larger values of dR= R R than other capsids in the same structure class, but NV has a sharp peak and a thick capsid.
The structure in F(s) and variability between trajectories leads to systematic uncertainties in determining F c . When there is a single major peak or clear plateau, we take the average height over trajectories, including some that are not shown in the figures. When there is a shoulder where some trajectories show sharp force drops, we use the shoulder height. The choice is less clear for cases with multiple peaks and leads to significant uncertainty for HBV (T = 4), TYMV (T = 3), and TSRV (pseudo T = 3). The quoted values correspond to the height of the lowest significant peak, while later peaks are a factor of two higher. The highest value is for HRV with F c~3 2 e=Å corresponding to about 4 nN. The smallest values for BMV, CMV and IBDV are about 20 times smaller.
We now ask: what do the values of k and F c depend on? A natural attribute of the capsids to consider is their size as characterized by R R. The size dependence of k and F c is shown in Figures 10 and 11, respectively. There is clearly no simple relation of mechanical properties to R R. The largest capsid, HK97, has half the stiffness of the smallest capsid, SPMV, but bigger variations in k are found between capsids with the same radius. Another important factor might be the structural classification. The pseudo T = 3 capsids tend to be strongest. The weakest are T = 2, T = 7 and some of the T = 3 capsids. T = 4 capsids span the full range of strength, while T = 1 and T = 3 span intermediate values. These observations indicate that there is some correlation to structure, but not a strong one.
In our previous study of CCMV and CPMV we noted that while most of their geometrical properties were very similar they had very different average coordination numbers SzT. When counting non-bonding contacts and the two covalently bound neighbors along the backbone, SzT is 6.36 and 7.40 for CCMV and CPMV, respectively. The minimum number of neighbors required for stability in 3 dimensions is 6 if, as here, there are no frictional or bond angle forces. Studies of rigidity percolation indicate that the elastic modulus is a strong function of SzT near the onset of rigidity [35][36][37][38]. We argued that this could explain why increasing the number of native contacts per atom by only 45% (and total bonds per atom by 30%) could lead to an order of magnitude increase in stiffness for CPMV relative to CCMV.
To explore the effect of coordination number, the results for k and F c are replotted against SzT in Figures 12 and 13, respectively. For the T = 3 and pseudo T = 3 capsids there is a clear tendency for k and F c to grow with increasing SzT. There is no clear trend for the other capsids when viewed as a group. However, they have a much larger range of sizes than T = 3 and pseudo T = 3. The dotted lines in the lower panel show that there is a clear trend for mechanical strength to grow with SzT if capsids with similar radii are considered separately.
To try to separate the effect of capsid dimensions from the local modulus describing the mechanical properties within the shell, we follow a dimensional analysis motivated by the elastic shell models of Gibbons and Klug [23]. In the thin shell limit, the stiffness is proportional to the Young's modulus E and the square of the shell thickness and inversely proportional to the radius. Up to a numerical prefactor that would depend on how thickness is defined, one can use this scaling to define an effective modulus.
that characterizes the local response in the shell (Tables 1 and  2). Figure 14 shows how E varies with SzT for all capsids studied. There is a much clearer correlation between these quantities than found in the previous plots. The majority of the capsids show a roughly parabolic dependence on the excess above the minimum coordination number for rigidity.
where E 0 & 0.05 [e/Å 3 ]. This correlation is quite good given that there is no correction for the fact that SzT will be lower for the significant fraction of C a that lie on the outer and inner surfaces of the capsid, and that the distribution of bonds may be nonuniform. For example, reduced local coordination along the boundaries between capsomers could greatly lower the global stiffness k and explain why bond breaking may localize there [16]. It is also interesting to note that the greatest outliers (IBDV, PhiX, NBV and STNV) are those that are farthest from the thin shell limit, having dR= R R greater than 0.24. The thin shell formula is questionable in this limit and the loss of coordination number due to surfaces would be smaller than for other capsids.
Studies of randomly linked systems have also found a power law relation between modulus and the excess coordination. Square and Kagome lattices are marginally stable when only nearestneighbors interact. Their shear modulus grows as the square of the number of next-nearest-neighbor bonds [37], which is consistent with Fig. 14. A different linear scaling is found for completely random networks near the jamming transition [38]. Capsids may be closer to the orderly structure of lattices because of the backbone connectivity and repeated structure, but it is not clear that they should fit into either class.
It should be noted that the coordination number for a given capsid is sensitive to the details of the determination of the contact map. Changing the details would alter the number of native contacts and thus SzT. However, we expect that the trend for local elastic properties to rise rapidly with SzT{6 would remain. The overlap criterion used is based on previous studies of protein folding and stretching, where it has been shown to provide reasonable results [20,22].

Discussion
A coarse-grained molecular model was used to study the mechanical response of 35 virus capsids. The full force-separation curves have a variety of shapes, but in general share two common features. In particular, there is a linear elastic response characterized by a spring constant k at small deformations and a sharp drop or plateau at a characteristic force F c that signals an irreversible instability. As found previously for CCMV [16], F c is usually associated with bond breaking at finite temperature. Because bond-breaking is thermally activated, there are run-to-run fluctuations in F c , and F c decreases with increasing temperature and increases with increasing indentation rate. Similar effects have been seen in previous simulations of mechanical unfolding of proteins where our model captures the breaking of bonds during experiments on unfolding. There is rarely any bond-breaking in capsids at zero temperature, indicating that the bonds are always metastable. Some capsids show sharp instabilities without bond breaking at zero temperature and large forces that is indicative of a buckling instability like that seen in elastic shell models [23,24].
The elastic response and onset of irreversibility have also been considered by Arkhipov et al. [39,40] using an even more coarsegrained model based on capsid geometry where a single bead represents 150 protein atoms, on average. Bonds are defined based on proximity of these larger beads and bond-angle terms are introduced to include the effect on rigidity of the many additional atoms that have been removed. All interactions are then fit to atomistic simulations of capsids. No bond-breaking is possible in this model because the bonds are strictly harmonic, but capsids undergo instabilities like those seen here at zero temperature. The model has been used to study the native state stability of several T = 1 and T = 3 capsids (STMV, SPMV, STNV, BMV) [39] and to get the pre-collapse behavior of the force-indentation curves for the T = 4 capsid HBV [40]. It is interesting that even with the greater coarse-graining, their results for HBV show a sinuous character that is similar to our results.
For the capsids studied in this paper, the values of k and F c vary by about a factor of 20. These variations are not correlated with virus symmetry (T) or size. Indeed, nearly the full range of values is sampled by T = 3 and pseudo T = 3 capsids with radii between 130 and 134Å . The greatest correlation was found with the coordination number SzT that describes the number of bonds constraining motion. To isolate the effective local elasticity from geometrical effects, we determined an effective Young's modulus E for regions within the capsid shell using the thin elastic shell formula for stiffness (Eq. 1). On average there is a trend for E to rise quadratically with SzT{6, where 6 is the minimum coordination for stability in three dimensions. The largest deviations are for the thickest capsids and lie below the trend. This may reflect deviations from the thin shell scaling or local fluctuations in the coordination number that produce weak spots that dominate the response. Even the largest values of E are an order of magnitude smaller than an fcc lattice with SzT~12 and the same interactions. It is interesting to ask whether the greater flexibility of all capsids and the variations in E for specific capsids are important to function. This is outside the scope of the current paper, but one may speculate that viruses may be more rigid if they do not need to reform during their life cycle or are exposed to more extreme environments.
Experimental values of k and F c are only available for 5, viruses [3]. Our calculations reproduce the trends in these quantities and good quantitative agreement with experiments is obtained if the interaction strength is set to the value obtained from fits to protein stretching e=Å~110 pN. Calculated values of k are consistently about a factor of 2 too small with this interaction strength. As discussed, the current coarse-grained model assumes that the separation between C a bonds determines the rate of change of forces. Better quantitative agreement could be obtained if the variation in force was related to the smaller separation between amino acids.
The variability in elastic properties of virus capsids has also been observed by Tama and Brooks [41]. They also considered only C a atoms and assigned Hookean springs between nearest-neighbors. They then analyzed the normal modes of the system and correlated them to structural changes of swelling rather than nanoindentation. It would be interesting to determine the frequency of normal modes that correspond most closely to indentation using this method. This would also enable rapid studies of the effect of capsid orientation on k, which has been found to be small in previous studies [16,39,40].
The studies presented here have focused on the experimentally accessible macroscopic response of capsids. Future studies should assess the variability in local response within capsids. Variations in local deformation may be correlated with changes in local coordination number and/or with the boundaries of proteins as in our earlier simulations and recent experiments on MVM [9]. It may also be possible to relate them to local magnitude variations of the eigenmodes obtained by normal mode analysis [41]. These studies could help explain the variations in E at a given SzT and will be a useful stepping stone towards modeling still larger biological systems.