Soft modes near the buckling transition of icosahedral shells
Abstract
Icosahedral shells undergo a buckling transition as the ratio of Young’s modulus to bending stiffness increases. Strong bending stiffness favors smooth, nearly spherical shapes, while weak bending stiffness leads to a sharply faceted icosahedral shape. Based on the phonon spectrum of a simplified mass-and-spring model of the shell, we interpret the transition from smooth to faceted as a soft-mode transition. In contrast to the case of a disclinated planar network where the transition is sharply defined, the mean curvature of the sphere smooths the transitition. We define elastic susceptibilities as the response to forces applied at vertices, edges and faces of an icosahedron. At the soft-mode transition the vertex susceptibility is the largest, but as the shell becomes more faceted the edge and face susceptibilities greatly exceed the vertex susceptibility. Limiting behaviors of the susceptibilities are analyzed and related to the ridge-scaling behavior of elastic sheets. Our results apply to virus capsids, liposomes with crystalline order and other shell-like structures with icosahedral symmetry.
I Introduction
Virus capsids Caspar and Klug (1962) and other structures such as colloidosomes Dinsmore et al. (2002) and liposomes Spector et al. (1996); Delorme et al. (2006) consist of thin shells of spherical topology that frequently exhibit icosahedral symmetry. A popular simplified model Lidmar et al. (2003); Nguyen et al. (2005); Vliegenthart and Gompper (2006); Hicks and Henley (2006) replaces the shell with a triangulated network of masses and springs (see Fig. 1). This network consists of five- and six-coordinated vertices, with the five-coordinated vertices aligned with the five-fold icosahedral symmetry axes. Five-coordinated vertices may be considered as disclinations within an otherwise six-coordinated lattice. These disclinations are absent in conventional continuum models of spherical shells Timoshenko (1940); Love (1944); Niordson (1985).
Elastic properties of the capsid can be mimicked by suitably adjusting the spring constants to obtain the desired Young’s modulus and by imposing a curvature energy to obtain the bending modulus . Strains associated with the disclinations cause the network to buckle Seung and Nelson (1988), transforming the shape from smooth and nearly spherical to strongly faceted and nearly icosahedral Lidmar et al. (2003). A dimensionless parameter controls the transformation. We define the Foppl-von Karman number
(1) |
where is a linear dimension of the shell. The buckling occurs when exceeds a value of order (see Fig. 2). For the virus HK97, which appears to facet as it matures Wikoff et al. (2006); Ivanovska et al. (2004), reaches a value of order according to the estimate of Ref. Lidmar et al. (2003). Varying the pH of solution can alter , with the range 100-900 reported for the virus CCMV Klug et al. (2006); Tama and Brooks (2002). At much larger values of (in excess of ) which should characterize liposomes with crystalline order, an interesting phenomenon known as “ridge scaling” emerges Witten and Li (1993); Lobkovsky et al. (1995); DiDonna and Witten (2001); DiDonna (2002); Wood (2002); Witten (2007).
Caspar and Klug Caspar and Klug (1962) classify icosahedral structures by a pair of integers . A pair of five-coordinated vertices is connected by a path consisting of edges in some given direction and edges in a direction to the left (e.g. between two blue vertices via a red vertex in Fig. 1). The -number of the network, gives the number of vertices as . There are always 12 five-coordinated vertices, so the number of six-coordinated vertices is . Structures with and both nonzero and are chiral, such that and are mirror images. Their symmetry group is the 60-element icosahedral rotation group . Structures with either or , or with are nonchiral. Their symmetry belongs to the 120-element group , which should not be confused with the 120-element icosahedral double group Widom (1986) .
We exploit the rotational symmetry group to analyze the normal modes of the network model by diagonalizing the Hessian matrix of the elastic energy. Eigenvectors represent characteristic modes of deformation, which transform according to irreducible representations of , and the corresponding eigenvalues measure the mechanical stability. Because the buckling occurs in a symmetric fashion, the corresponding modes must exhibit full icosahedral symmetry. Nondegenerate modes transform as the unit representation. Tracking these nondegenerate eigenvalues reveals a softening and also a mixing of modes as passes through .
Other studies consider more microscopis elastic network models Rader et al. (2005); Tama and Brooks (2005) that place nodes at every atom in the amino acid chains. These studies find that the displacements during maturation (i.e. as the virus goes through the buckling transition) can be accurately represented using a superposition of only the lowest few nondegenerate modes, consistent with our expectations.
Section II of this paper reviews the continuum-elastic theory for deformations of planes and spheres, to establish notation and for comparison with our later numerical results. Our network model is defined in section III and applied to the special cases of disclination-free triangular lattices, single disclinations of positive charge, and icosahedral structures of spherical topology containing twelve disclinations. Low-lying eigenvalue spectra reveal a sharp buckling transition in the case of a single disclination but a broadly smeared transition for the icosahedral case. Following Ref. Lidmar et al. (2003) we find that the positive curvature of the sphere plays a symmetry-breaking role analagous to an applied magnetic field at a ferromagnetic phase transition.
The final section (IV) applies forces to selected points on a plane or a shell to probe the elastic response of the network as a whole. The resulting displacements define susceptibilities which diverge in the case of the single disclination. In the case of the icosahedron, we find the effective stiffness (inverse of the susceptibility) drops most rapidly at for forces applied at five-fold symmetry axes, but the stiffness falls off more rapidly for forces applied at two- and three-fold symmetry axes for . We analyze these susceptibilities in limiting cases of small and large .
Ii Continuum-elastic theory
The general elasticity theory of membranes can be expressed in coordinate-free form Peterson (1984); Niordson (1985). Let be a manifold (a two-dimensional smooth surface embedded in three dimensional space) assumed to be in mechanical equilibrium. Now impose tangential deformation and normal deformation corresponding to displacements of points on the surface. Let and be the metric and curvature tensors respectively of after distortion. Greek indices take values 1 and 2 corresponding to the dimensions of . Define the strain tensor
(2) |
where and indicates covariant differentiation with respect to . The trace measures dilation, while
(3) |
measures shear strain. Bending of is characterized by mean curvature and Gaussian curvature .
The free energy density at contains dilation, shear and bending contributions,
(4) |
and can be integrated over to obtain the total free energy
(5) |
The separate contributions are
(6) | |||
The elastic constants and are the Lame constants Landau and Lifshitz (1986). The 2D area (bulk) modulus , while itself is the shear modulus, and the 2D Young’s modulus . Upon integration over the surface , the Gaussian curvature term becomes constant, and we neglect this term henceforth. Likewise, we set the spontaneous curvature , thus assuming the manifold would be flat in the absence of constraints associated with the spherical topology. Effects of are discussed in Ref. Nguyen et al. (2005).
Given the elastic free energy we obtain the stress tensor
(7) |
whose divergence yields the tangential force
(8) |
The normal force is given by
(9) |
In mechanical equilibrium the stress tensor and normal force vanish. Slightly out of equilibrium, to first order in the displacements, special forms of and known as normal modes solve the eigenvalue equation
(10) |
When displaced from equilibrium according to the normal mode , the free energy increases by
(11) |
According to the equipartition theorem the modes fluctuate with thermal energy and amplitude .
Time dependence of the strains depends on the equations of motion. In the overdamped case we write
(12) |
where we take proportional to an inverse viscosity as in a Stokes-Einstein relation. In this case a normal mode decays in time with a decay rate . Ref Levine and MacIntosh (2002) carries out a more thorough investigation of flat membranes coupled to fluid flow. In the absence of damping we write
(13) |
with the 2D mass density. A normal mode now oscillates in time at frequency .
ii.1 Deformations of a Plane
An infinite flat elastic sheet in equilibrium has no curvature, so for small perturbations the energy decouples into contributions from the in-plane strain and perpendicular displacement .
(14) |
Here is the usual 2D laplacian operator and the usual gradient. By differentiating the energy we obtain the forces
(15) |
and
(16) |
Because the in-plane and out-of-plane displacements decouple, we solve them separately. The solutions are based on the plane wave function
(17) |
which is an eigenfunction of the Laplacian operator . In-plane normal modes are expressed as longitudinal waves
(18) |
and transverse waves
(19) |
Note the identities and , as expected for longitudinal and transverse waves. These waves are eigenvectors of the in-plane force eq. (15) provided their eigenvalues obey the longitudinal and transverse dispersion relations, respectively
(20) |
and
(21) |
Perpendicular out-of-plane waves
(22) |
obey eq. (16) subject to the perpendicular wave dispersion relation
(23) |
For future reference we recast the normal modes in plane-polar coordinates , replacing the plane-wave function with cylindrical Bessel functions
(24) |
The Laplacian operator takes the form
(25) |
Like the plane wave function , waves of type (24) are eigenfunctions of the Laplacian operator, . Upon defining normal modes as in eqs. (18), (19) and (22) the longitudinal, transverse and perpendicular dispersion relations given in eq. (20), (21) and (23) result. These polar forms generalize nicely to conical and spherical geometries.
ii.2 Deformations of a Sphere
Now we redo the prior calculation of section II.1 for the case of small perturbations around a sphere of equilibrium radius . In this case the unperturbed manifold has constant mean curvature . Consequently the free energy acquires a term coupling the in-plane and normal strains through the dilation energy.
(26) | |||
In the above, the Laplacian operator takes the form
(27) |
Notice we include the integration measure along with the bending energy , because it contributes the term . The factor is not needed in or because these are already second order in the deformation.
Taking functional derivatives of yields the stress tensor, in-plane and normal force
(28) | |||
where we define
(29) |
The extra in eq. (28) for comes from commutation of covariant derivatives. The final, second derivative, term in (29) comes from integrating by parts the square of the first derivative in .
Take the spherical harmonic as the basic deformation, analagous to the plane wave in eq. (17) or the cylindrical wave in eq. (24). The spherical harmonic is an eigenfunction of with eigenvalue and an eigenfunction of with eigenvalue . By analogy with the procedure for plane waves in flat space, we take derivatives as
(30) | |||
where is the alternating tensor. We also define
(31) |
These functions are linear combinations of the “Vector Spherical Harmonics” , and , which form a complete set of orthogonal functions for expanding vector fields on the surface of a sphere Hill (1954); Widom (1986). Notice that the transverse mode is proportional to the angular momentum operator acting on , thus identifying it with the vector spherical harmonic . The longitudinal and perpendicular modes, and , are linear combinations of and . Note that the longitudinal and transverse modes and exist only for , while exists for .
The transverse mode is divergenceless () and hence creates no perpendicular force and no longitudinal force (the gradient part of ). In fact, it is an eigenfunction of the force (28). Upon taking into account the commutation of covariant derivatives, we find from which we obtain the eigenvalue
(32) |
As expected, for because these modes correspond to rigid rotations.
In contrast to the transverse modes, the longitudinal and perpendicular modes and are coupled in both the tangential force and perpendicular force . In matrix form,
(33) |
Setting as in eq. (30) and setting as the radial component of in eq. (31), the matrix elements of become
(34) |
The eigenvalues of this matrix, , are the desired normal mode eigenvalues . For the special case the eigenvalues are and . The vanishing eigenvalue corresponds to rigid translation (for example, the north and south pole displace upwards perpendicular to the shell while the the equator displaces upwards tangent to the shell). The finite eigenvalue corresponds to an “optical” mode in which polar and equatorial regions displace in opposite directions (for example, the north and south poles displace upwards while the equator displaces downwards).
The spherical solution should go smoothly to the flat space solution in polar coordinates as the sphere radius . This correspondence can be verified by holding , and fixed, and noting Abromowitz and Stegun (1970)
(35) |
In addition, the eigenvalues should approach their proper limits. Clearly approaches its flat space limit (21). To check , note that the eigenvalues of the matrix (34) approach and in the limit of large sphere radius , yielding the flat space limits Eqs. (20) and (23).
Iii Mass and spring model
We now introduce the discrete mass and spring model for which numerical calculations will be performed. This model is also closer to reality for liposomes and colloidosomes, which consist respectively, of discrete lipid molecules and colloidal particles, and also for viruses, which consist of an aggregation of discrete protein subunits known as capsomers. Let be the position of mass and be the orientation of plaquette . A plaquette is a set of three masses joined to each other by springs, and we take the normal in the outward direction. Following Lidmar et al. (2003) we define
(36) |
and
(37) |
and set the unstretched spring length . Here denote pairs of nearest-neighbor vertices, and denote pairs of adjacent (edge-sharing) plaquettes. In the continuum limit the discrete model reproduces the continuum system with elastic constants
(38) |
Foppl-von Karman number
(39) |
Lame coefficients and bulk modulus
(40) |
and 2D mass density (taking the vertex mass )
(41) |
iii.1 Deformations from flat-space
Consider a regular six-coordinated triangulated network of masses and springs. As before we start with the plane wave function (17) and take its gradient to obtain the longitudinal sound wave. The dispersion relation is simplest for wavevector in the direction (chosen to lie midway between two near-neighbor bonds)
(42) |
Taking the cross product with yields the transverse sound wave with dispersion relation
(43) |
Finally, taking the perpendicular displacements as the planewave yields the perpendicular modes with dispersion relation
(44) |
In the continuum limit these dispersion relations revert to the prior results of continuum elastic theory.
iii.2 Buckling of a Plane into a Cone
Upon introducing a five-fold disclination into the flat triangulated network discussed previously, strain energy accumulates Seung and Nelson (1988) and grows without bound as the radius of the network increases. At some specific “buckling radius” it becomes energetically favorable to buckle out of plane, trading a reduction in strain energy for a cost in bending energy. The trade-off is measured by the Foppl-von Karman number . Small favors flat networks, while larger favors buckling into a conical shape.
In the following we analyze the vibrational spectrum of the network as it passes from flat to conical. Rather than vary the radius, we hold fixed and vary the bending stiffness. Large opposes buckling and the network lies flat, while small allows buckling out of plane into a cone. For the network of radius analyzed below, buckling occurs for . As increases the threshold value of grows as so that approaches the limiting value Seung and Nelson (1988); Lidmar et al. (2003).
Eigenvectors of the Hessian matrix form basis functions for representations of the symmetry group of the structure Tinkham (1964). Eigenvectors sharing a common eigenvalue form the basis for an irreducible representation. Thus the patterns of degeneracy follow the dimensionalities of the irreducible representations, as can be seen in Fig. 3a. Likewise the eigenvectors exhibit special symmetry properties associated with subgroups of the full symmetry group.
The symmetry point group of the cone is in general, corresponding to five-fold rotations around an axis passing through the five-coordinated vertex, together with reflections in vertical planes passing through this axis (see Table 1). For the specific case of the flat network, the group is even higher, , adding reflections in the horizontal plane, and two-fold rotations around axes lying within the plane. For both groups all irreducible representations are either 1- or 2-dimensional, so all nonzero eigenvalues must be nondegenerate or two-fold degenerate. Of course, there must be a sixfold degeneracy of zero eigenvalues, corresponding to rigid translations and continuous rotations (not belonging to the finite point group) that leave the energy invariant.
5 | |||||
---|---|---|---|---|---|
0 | 1 | 1 | 1 | 1 | |
0 | 1 | 1 | 1 | -1 | |
1 | 2 | 0 | |||
2 | 2 | 0 |
For the group , the irreducible representations are based on those of supplemented with an additional label according to whether they are even or odd under reflection through the horizontal plane . The requirement that each irreducible representation be either even or odd under requires that each mode be polarized either fully in-plane or fully perpendicular.
Let be the lowest nonzero eigenvalue. Its eigenvector is polarized strictly perpendicular to the sheet and transforms as the irreducible representation . Its value is nonzero at the origin. The energy of mode varies as where measures the amplitude of the mode. Mechanical equilibrium thus demands that all eigenvalues (other than the six zero modes) be strictly positive. In particular it requires . However, if we monitor the value of as a function of (Fig. 3a) we find it crosses through zero at .
For small deformations we express the energy as
(45) |
Now set . The mechanically stable minimum energy structure is perfectly flat () for , but it buckles out of plane for , in a shape described by the eigenvector , with amplitude growing as . Meanwhile the energy drops as . These effects can be seen in figure 3b.
For , figure 3a shows the spectrum of vibrations around the mechanically stable, buckled structure. Note that (the lowest nondegenerate eigenvalue) becomes positive again.
iii.3 Buckling of Spherical Shells
iii.3.1 icosahedron
Table 2 presents the character table of the 60-element icosahedral rotational symmetry group , which has 5 irreducible representations. The conjugacy classes are labeled , where is the order of elements in the class, so that an element of corresponds to a rotation by . Recall that the spherical harmonics form basis functions for the irreducible representations of the continuous rotation group , and therefore they induce representations (in general reducible) of . For a given angular momentum and rotation angle , the character is
(46) |
Irreducible representations of are labeled in Table 2 according to the lowest angular momentum under which they transform. Of particular interest is the representation corresponding to angular momentum . This is the representation under which three-dimensional vectors transform.
0 | 1 | 1 | 1 | 1 | 1 | |
---|---|---|---|---|---|---|
1 | 3 | -1 | 0 | |||
(3) | 3 | -1 | 0 | |||
(3) | 4 | 0 | 1 | -1 | -1 | |
2 | 5 | 1 | -1 | 0 | 0 | |
12 | 0 | 0 | 2 | 2 | ||
36 | 0 | 0 |
The simple icosahedron has 12 vertices, 20 faces and 30 edges. Since we place masses on the vertices, our eigenstates are functions defined only at vertex positions. Arbitrary scalar-valued functions can be expressed as linear combinations of the basis functions of the “regular representation”, one of which is concentrated at each icosahedron vertex. The characters of the regular representation equal the number of vertices that remain stationary under a given symmetry operation. Our vibrational modes are vector-valued functions on the set of vertices and thus can be expressed as linear combinations of the product of the regular representation times the representation corresponding to a three dimensional vector. We call the resulting product representation the “total vibrational representation” Widom (1986), and its characters .
Reducible representations can be decomposed into their irreducible components using orthogonality properties of character tables. In particular we obtain the decomposition
(47) |
Of the three occurrences of the vector representation we know that two must correspond to rigid global translations and rotations. These leave the energy invariant and hence are zero frequency modes. The nondegenerate mode transforming as the unit representation must correspond to a “breathing mode” in which all vertices displace equally in the radial direction. We find that the remaining modes have specific interpretations in terms of vector spherical harmonics, as listed in table 3.
Formula | Irrep | Comments | ||
---|---|---|---|---|
0.00000 | 0 | 6 | Translations + rotations | |
0.58579 | 5 | Mixed contains and | ||
0.76393 | 3 | Radial contains | ||
1.00000 | 1 | 5 | Tangent | |
1.80901 | 4 | Tangent | ||
2.76393 | 1 | Radial breathing mode | ||
3.00000 | 3 | 3 | Mixed | |
3.41421 | 5 | Mixed | ||
3.42705 | 4 | Tangent |
iii.3.2 Higher Order Icosahedra
As the icosahedron is subdivided and the total number of vertices grows, the classification of modes into irreducible representations remains similar, but each irreducible representation now occurs many times. Fig 4a shows the lowest frequency modes for a icosahedron with vertices. To obtain this figure, we set , and varied . For each value of we relaxed the structure to mechanical equilibrium using steepest descents, evaluated the Hessian matrix by numerical differentiation, then diagonalized the matrix. The Foppl-von Karman number is defined as in eq. (39), where now is defined as the root-mean-square radius (defining instead as the mean radius Lidmar et al. (2003) has little impact below or near the buckling transition and results only in a slight rescaling as grows large) and takes values in the range 6.6-7.6 for the icosahedron.
Owing to rotation and translation invariance of the total energy, we always have a 6-fold degenerate mode of zero eigenvalue. The remaining eigenvalues fall into the classification of icosahedral symmetry introduced in Table 2.
At low , when the shape is spherical in the continuum limit of large radius, and the energy cost of bending dominates over the energy cost of stretching or shearing, the lowest frequency nondegenerate mode is a “breathing” mode, corresponding to a sphere with oscillating radius. Perturbing the radius by an amount (i.e. adding mode ) increases the energy by while displacing vertices by . Identifying the energy with , and noting the area modulus , we find eigenvalue which fits well to the data in Fig. 4.
At higher frequencies, where the wavelength of the modes becomes small compared to the radius of curvature, we expect that the eigenvalues should revert to their flat space limits as discussed in section II.2. The validity of this hypothesis is demonstrated in the dispersion relations shown in fig. 5. Here we plot the vibrational frequencies (i.e. the square roots of Hessian eigenvalues) as functions of the equivalent wave number, defined as the angular momentum index divided by the radius . The radii of the circles represent the projections of the eigenvectors onto the vector spherical harmonics (top), and the longitudinal and transverse eigenfunctions and (middle and bottom). The solid lines are the predictions of continuum elastic theory for the plane, eqs (20-22).
Soft-mode behavior at the buckling transition is less pronounced than in the case of the cone. The crossover from spherical to faceted shape, which occurs gradually for , preserves the icosahedral symmetry. As such, the displacements respect icosahedral symmetry. If the transition is due to a “soft mode”, this mode itself must be invariant under operations of the icosahedral symmetry group. That is, it must transform as the unit representation and therefore must be nondegenerate. The soft mode is best seen in figure 4b, where only the nondegenerate modes are shown. Always the lowest frequency nondegenerate mode is an “breathing” mode, and as just discussed its frequency does not depend significantly on . However, the next occurrence of the unit representation, at , contains a mode, of type and labeled , that does indeed soften and mixes with the breathing mode in an avoided crossing around . Another mode, of type and labeled , is prevented by symmetry from mixing with the mode. A series of other nondegenerate modes () soften at higher values and mix with the other soft modes.
Around the buckling mode consists predominantly of and spherical harmonics, with a small admixture of and higher harmonics. The weight of this mode is concentrated in the vicinities of the icosahedron vertices, and it has strong overlap with the displacements of vertices under the buckling transition.
The forbidden crossing of the buckling and breathing mode smears of the buckling transition, because prevents the eigenvalue of the buckling mode from actually crossing zero. This contrasts with the case of the disclinated flat sheet buckling into a cone, where the eigenvalue does indeed cross zero. For the sheet-to-cone transition the analogue of the breathing mode is just a zero energy translation, rather than a finte frequency radial displacement. Also, up-down symmetry of the plane allows the crossing of the buckling mode (which is odd) with this translation. On a sphere the symmetry breaking between inside and outside the sphere causes the breathing mode to mix with the buckling.
Iv Susceptibilities
iv.1 Cones
The soft-mode transition is a genuine sharp phase transition for the buckling of a disclinated sheet into a cone. We already discussed the order parameter (height) and energy variation through the transition, in section III.2. Now we consider the susceptibility, namely the response of the order parameter to an applied field. In this case we examine the response of the buckling height to a point force applied at the disclination.
Assume the height of the cone (i.e. the vertical displacement of the 5-coordinated particle at the center) is given by , where again is the amplitude and measures the projection of the mode onto the height variable. Then in the presence of an applied force we express the energy as
(48) |
which differs from (45) by the work done against the applied force. Minimizing the energy yields resulting in total height and susceptibility
(49) |
Thus a vanishing eigenvalue, say , passing linearly through zero at , translates immediately into a diverging susceptibility. This divergence is evident in Fig. 3b. Note that the amplitudes differ on the two sides of because in one case we perturb around a flat network while in the other we perturb around a buckled cone.
iv.2 Spheres
Now consider the analogous response for the case of icosahedrally symmetric triangulated spheres and faceted icosahedra. We consider the inverse of the susceptibility as an effective spring constant , and make the spring constant dimensionless by dividing by the Young’s modulus . We first present numerical results for symmetric forces over a wide range of values, then in later subsections consider the limiting cases of small and large .
The data was generated from a sequence of icosahedra of varying sizes and elastic properties. We consider nonchiral icosahedra with ranging from 4 to 512 (i.e. ranging from 482 to 7864322). To speed calculation, the full 120-element icosahedral symmetry group was employed, resulting in a speedup of nearly 120 times. Beyond the buckling transition the five-coordinated vertices sharpen, with a radius of curvature related to the buckling radius
(50) |
by a geometrical factor of order 1. In order to approximate the continuum limit we chose to hold fixed and keep , resulting in .
Each structure was relaxed using a conjugate-gradient method. We found that the necessary number of relaxation steps diverges with increasing size, consistent with the vanishing relaxation rate predicted by eq. (23). To ensure sufficient accuracy in the susceptibility, we used 128-bit real arithmetic in the final stages of all relaxations. For , complete relaxation requires approximately two months on a 3.0GHz Intel Xeon computer. For studies such as ours which seek the continuum limit, a finite element aproach DiDonna and Witten (2001); DiDonna (2002); Klug et al. (2006) might be more computationally efficient than our discrete mass and spring model.
Once the structure was relaxed without applied stress, we re-relaxed with a radially inward force applied symmetrically at all vertices, all faces or all edges. The effective spring constant was defined as the applied force divided by the displacement of the mass to which the force was applied. We actually consider because we define as the derivative of with respect to all simultaneous applied forces . Small applied force was required to achieve linear response in cases where became small.
Figure 6 shows numerical data for symmetric forces applied to vertices, edges or faces. In the limit of small the three data sets converge to a -independent value. As increases, the vertices weaken more quickly than the faces or edges, consistent with our picture of the buckling transition as concentrating at the disclinations which are located at vertices. However, beyond the vertex stiffness falls off very slowly, while both face and edge stiffness continue their rapid decline.
iv.2.1 Small limit
The following discussion first considers the limit of small , in which the shapes are nearly spherical and calculations can be done exactly. The response depends on whether the stress is applied in a uniaxial manner (e.g. at diametrically opposed points) or in a more symmetric manner (e.g. applied simultaneously at all vertices or faces or edges, or even an isotropic pressure).
For an applied pressure the deformation is purely radial, with amplitude as in the breathing mode discussed previously. This increases the energy by , while doing work against the pressure. Balancing the two yields , susceptibility and spring constant . In the case of symmetrically applied point forces, we identify yielding , where we used eqs. (38) and (40). The stiffness is independent of , consistent with the numerical result shown.
For uniaxial stress, let the displacement at the two poles be and assume this displacement persists over a polar region of size (see section 15 of Ref. Landau and Lifshitz (1986)). The bending energy density , and integrating over the polar region yields total bending energy . Meanwhile the strain tensor yields a total stretching energy (see eq. (28)) . Minimizing the sum to find the optimal shape yields and . Equating this to , the work done against the applied force, we find and . The elastic constant , independent of the axis along which the force is applied, consistent with our numerical results (see Fig. 6, inset).
iv.2.2 Large limit
For the radius of curvature at the icosahedron vertices quickly approaches (eq. 50) and remains fixed independent of the icosahedron radius . Forces applied at icosahedron vertices get transfered through the curved vertex region to the flat facets in a primarily longitudinal manner. According to the theory of longitudinal deformation of plates (see section 13 of Ref. Landau and Lifshitz (1986)) the displacement at large distances from the applied force varies as with some fixed length. Upon setting and choosing proportional to , we find that . The numerical data shown in Fig. 6 fits well to this form with values =17.3 and =0.79. The curve shown for comparison illustrates imposing a uniform displacement for visual clarity.
For forces applied to the icosahedron edges we expect to see the onset of ridge scaling behavior Witten and Li (1993); Lobkovsky et al. (1995); DiDonna and Witten (2001); DiDonna (2002); Wood (2002); Witten (2007) as approaches . Unfortunately the diverging relaxation time prevents us from exploring larger within our current calculational method, preventing us from observing this behavior cleanly. We briefly review the predictions of ridge scaling.
Let (which is proportional to ) be the length of an icosahedron edge. At each end of this edge the facets join at a fixed angle of . At the middle the edge sags inward by an amount , creating a saddle shaped ridge with a small radius of curvature across the ridge and a large (and negative) radius of curvature along the ridge Witten and Li (1993). The strain along the ridgeline is of order . Because the facets on either side of the ridge approach angle , the radius is proportional to the sag Witten and Li (1993). Assuming that the bending and strain energy persist along the length of the ridge and extend a distance to either side, we estimate the energy as
(51) |
where the final term represents the action of a force acting at mid-edge.
Upon setting and varying to minimize the energy, we find, in the absence of force ,
(52) |
In the presence of weak applied force , the small radius increases by an amount of order
(53) |
Recalling that and converting this to an effective spring constant yields . Indeed, the edge elasticity in Fig. 6 seems to show a crossover towards slope on our log-log plot.
Meanwhile, the icosahedron faces become almost planar in the limit of large . Timoshenko Timoshenko (1940) discusses the deflection of an equilateral triangular plate under load applied at the center. The deflection is proportional to , from which we conclude, using eq. (1), that . However, in Fig. 6 the face elasticity seems to follow a power law closer to -0.8 than -1. Perhaps residual stresses in the faces or on their boundaries are responsible for this difference.
V Conclusions
In summary, we investigated the eigenvalue spectrum of a simple mass-and-spring model of a virus capsid as it passes through its buckling transition. The buckling of a spherical shell occurs in a smooth, nonsingular fashion, in contrast to the buckling of a disclinated planar network. The smearing can be attributed to symmetry-breaking between the interior and exterior of the shell and is caused by the forbidden crossing of the buckling mode with a lower frequency breathing mode.
Symmetries of the icosahedron and analogies with continuum elastic theory were used to classify the normal modes. Modes of full icosahedral symmetry, transforming as the unit representation, soften as the Foppl-von Karman number passes through the buckling transition. Displacements during buckling, which resemble the maturation of real virus capsids, can be well represented as a superposition of the two lowest icosahedrally symmetric modes.
Susceptibilities to applied forces diverge at the buckling transition for planar networks. For spherical topology they evolve smoothly, with anomalies in the vicinity of . Susceptibility to forces applied at icosahedron vertices dominates near , but icosahedron edges and faces are much softer for large . In the limit of small the effective spring constant approaches the behavior of a spherical continuum.
Beyond the buckling transition the faces have the softest linear response, so this is where one might expect rupture in response to an isotropic osmotic pressure. Relative softness of icosahedron faces as compared to vertices has been reported experimentally in liposomes Delorme et al. (2006). We verified this numerically by calculating the parameter which measures the distortion from a sphere to an icosahedron Lidmar et al. (2003). Below isotropic pressure weakly increases the value of , while above pressure strongly decreases , bending the facets to make the shape more nearly spherical.
Acknowledgements.
Work by MW was supported in part by NSF grant DMR-0111198. Work by DRN was supported by NSF through grant DMR-0231631 and through the Harvard Materials Research Science and Engineering Center via grant DMR-0213805.References
- Caspar and Klug (1962) D. L. Caspar and A. Klug, Cold Spring Harbor Symposia Quant. Bio. 27, 1 (1962).
- Dinsmore et al. (2002) A. D. Dinsmore, M. F. Hsu, M. G. Nikolaides, M. Marques, A. R. Bausch, and D. A. Weitz, Science 298, 1006 (2002).
- Spector et al. (1996) M. S. Spector, J. A. Zasadzinski, and M. B. Sankaram, Langmuir 12, 4704 (1996).
- Delorme et al. (2006) N. Delorme, M. Dubois, S. Garnier, A. Laschewsky, R. Weinkamer, T. Zemb, and A. Fery, J. Phys. Chem. B 110, 1752 (2006).
- Lidmar et al. (2003) J. Lidmar, L. Mirny, and D. R. Nelson, Phys. Rev. E 68, 051910 (2003).
- Nguyen et al. (2005) T. T. Nguyen, R. F. Bruinsma, and W. M. Gelbart, Phys. Rev. E 72, 051923 (2005).
- Vliegenthart and Gompper (2006) G. A. Vliegenthart and G. Gompper, Biophysical J. 91, 834 (2006).
- Hicks and Henley (2006) S. D. Hicks and C. L. Henley, Phys. Rev. E 74, 031912 (2006).
- Timoshenko (1940) S. Timoshenko, Theory of plates and shells (McGraw-Hill, 1940).
- Love (1944) A. E. H. Love, A Treatise on the mathematical Theory of Elasticity (Dover, 1944).
- Niordson (1985) F. I. Niordson, Shell Theory (North-Holland, 1985).
- Seung and Nelson (1988) S. Seung and D. R. Nelson, Phys. Rev. A 38, 1005 (1988).
- Wikoff et al. (2006) W. R. Wikoff, J. F. Conway, J. Tang, K. K. Lee, L. Gan, N. Cheng, R. L. Duda, R. W. Hendrix, A. C. Steven, and J. E. Johnson, J. Struct. Biol. 153, 300 (2006).
- Ivanovska et al. (2004) I. L. Ivanovska, P. J. de Pablo, B. Ibarra, G. Sgalari, F. C. MacKintosh, J. L. Carrascosa, C. F. Schmidt, and G. J. L. Wuite, Proc. Nat. Acad. Sci. 101, 7600 (2004).
- Klug et al. (2006) W. S. Klug, R. F. Bruinsma, J.-P. Michel, C. M. Knobler, I. L. Ivanovska, C. F. Schmidt, and G. J. L. Wuite, Phys. Rev. Lett. 97, 228101 (2006).
- Tama and Brooks (2002) F. Tama and C. L. Brooks, J. Mol. Biol. 318, 733 (2002).
- Witten and Li (1993) T. A. Witten and H. Li, Europhys. Lett. 23, 51 (1993).
- Lobkovsky et al. (1995) A. Lobkovsky, S. Gentges, H. Li, D. Morse, and T. A. Witten, Science 270, 1482 (1995).
- DiDonna and Witten (2001) B. A. DiDonna and T. A. Witten, Phys. Rev. Lett. 87, 206105 (2001).
- DiDonna (2002) B. A. DiDonna, Phys. Rev. E 66, 016601 (2002).
- Wood (2002) A. J. Wood, Physica A (2002).
- Witten (2007) T. A. Witten, Rev. Mod. Phys. 79, 643 (2007).
- Widom (1986) M. Widom, Phys. Rev. B 34, 756 (1986).
- Rader et al. (2005) A. J. Rader, D. Vlad, and I. Bahar, Structure 13, 413 (2005).
- Tama and Brooks (2005) F. Tama and C. L. Brooks, J. Mol. Biol. 345, 299 (2005).
- Peterson (1984) M. A. Peterson, J. Math. Phys (1984).
- Landau and Lifshitz (1986) L. D. Landau and E. M. Lifshitz, Theory of Elasticity (Pergamon, 1986).
- Levine and MacIntosh (2002) A. J. Levine and F. C. MacIntosh, Phys. Rev. E 66, 061606 (2002), see also E. Frey and D. R. Nelson, J. Phys I France, 1715 (1991).
- Hill (1954) E. L. Hill, Am. J. Phys. 22, 211 (1954).
- Abromowitz and Stegun (1970) M. Abromowitz and I. A. Stegun, Handbook of Mathematical Functions (Dover, 1970).
- Tinkham (1964) M. Tinkham, Group theory and quantum mechanics (McGraw-Hill, 1964).