Jun 28, 2007 - Strains associated with the disclinations cause the network to buckle [12], transforming the .... Given the elastic free energy F we ob...

0 downloads 4 Views 1MB Size

arXiv:0706.4291v1 [physics.bio-ph] 28 Jun 2007

Carnegie Mellon University Pittsburgh, PA 15213 USA Department of Computational Biology, School of Medicine University of Pittsburgh, Pittsburgh, PA 15213 USA

J. Lidmar Department of Physics, Royal Institute of Technology AlbaNova, SE-106 91 Stockholm Sweden

David R. Nelson Department of Physics, Harvard University Cambridge, MA 02139 USA (Dated: May 25, 2018) 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

2 icosahedral symmetry. PACS numbers:

I.

INTRODUCTION

Virus capsids [1] and other structures such as colloidosomes [2] and liposomes [3, 4] consist of thin shells of spherical topology that frequently exhibit icosahedral symmetry. A popular simplified model [5, 6, 7, 8] 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. Fivecoordinated vertices may be considered as +2π/6 disclinations within an otherwise sixcoordinated lattice. These disclinations are absent in conventional continuum models of spherical shells [9, 10, 11]. Elastic properties of the capsid can be mimicked by suitably adjusting the spring constants to obtain the desired Young’s modulus Y and by imposing a curvature energy to obtain the bending modulus κ. Strains associated with the disclinations cause the network to buckle [12], transforming the shape from smooth and nearly spherical to strongly faceted

FIG. 1: Triangulated network of P = Q = 2. Colors identify local environments with 5-fold vertices shown in blue.

3

FIG. 2: Shell shape above the buckling transition for P = 128, Q = 0 shell with ks = 1 and kb = 16 yielding Foppl-von Karman number γ=930. Color coding is logarithmic according to total elastic energy (violet=low, red=high) .

and nearly icosahedral [5]. A dimensionless parameter controls the transformation. We define the Foppl-von Karman number γ=

Y R2 κ

(1)

where R is a linear dimension of the shell. The buckling occurs when γ exceeds a value γb of order 102 (see Fig. 2). For the virus HK97, which appears to facet as it matures [13, 14], γ reaches a value of order 103 according to the estimate of Ref. [5]. Varying the pH of solution can alter γ, with the range 100-900 reported for the virus CCMV [15, 16]. At much larger values of γ (in excess of 106 ) which should characterize liposomes with crystalline order, an interesting phenomenon known as “ridge scaling” emerges [17, 18, 19, 20, 21, 22]. Caspar and Klug [1] classify icosahedral structures by a pair of integers (P, Q). A pair of five-coordinated vertices is connected by a path consisting of P edges in some given direction and Q edges in a direction 60◦ to the left (e.g. between two blue vertices via a red vertex in Fig. 1). The T -number of the network, T = P 2 + P Q + Q2 gives the number of vertices as Nv = 10T + 2. There are always 12 five-coordinated vertices, so the number of six-coordinated vertices is 10(T − 1). Structures with P and Q both nonzero and P 6= Q are chiral, such that (P, Q) and (Q, P ) are mirror images. Their symmetry group is the 60-element icosahedral rotation group Y . Structures with either P or Q = 0, or with P = Q

4 are nonchiral. Their symmetry belongs to the 120-element group Yh = Y × Z2 , which should not be confused with the 120-element icosahedral double group [23] Y ′ .

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 Y , 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 γb. Other studies consider more microscopis elastic network models [24, 25] that place nodes at every C α 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 disclinationfree 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. [5] 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 γb for forces applied at five-fold symmetry axes, but the stiffness falls off more rapidly for forces applied at twoand three-fold symmetry axes for γ > γb . We analyze these susceptibilities in limiting cases of small and large γ.

5 II.

CONTINUUM-ELASTIC THEORY

The general elasticity theory of membranes can be expressed in coordinate-free form [11, 26]. Let M be a manifold (a two-dimensional smooth surface embedded in three dimensional space) assumed to be in mechanical equilibrium. Now impose tangential deformation u(x) and normal deformation ζ(x) corresponding to displacements of points x on the surface. Let gαβ and Cαβ be the metric and curvature tensors respectively of M after distortion. Greek indices take values 1 and 2 corresponding to the dimensions of M. Define the strain tensor Uαβ = uαβ + ζCαβ

(2)

where uαβ = 21 (Dα uβ + Dβ uα ) and Dα indicates covariant differentiation with respect to xα . The trace Uγγ measures dilation, while 1 Sαβ = Uαβ − gαβ Uγγ 2

(3)

measures shear strain. Bending of M is characterized by mean curvature H = 21 Cγγ and Gaussian curvature K = det C. The free energy density at x contains dilation, shear and bending contributions, f (x) = fd + fs + fb and can be integrated over M to obtain the total free energy Z p F = f (x) det gd2 x.

(4)

(5)

The separate contributions are 1 fd (x) = (λ + µ)(Uγγ )2 2 fs (x) = µS αβ Sαβ 1 fb (x) = κ(2H − c0 )2 + κG K. 2

(6)

The elastic constants λ and µ are the Lame constants [27]. The 2D area (bulk) modulus B = λ + µ, while µ itself is the shear modulus, and the 2D Young’s modulus Y = 4µ(λ + µ)/(λ + 2µ). Upon integration over the surface M, the Gaussian curvature term becomes constant, and we neglect this term henceforth. Likewise, we set the spontaneous curvature

6 c0 = 0, thus assuming the manifold M would be flat in the absence of constraints associated with the spherical topology. Effects of c0 6= 0 are discussed in Ref. [6]. Given the elastic free energy F we obtain the stress tensor σ αβ (x) =

δF δUαβ (x)

(7)

whose divergence yields the tangential force F β = Dα σ αβ .

(8)

The normal force is given by N(x) = −

δF . δζ(x)

(9)

In mechanical equilibrium the stress tensor and normal force vanish. Slightly out of equilibrium, to first order in the displacements, special forms of u(x) and ζ(x) known as normal modes solve the eigenvalue equation (F, N) = −Λ(u, ζ).

(10)

When displaced from equilibrium according to the k th normal mode (uk , ζk ), the free energy increases by 1 ∆F = Λk 2

Z

(|uk |2 + ζj2)dx.

(11)

According to the equipartition theorem the modes fluctuate with thermal energy ∆F = R kB T /2 and amplitude (|uk |2 + ζk2 )dx = 2kB T /Λk . Time dependence of the strains depends on the equations of motion. In the overdamped

case we write ˙ u˙ β (x) = ΓF β (x), ζ(x) = ΓN(x)

(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 [28] carries out a more thorough investigation of flat membranes coupled to fluid flow. In the absence of damping we write ¨ ρ¨ uβ (x) = F β (x), ρζ(x) = N(x) with ρ the 2D mass density. A normal mode now oscillates in time at frequency ω =

(13) p

Λ/ρ.

7 A.

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 u and perpendicular displacement ζ. 1 1 f = λ(uγγ )2 + µuαβ uαβ + κ(∆ζ)2 2 2

(14)

Here ∆ = Dα D α = ∇2 is the usual 2D laplacian operator and ∇ the usual gradient. By differentiating the energy we obtain the forces F = (λ + µ)∇∇ · u + µ∆u

(15)

N = −κ∆2 ζ

(16)

and

Because the in-plane and out-of-plane displacements decouple, we solve them separately. The solutions are based on the plane wave function ψk (r) = eik·r

(17)

which is an eigenfunction of the Laplacian operator ∆ψk (r) = −k 2 ψk (r). In-plane normal modes are expressed as longitudinal waves uL (r) = ∇ψk (r) = ikeik·r

(18)

uT (r) = zˆ × uL = i(kx yˆ − ky xˆ)eik·r .

(19)

and transverse waves

Note the identities ∇ × uL = 0 and ∇ · uT = 0, 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 ΛL = (λ + 2µ)k 2

(20)

ΛT = µk 2 .

(21)

uP (r) = zˆψk (r)

(22)

and

Perpendicular out-of-plane waves

8 obey eq. (16) subject to the perpendicular wave dispersion relation ΛP = κk 4 .

(23)

For future reference we recast the normal modes in plane-polar coordinates (r, φ), replacing the plane-wave function ψk (r) with cylindrical Bessel functions ψkm (r, φ) = Jm (kr)eimφ .

(24)

The Laplacian operator takes the form 1 ∂ ∆= r ∂r

∂ r ∂r

+

1 ∂2 . r 2 ∂φ2

(25)

Like the plane wave function ψk (r), waves of type (24) are eigenfunctions of the Laplacian operator, ∆ψkm (r, φ) = −k 2 ψkm (r, φ). Upon defining normal modes uL , uT , uP 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.

B.

Deformations of a Sphere

Now we redo the prior calculation of section II A for the case of small perturbations around a sphere of equilibrium radius R. In this case the unperturbed manifold has constant mean curvature H0 = 1/R. Consequently the free energy acquires a term coupling the in-plane and normal strains through the dilation energy. 1 fd = (λ + µ)(uγγ + 2ζ/R)2 2 fs = µ(uαβ uαβ − (uγγ )2 ) p 2 2 1 2 2 fb det g = κ ( − ∆ζ) − 2 (Dα ζ) 2 R R

(26)

In the above, the Laplacian operator takes the form 1 ∂2 ∂ 1 ∂ ∆= 2 sin θ + 2 2 . (27) R sin θ ∂θ ∂θ R sin θ ∂φ2 √ Notice we include the integration measure det g along with the bending energy fb , because √ it contributes the term (Dα ζ)2. The det g factor is not needed in fd or fs because these are already second order in the deformation.

9 Taking functional derivatives of F yields the stress tensor, in-plane and normal force 1 2 ζ) + 2µ(uαβ + g αβ ζ) R R 1 2 F β = (λ + µ)D β (uγγ + ζ) + µ(∆ + 2 )uβ R R 2 γ 4 N = −(λ + µ)( uγ + 2 ζ) − κLζ R R σ αβ = λg αβ (uγγ +

(28)

where we define L = Dα D α D β D β +

2 Dα D α . 2 R

(29)

The extra µuβ /R2 in eq. (28) for F β 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 fb det g. Take the spherical harmonic Ylm (θ, φ) as the basic deformation, analagous to the plane wave eik·r in eq. (17) or the cylindrical wave Jm (kr)eimφ in eq. (24). The spherical harmonic is an eigenfunction of ∆ with eigenvalue −l(l+1)/R2 and an eigenfunction of L with eigenvalue

l(l − 1)(l + 1)(l + 2)/R4 . By analogy with the procedure for plane waves in flat space, we take derivatives as uL = R∇Ylm uT = ˆr × uL

uαL = RD α Ylm

(30)

uαT = Rǫαβ D β Ylm

where ǫ is the alternating tensor. We also define uP = ˆ rYlm

(31)

These functions are linear combinations of the “Vector Spherical Harmonics” Vlm , Wlm and Xlm , which form a complete set of orthogonal functions for expanding vector fields on the surface of a sphere [23, 29]. Notice that the transverse mode uT is proportional to the angular momentum operator acting on Ylm , thus identifying it with the vector spherical harmonic Xlm . The longitudinal and perpendicular modes, uL and uP , are linear combinations of Vlm and Wlm . Note that the longitudinal and transverse modes uL and uT exist only for l ≥ 1, while uP exists for l ≥ 0.

The transverse mode uT is divergenceless (uγγ = 0) and hence creates no perpendicular

force N and no longitudinal force (the gradient part of F β ). In fact, it is an eigenfunction

10 of the force (28). Upon taking into account the commutation of covariant derivatives, we find F β = [µ(l − 1)(l + 2)/R2 ]uβT from which we obtain the eigenvalue ΛT = µ

(l − 1)(l + 2) R2

(32)

As expected, ΛT = 0 for l = 1 because these modes correspond to rigid rotations. In contrast to the transverse modes, the longitudinal and perpendicular modes uL and uP are coupled in both the tangential force F β and perpendicular force N. In matrix form, α α u M MLP F L . = LL (33) ζ MP L MP P N Setting uL as in eq. (30) and setting ζ as the radial component of uP in eq. (31), the matrix elements of M become MLL = (λ + µ) l(l+1) + µ (l−1)(l+2) R2 R2 MLP = (λ + µ) R22 MP L = (λ + µ) 2l(l+1) R2

(34)

MP P = (λ + µ) R42 + κ (l−1)l(l+1)(l+2) R4 The eigenvalues of this matrix, λ± , are the desired normal mode eigenvalues Λ. For the special case l = 1 the eigenvalues are λ− = 0 and λ+ = 6(λ+µ)/R2 . 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 R → ∞. This correspondence can be verified by holding r = Rθ, k = l/R and m fixed, and noting [30] r 4π Ylm (θ, φ) = (−1)m Jm (kr)eimφ . lim l→∞ 2l + 1

(35)

In addition, the eigenvalues should approach their proper limits. Clearly ΛT approaches its flat space limit (21). To check ΛL,P , note that the eigenvalues λ± of the matrix (34) approach (λ + 2µ)l(l + 1)/R2 and κ(l − 1)l(l + 1)(l + 2)/R4 in the limit of large sphere radius R, yielding the flat space limits Eqs. (20) and (23).

11 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 ri be the position of mass i and n ˆ I be the orientation of plaquette I. A plaquette is a set of three masses joined to each other by springs, and we take the normal in the outward direction. Following [5] we define Hs =

ks X (|ri − rj | − a)2 2

(36)

kb X |ˆ nI − n ˆ J |2 2

(37)

hiji

and Hb =

hIJi

and set the unstretched spring length a = 1. Here hiji denote pairs of nearest-neighbor vertices, and hIJi denote pairs of adjacent (edge-sharing) plaquettes. In the continuum limit the discrete model reproduces the continuum system with elastic constants √ 3 2 kb , κ= Y = √ ks 2 3 Foppl-von Karman number γ=

4ks R2 Y R2 , = κ 3kb

(38)

(39)

Lame coefficients and bulk modulus √

3 λ=µ= ks 4

B=

√

3 ks , 2

(40)

and 2D mass density (taking the vertex mass m = 1) √ ρ = 2/ 3.

A.

(41)

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

12 sound wave. The dispersion relation is simplest for wavevector k in the yˆ direction (chosen to lie midway between two near-neighbor bonds) ΛL = 3ks (1 − cos (

√

3 ky a)) 2

(42)

Taking the cross product with zˆ yields the transverse sound wave with dispersion relation √ 3 ky a)) (43) ΛT = ks (1 − cos ( 2 Finally, taking the perpendicular displacements as the planewave yields the perpendicular modes with dispersion relation ΛP = kb (2 − 2 cos (

√

3 ky a))2 2

(44)

In the continuum limit ka << 1 these dispersion relations revert to the prior results of continuum elastic theory.

B.

Buckling of a Plane into a Cone

Upon introducing a five-fold +2π/6 disclination into the flat triangulated network discussed previously, strain energy accumulates [12] and grows without bound as the radius R of the network increases. At some specific “buckling radius” Rb 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 R fixed and vary the bending stiffness. Large kb opposes buckling and the network lies flat, while small kb allows buckling out of plane into a cone. For the network of radius R = 8a analyzed below, buckling occurs for kb ≈ 0.71. As R increases the threshold value of kb grows as R2 so that γ approaches the limiting value γb ≈ 154 [5, 12]. Eigenvectors of the Hessian matrix form basis functions for representations of the symmetry group of the structure [31]. Eigenvectors sharing a common eigenvalue form the basis for an irreducible representation. Thus the patterns of degeneracy follow the dimensionalities of

13

Eigenvalue

0.075

2x

0.05

1x 6x 0.025

0

0

0.2

0.4

0.6

0.8

1

5

150 Susceptibility 100

3 Energy Height E=0.7361

2

50

Susceptibility

Energy, Height

4

1 0

0

0.2

0.4

0.6

0.8

1

0

kb

FIG. 3: Triangulated network of radius 8a and spring constant ks = 1 with a single 5-fold disclination at center. (top) Eigenvalue spectrum color coded according to degeneracy. Note the nondegenerate 1x mode that goes to zero at the buckling transition. (bottom) energy, cone height and susceptibility.

14 C5v m 1C0 2C5 2C52 5σv A1 0

1

1

1

1

A2 0

1

1

1

-1

E1 1

2 τ −1 −τ

E2 2

2

0

−τ τ −1 0

TABLE I: Character table of C5v . Cn denotes conjugacy class of order n. Values of m denote √ in-plane angular momenta. τ = ( 5 + 1)/2 is the Golden Mean.

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 C5v 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 I). For the specific case of the flat network, the group is even higher, D5h , 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. For the group D5h , the irreducible representations are based on those of C5v supplemented with an additional label g, u according to whether they are even (g) or odd (u) under reflection through the horizontal plane σh . The requirement that each irreducible representation be either even or odd under σh requires that each mode be polarized either fully in-plane or fully perpendicular. Let Λ1 be the lowest nonzero eigenvalue. Its eigenvector e1 is polarized strictly perpendicular to the sheet and transforms as the irreducible representation A2u . Its value is nonzero at the origin. The energy of mode i varies as Λi a2i where ai 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 Λ1 > 0. However, if we monitor the value of Λ1 as a function of γ (Fig. 3a) we find it crosses through zero at γb .

15 For small deformations we express the energy as E=

X1 i

2

Λi a2i + O(a4i )

(45)

Now set Λ1 = c(γb − γ). The mechanically stable minimum energy structure is perfectly flat (ai = 0) for γ < γb , but it buckles out of plane for γ > γb , in a shape described by the √ eigenvector e1 , with amplitude growing as γ − γb. Meanwhile the energy drops as (γ −γb )2 . These effects can be seen in figure 3b. For γ > γb , figure 3a shows the spectrum of vibrations around the mechanically stable, buckled structure. Note that Λ′1 (the lowest nondegenerate eigenvalue) becomes positive again.

C.

Buckling of Spherical Shells 1.

P = 1, Q = 0 icosahedron

Table II presents the character table of the 60-element icosahedral rotational symmetry group Y , which has 5 irreducible representations. The conjugacy classes are labeled Cn , where n is the order of elements in the class, so that an element of Cn corresponds to a rotation by 2π/n. Recall that the spherical harmonics Ylm form basis functions for the irreducible representations of the continuous rotation group SO(3), and therefore they induce representations (in general reducible) of Y . For a given angular momentum l and rotation angle θ, the character is χl (θ) =

sin (l + 1/2)θ . sin θ/2

(46)

Irreducible representations of Y are labeled in Table II according to the lowest angular momentum l under which they transform. Of particular interest is the representation F1 corresponding to angular momentum l = 1. This is the representation under which threedimensional vectors transform. 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 scalarvalued 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 χR of the regular representation equal the number of vertices that remain stationary under

16 Y

l 1C0 15C2 20C3 12C5 12C52

A

0

1

1

1

1

1

F1 1

3

-1

0

τ

−τ −1

F2 (3) 3

-1

0

−τ −1

τ

G (3) 4

0

1

-1

-1

5

1

-1

0

0

R

12

0

0

2

2

V

36

0

0

2τ

−2τ −1

H

2

TABLE II: Character table of Y . Cn denotes conjugacy class of order n. Values of l denote angular momenta. R is the “regular representation” and V the “total vibrational representation” discussed in subsection III C 1

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 R times the representation F1 corresponding to a three dimensional vector. We call the resulting product representation the “total vibrational representation” [23], and its characters χV = χR χF1 . Reducible representations can be decomposed into their irreducible components using orthogonality properties of character tables. In particular we obtain the decomposition V = A ⊕ 3F1 ⊕ F2 ⊕ 2G ⊕ 3H.

(47)

Of the three occurrences of the vector representation F1 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 A 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 III.

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

17 Λ

Formula

0.00000 0

Irrep

g Comments

2 × F1 6 Translations + rotations

√

0.58579 2 − 3 Ha √ 0.76393 5/R2 − 2 F2

5 Mixed contains V2m and W2m 3 Radial contains rY3m

1.00000 1

H

5 Tangent X2m

1.80901 1 + τ /2 √ 2.76393 5/R2

Ga

4 Tangent X3m

A

1 Radial rY00 breathing mode

3.00000 3

F1

3 Mixed V1m

Hb

5 Mixed

Gb

4 Tangent

3.41421 2 +

√

2

3.42705 1 + 3τ /2

TABLE III: Vibrational eigenvalues for P = 1, Q = 0 icosahedron with a = ks = 1, kb = 0. Λ is √ eigenvalue and g is degeneracy. R = 1 + τ 2 /2 = 0.95106 is the radius of the icosahedron.

now occurs many times. Fig 4a shows the lowest frequency modes for a P = 8, Q = 0 icosahedron with Nv = 642 vertices. To obtain this figure, we set ks = 1, and varied kb . For each value of kb 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 R is defined as the root-mean-square radius (defining R instead as the mean radius [5] 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 P = 8, Q = 0 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 II. 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 uP = ˆ rζ) increases the energy by 8πBζ 2 while displacing Nv vertices by ζ. Identifying the energy with 12 Nv Λbreathing , √ √ and noting the area modulus B = 3/2, we find eigenvalue Λbreathing = 8π 3ks /Nv which fits well to the data in Fig. 4.

18

Eigenvalue

0.2

1x

0.1

4x 3x

Λbreathing

5x 6x

0 1

Eigenvalue

0.8 0.6

100

1000

All 1x modes Breathing mode ~ Y00 Buckling mode ~ Y6m Tangent mode ~ X6m Mixing mode ~ Y10m Mixing mode ~ Y12m

0.4 0.2

Λbreathing

0

100

γ

1000

FIG. 4: Lowest frequency modes of P = 8, Q = 0 icosahedron with Nv = 642 vertices. (top) Color coded according to degeneracy. (bottom) Nondegenerate modes only. Arrows locate eigenvalue √ Λbreathing = 8π 3/Nv . Note the buckling mode (red) dips close to zero near the buckling transition.

19 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 B. 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 l divided by the radius R. The radii of the circles represent the projections of the eigenvectors onto the vector spherical harmonics Xlm (top), and the longitudinal and transverse eigenfunctions uL and uT (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 γ ∼ 100 − 1000, 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 l = m = 0 “breathing” mode, and as just discussed its frequency does not depend significantly on kb . However, the next occurrence of the unit representation, at l = 6, contains a mode, of type uP and labeled Y6m , that does indeed soften and mixes with the breathing mode in an avoided crossing around γ ≈ 400. Another l = 6 mode, of type uT and labeled X2m , is prevented by symmetry from mixing with the uP mode. A series of other nondegenerate modes Ylm (l = 10, 12, 20, ...) soften at higher γ values and mix with the other soft modes. Around γb the buckling mode consists predominantly of l = 0 and l = 6 spherical harmonics, with a small admixture of l = 10 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 Λbuckling > Λbreathing > 0) 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

20

1.5

Frequency

Projection onto XLM sqrt(3ks/8)*q

1

0.5

0 0

0.5

1

1.5

1

1.5

2.5

Frequency

2

Projection onto LongLM sqrt(9ks/8)*q sqrt(3kb/4)*q

2

1.5 1 0.5 0 0

0.5

2.5

Frequency

2

Projection onto TransLM sqrt(9ks/8)*q sqrt(3kb/4)*q

2

1.5 1 0.5 0 0

0.5

1 Wavenumber q=L/R

1.5

FIG. 5: Vibrational frequencies plotted versus wave number q = l/R, where l is the angular momentum index. Data is for P = 8, Q = 0 icosahedron with ks = 1, kb = 0.1, γ = 653. The radii of the circles indicate the sizes of the various projections.

21 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. Owing to the smearing, the value of γb is not uniquely defined for the sphere-toicosahedron transition. Reported values range from 130 to 260 based on fitted energy models [5, 6]. We observe the avoided crossing around γ ≈ 400. IV.

SUSCEPTIBILITIES A.

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 B. 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 P at the center) is given by h = i Pi ai , where again ai is the amplitude and Pi measures the projection of the mode i onto the height variable. Then in the presence of an applied force

we express the energy as E=

X 1 ( Λi a2i − F Pi ai ), 2 i

(48)

which differs from (45) by the work done against the applied force. Minimizing the energy P yields ai = F Pi /Λi resulting in total height h = F i Pi2 /Λi and susceptibility χ=

X P2 ∂h i = . ∂F Λi i

(49)

Thus a vanishing eigenvalue, say Λ1 , passing linearly through zero at γ = γb , translates immediately into a diverging susceptibility. This divergence is evident in Fig. 3b. Note that the amplitudes differ on the two sides of γb because in one case we perturb around a flat network while in the other we perturb around a buckled cone.

22

100 12 π

kb=0.5 kb=1 kb=2 kb=4 kb=8 kb=64 vertex edge face

NK/Y

10

1

K/Y

0.1

1

0.1

0.01

1

10

100

100

1000

1000 10000 2 γ = Y R /κ

1e+05

FIG. 6: Stiffness (inverse susceptibility) to forces applied with icosahedral symmetry at vertices, √ edges and faces. Dashed blue lines show 15/ log (0.79 γ) (vertex, offset from the best fit for √ clarity), 100/ γ (edge, scaling expected for γ > 106 ), and 20000/γ (face, scaling for bending of flat facet [9]). Inset: diametrically opposed forces applied with uniaxial symmetry. Dashed blue √ line shows 5/ γ. B.

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 K = 1/χ, and make the spring constant dimensionless by dividing by the

1e+06

23 Young’s modulus Y . 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 P = Q ranging from 4 to 512 (i.e. Nv ranging from 482 to 7864322). To speed calculation, the full 120-element icosahedral symmetry group Yh was employed, resulting in a speedup of nearly 120 times. Beyond the buckling transition the five-coordinated vertices sharpen, with a radius of curvature Rv related to the buckling radius Rb =

p p γb κ/Y ≈ 12.4 κ/Y

(50)

by a geometrical factor of order 1. In order to approximate the continuum limit we chose to hold ks = 1 fixed and keep kb ≥ 1/2, resulting in Rb ≥ 7.6a. 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 1/R4 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 P = Q = 512, 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 [15, 19, 20] 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 F applied symmetrically at all N = 12 vertices, all N = 20 faces or all N = 30 edges. The effective spring constant K = F/ζ was defined as the applied force F divided by the displacement ζ of the mass to which the force was applied. We actually consider NK/Y because we define K as the derivative of ζ with respect to all N simultaneous applied forces F . Small applied force F = 0.001 was required to achieve linear response in cases where K 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 γb the vertex stiffness falls off very slowly, while both face and edge stiffness continue their rapid decline.

24 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 P the deformation is purely radial, with amplitude ζ as in the breathing mode discussed previously. This increases the energy by 8πBζ 2, while doing work 4πR2 P ζ against the pressure. Balancing the two yields ζ = R2 P/4B, susceptibility χ = ∂ζ/∂P = R2 /4B and spring constant K/Y = 4B/Y R2 . In the case of N symmetrically applied point forces, we identify P = NF/4πR2 yielding NK/Y = 16πB/Y = 12π = 37.7, 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 d (see section 15 of Ref. [27]). The bending energy density fb ∼ κ(ζ/d2)2 , and integrating over the polar region yields total bending energy

Eb = κζ 2/d2 . Meanwhile the strain tensor uαβ ∼ ζ/R yields a total stretching energy (see

eq. (28)) Es ∼ Y (ζ/R)2d2 . Minimizing the sum Es + Eb to find the optimal shape yields √ d4 ∼ (κ/Y )R2 and Es + Eb ∼ κY ζ 2/R. Equating this to F ζ, the work done against the √ √ √ applied force, we find ζ ∼ (R/ κY )F and χ = R/ κY . The elastic constant K/Y ∼ 1/ γ, independent of the axis along which the force is applied, consistent with our numerical results (see Fig. 6, inset).

2.

Large γ limit

For γ > γb the radius of curvature at the icosahedron vertices quickly approaches Rv (eq. 50) and remains fixed independent of the icosahedron radius R. 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. [27]) the displacement at large distances r from the applied force varies as u(r) ∼ (F/Y ) log r/r0 with r0 some fixed length. Upon setting K = F/u(R)

25 √ and choosing r0 proportional to Rv , we find that K/Y ∼ c/ log (b γ). The numerical data shown in Fig. 6 fits well to this form with values c=17.3 and b=0.79. The curve shown for comparison illustrates c = 15 imposing a uniform displacement for visual clarity. For forces applied to the icosahedron edges we expect to see the onset of ridge scaling behavior [17, 18, 19, 20, 21, 22] as γ approaches 106 . 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 L (which is proportional to R) be the length of an icosahedron edge. At each end of this edge the facets join at a fixed angle of θ = 138.2◦. At the middle the edge sags inward by an amount ζ, creating a saddle shaped ridge with a small radius of curvature R1 across the ridge and a large (and negative) radius of curvature R2 along the ridge [17]. The strain along the ridgeline is of order (ζ/L)2. Because the facets on either side of the ridge approach angle θ, the radius R1 is proportional to the sag ζ [17]. Assuming that the bending and strain energy persist along the length L of the ridge and extend a distance R1 to either side, we estimate the energy as E = R1 L[Y (ζ/L)4 + κ(1/R1 )2 ] − F ζ

(51)

where the final term represents the action of a force F acting at mid-edge. Upon setting ζ ∼ R1 and varying R1 to minimize the energy, we find, in the absence of force F , R1 ∼ (κ/Y )1/6 L2/3 ∼

p

κ/Y γ 1/3 ∼ Lγ −1/6 .

(52)

In the presence of weak applied force F , the small radius R1 increases by an amount of order L ∆R1 ∼ √ F (53) κY Recalling that ζ ∼ R1 and converting this to an effective spring constant K = dF/dζ yields p √ K/Y ∼ κ/Y L2 ∼ 1/ γ. Indeed, the edge elasticity in Fig. 6 seems to show a crossover towards slope −1/2 on our log-log plot.

Meanwhile, the icosahedron faces become almost planar in the limit of large γ. Timoshenko [9] discusses the deflection of an equilateral triangular plate under load applied at the center. The deflection is proportional to R2 /κ, from which we conclude, using eq. (1), that K/Y ∼ 1/γ. 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.

26 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 γb . Susceptibility to forces applied at icosahedron vertices dominates near γb , 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 [4]. We verified this numerically by calculating the Q6 parameter which measures the distortion from a sphere to an icosahedron [5]. Below γb isotropic pressure weakly increases the value of Q6 , while above γb pressure strongly decreases Q6 , bending the facets to make the shape more nearly spherical.

Acknowledgments

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

27 Science and Engineering Center via grant DMR-0213805.

[1] D. L. Caspar and A. Klug, Cold Spring Harbor Symposia Quant. Bio. 27, 1 (1962). [2] A. D. Dinsmore, M. F. Hsu, M. G. Nikolaides, M. Marques, A. R. Bausch, and D. A. Weitz, Science 298, 1006 (2002). [3] M. S. Spector, J. A. Zasadzinski, and M. B. Sankaram, Langmuir 12, 4704 (1996). [4] N. Delorme, M. Dubois, S. Garnier, A. Laschewsky, R. Weinkamer, T. Zemb, and A. Fery, J. Phys. Chem. B 110, 1752 (2006). [5] J. Lidmar, L. Mirny, and D. R. Nelson, Phys. Rev. E 68, 051910 (2003). [6] T. T. Nguyen, R. F. Bruinsma, and W. M. Gelbart, Phys. Rev. E 72, 051923 (2005). [7] G. A. Vliegenthart and G. Gompper, Biophysical J. 91, 834 (2006). [8] S. D. Hicks and C. L. Henley, Phys. Rev. E 74, 031912 (2006). [9] S. Timoshenko, Theory of plates and shells (McGraw-Hill, 1940). [10] A. E. H. Love, A Treatise on the mathematical Theory of Elasticity (Dover, 1944). [11] F. I. Niordson, Shell Theory (North-Holland, 1985). [12] S. Seung and D. R. Nelson, Phys. Rev. A 38, 1005 (1988). [13] 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). [14] 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). [15] 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). [16] F. Tama and C. L. Brooks, J. Mol. Biol. 318, 733 (2002). [17] T. A. Witten and H. Li, Europhys. Lett. 23, 51 (1993). [18] A. Lobkovsky, S. Gentges, H. Li, D. Morse, and T. A. Witten, Science 270, 1482 (1995). [19] B. A. DiDonna and T. A. Witten, Phys. Rev. Lett. 87, 206105 (2001). [20] B. A. DiDonna, Phys. Rev. E 66, 016601 (2002). [21] A. J. Wood, Physica A (2002). [22] T. A. Witten, Rev. Mod. Phys. 79, 643 (2007). [23] M. Widom, Phys. Rev. B 34, 756 (1986).

28 [24] A. J. Rader, D. Vlad, and I. Bahar, Structure 13, 413 (2005). [25] F. Tama and C. L. Brooks, J. Mol. Biol. 345, 299 (2005). [26] M. A. Peterson, J. Math. Phys (1984). [27] L. D. Landau and E. M. Lifshitz, Theory of Elasticity (Pergamon, 1986). [28] 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). [29] E. L. Hill, Am. J. Phys. 22, 211 (1954). [30] M. Abromowitz and I. A. Stegun, Handbook of Mathematical Functions (Dover, 1970). [31] M. Tinkham, Group theory and quantum mechanics (McGraw-Hill, 1964).