A general phase plot is proposed for discrete particle shells that allows for thermal fluctuations of the shell geometry and of the inter-particle connectivities. The phase plot contains a first-order melting transition, a buckling transition, and a collapse transition and is used to interpret the thermodynamics of microbiological shells.
The development of shells that protect microbiological systems from a hostile environment yet still allow for exchange of key nutrients was an essential step in the evolution of life . These shells are composed of molecules decorated with both hydrophobic and hydrophilic groups (“amphiphiles”) in such a way that in an aqueous environment they assemble into closed, semipermeable shells. An important example are the amphiphilic protein shells that surround viruses  as well as many bacteria and most archaea . Cryogenic-based microscopy studies  had indicated that these “capsids” in general are strictly organized, crystallographic structures (usually icosahedral or helical) , but this view is being challenged. Thermodynamic studies of viral self-assembly indicate, for the case of the assembly of viruses in solutions with higher protein concentrations, that the interaction energies between capsid proteins (“subunits”) should be no more than a few times the thermal energy at room temperature in order to avoid the production of malformed capsids [5–8]. Finite temperature studies also showed that, due to thermal fluctuations, at least some viral capsids are dynamical in nature and that the dynamics plays a role in the life-cycle of the virus [10,11]. Some capsids are even in a molten or “pleomorphic” state [12,13]. Finally, all-atom molecular dynamics (MD) simulations of capsids revealed that they can collapse under the action of thermal fluctuations .
The study of the melting and thermal collapse of a shell with a limited number of constituent components () is an interesting statistical physics problem in its own right. The geometry of the shell over which the components are moving itself is defined by the position vectors of these same particles and hence subject to thermal fluctuations . Here, we propose a generic phase diagram for the melting and collapse of discrete shells obtained by comparing MD simulations of a coarse-grained model of capsids with the continuum theory of thermally fluctuating surfaces.
We first discuss the MD simulations. The coarse-grained model is based on the so-called oriented particle system (or “OPS”) . An OPS is defined as a cluster of orientable and interacting point particles located at . An orientation-dependent pair interaction acts between particle pairs and with a separation vector and unit vectors and describing their orientations [see Fig. 1 (left)]. The oriented pair interaction used in the simulations was
An OPS can be viewed as a coarse-grained representation of a viral capsid by having the particle locations correspond to the centers of the “capsomers” of viral capsids. The latter are disklike groups of either six or five subunits that frequently act as the basic building blocks of a capsid [4,24]. The orientational degrees of freedom correspond to the normals to the capsomers, the depth of the Morse potential to the capsomer binding energy (a few [25,26]), and the length scale to the diameter of a capsomer (a few nanometers ). Because the range of the hydrophobic attraction between capsomers is short compared to their diameter, the dimensionless parameter characterizing the width of the Morse potential needs to be significantly larger than one. We used . Next, is a measure of the bending stiffness of the shell (estimated to be in the range of [27,28]). Finally, because a large energy penalty is known to obstruct the removal of single capsomers from assembled shells , evaporation of particles from the OPS shell is suppressed by a soft fixed-area constraint imposed via the augmented Lagrange multiplier method (see , Sec. I).
Figure 1 (middle) shows the minimum energy state of an OPS for the case that is large compared to one. The shell has icosahedral symmetry with the twelve blue disks indicating the five-fold symmetry sites of the icosahedron (for actual capsids, these disks would correspond to pentameric protein capsomers while the remaining 20 red disks would correspond to hexameric capsomers). This structure should be compared to that of the “” icosahedral pattern  of the 32 capsomers of the capsid of the Cowpea Chlorotic Mottle virus  shown in the right of Fig. 1.
Next, we carried out Brownian Dynamics simulations of OPS systems using computational methods discussed further in the Supplemental Material  (Sec. I). The phase behavior was determined in terms of the two thermodynamic parameters , a dimensionless measure of temperature in units of the depth of the pair interaction, and a dimensionless measure of the inverse of the bending stiffness ( is the shell radius). In continuum theory is known as the Föppl-von Kármán (FvK) number . For different values of these two parameters we encountered ordered, molten, buckled, and collapsing shells. Representative realizations are shown in Figs. 2 and 3.
The degree of fluidity of a shell was monitored using a dynamical method based on plots of the mean square of the particle separations as a function of time , averaged over all pairs of particles that were nearest neighbors in the initial configuration . For the present case, if in the long time limit saturated (on average) to a constant value much smaller than then the shell was assigned to be in a solid state. If, on the other hand, increases linearly in time until it reaches —which is consistent with particle diffusion—then the shell was assigned to be in a fluid state. Finally, when plots of showed a random sequence of alternating time intervals of saturation and linear growth for a given simulation run with drastic variations between different runs then the shell was assigned to be in a fluid-solid coexistence state. Examples of these three cases are shown in Fig. 2 for . In the low-temperature solid state, with , the shell shape is spherical while in the coexistence regime, with , significant large-scale shape fluctuations are visible. The particle array still maintains local positional order but this has largely disappeared for (fluid state).
The collapse of shells is a pronounced feature of the phase behavior of shells with larger values of the FvK number . Collapse was monitored by computing the volume of a continuous and differentiable surface that interpolates between the particle locations, which was constructed using the Loop shell subdivision method [41,42]. Figure 3(top) shows an example of irreversible collapse induced by thermal crumpling for , as indicated by a drastic reduction in volume over time and the production of very irregular shell shapes. In Fig. 3(bottom) this simulation was repeated at a reduced temperature. After a small initial reduction, the shell volume reached a steady state . The low-temperature shell shape is now icosahedral . In Fig. 3(middle), simulations were performed near the critical temperature for collapse. In this case shell volumes exhibit large fluctuations, with some shells undergoing a first-order-like collapse transition (blue time trace), while other shells remained stable over the simulated time interval.
By collecting simulation runs for different values of the and parameters, the phase plot of Fig. 4 was produced. The vertical bars indicate temperature intervals over which fluid-solid coexistence was observed following the criterion discussed above, with the solid blue dots indicating midpoints. The large coexistence interval for small indicates that the melting transition should be first order on a rigid spherical surface, which is consistent with simulations of flat sheets of particles interacting via Morse potential . The onset of irreversible collapse of shells for increasing temperature is indicated by orange triangles. Shells with collapsed before the particle array could melt. On the other hand, shells with FvK numbers in the range of would melt with increasing temperature before they collapsed.
To interpret the phase plot, we compared it with the thin-shell elasticity theory (TSET)  in which a curved and stretched layer is assigned a bending and stretching energy given by
According to TSET, determines the ground-state shape of a thin elastic shell  such that for less than a critical value (the “buckling threshold”) that is in the range of , the shell has an approximately spherical shape [such as the spherical shape of Fig. 2(bottom)] while for above that threshold, the ground-state shape is approximately polyhedral [such as the icosahedral shape of Fig. 3(bottom)]. As indicated in the phase plot, we find a buckling threshold around . The discrepancy between the computed and predicted values of the buckling threshold, which has been noted before , is a measure of the importance of discreteness effects. To apply TSET to melting, one can combine it with the Lindemann melting criterion (LMC), which states that melting occurs when the rms of the fluctuations of particle displacements exceeds a certain fraction of the equilibrium interparticle spacing . The LMC is known to work well for melting on flat two-dimensional surfaces . In the regime of harmonic fluctuations, the mean square of the in-plane fluctuations and the mean square of the out-of-plane fluctuations can be shown to have the scaling form and , respectively. The scaling functions and are discussed in the Supplemental Material , Sec. II. Formally, TSET theory corresponds to the limit of particle shells with large and small but with fixed radius . To include discreteness effects, we expanded the three displacement fields in terms of a series of spherical harmonics and demanded that the total number of out-of-plane modes —with the maximum value of the quantum number —equals the number of particles minus two . For a shell of particles and are reasonable choices (since 72 is in the interval between and ). The corrected scaling function has the same mathematical form as the TSET case but it is larger by an overall constant scale factor in the range of 10–100, depending on the value of : discreteness strongly amplifies the effects of thermal fluctuations. The resulting LMC melting temperatures are plotted in Fig. 4 (orange line) with the melting temperature for treated as a fitting parameter . The resulting fit is reasonable for the range of values over which melting was observed.
According to TSET, elastic shells with large should undergo a collapse transition with increasing temperature [48–50]. Physically, this is due to the fact that crumpled shells have a much larger configurational space for shape fluctuations than (nearly) spherical shells so their entropy is much larger as well, while the volume of a crumpled shell with fixed area is correspondingly reduced. For larger , the enthalpic cost of crumpling the surface is diminished so a crumpling or collapse transition is to be expected. The collapse of a shell requires overcoming a free energy barrier that according to TSET vanishes when is about (green dashed line in Fig. 4[49,50]). For the crumpling or collapse transition that we observed (red stars), solid shells (but not liquid shells) roughly obeyed this scaling relation except that the value of had to be decreased by a factor of about 10 (purple line in Fig. 4). We interpret this as a discreteness effect similar to the one encountered for the melting transition. In the range of the fitted thin-shell elasticity theory (purple line) provides a reasonable estimate for the collapse transition. Small discrepancies with the molecular dynamics simulation (red stars) are attributed to thermally activated escape events over the energy barrier . This is reflected in the middle panel of Fig. 3, where shell volumes exhibit large fluctuations, with some shells undergoing a first-order-like collapse transition (blue time trace), while other shells remained stable over the simulated time interval. In the range of there is a much larger discrepancy between the thin-shell elasticity and molecular dynamics simulations. For these values of the collapse transition is in the same range as the melting transition. The unbound dislocation pairs that are forming during melting may significantly affect the collapse transition and this effect is not captured in the TSET.
As an example how the phase plot Fig. 4 could be applied to viral capsids, we compared the OPS with what is known about the phase properties of viral capsids having 72 capsomers. It should be noted that our phase plot corresponds to empty viral capsids, while the effect of osmotic pressure due to the packaged DNA would increase both the melting and collapse temperatures due to the suppression of radial fluctuations as discussed in . In the Caspar-Klug classification of capsids, icosahedral shells with 72 capsomers are known as “” structures . Medically important viruses are the human polyoma and papilloma unenvelopeded double-stranded DNA viruses, which both have a diameter of about 50 nm . Because the polyoma and papilloma capsids are quite spherical, the value of should, for these two cases, lie below the buckling threshold in Fig. 4. In contrast, two forms of the capsid of the thermostable DNA bacteriophage P23-45 with diameters of 66 and 82 nm are, respectively, weakly and strongly polyhedral . This progression of shapes straddling the buckling threshold could be understood within TSET by noting that scales as the square of the shell diameter. According to Fig. 4, in the relevant regime of [5,6], capsids should be rather unstable in this part of the phase plot: prone both to melting and collapse. In actuality these viruses are known to be quite stable but this is because of a post-assembly maturation process, during which the capsid subunits are linked together by covalent bonds. The initial procapsids of the viruses, which are not yet bonded together, still could be thermodynamically unstable. Interestingly, the assembly of the procapsids of phages typically takes place on a preformed spherical scaffold that is disassembled during maturation . It is very suggestive that one purpose of the scaffold is to prevent this instability. This could be checked experimentally by a study of the thermodynamic stability of self-assembling mutant empty shells for which maturation and/or scaffold formation is blocked. Instability of the procapsid of the related, but much larger, Herpes-Simplex virus should be even more pronounced.
We would like to thank Luigi Perotti, Jeff Eldredge and Bill Gelbart for useful discussions and the UCLA Mechanical and Aerospace Engineering Department for continued support. This work was supported by NSF through DMR Grants No. 1610384 and No. 1836404, and the CAREER Award No. DMR-1752100.
This is due to kinetic trapping. For lower protein concentrations, the probability of kinetic trapping is reduced and interaction energies can be higher .