RESEARCH ARTICLE Shape selection and mis-assembly in viral capsid formation by elastic frustration Carlos I Mendoza1*, David Reguera2,3 1 Instituto de Investigaciones en Materiales, Universidad Nacional Autónoma de México, México, Mexico; 2Departament de Fı́sica de la Matèria Condensada, Universitat de Barcelona, Barcelona, Spain; 3Universitat de Barcelona Institute of Complex Systems (UBICS), Universitat de Barcelona, Barcelona, Spain Abstract The successful assembly of a closed protein shell (or capsid) is a key step in the replication of viruses and in the production of artificial viral cages for bio/nanotechnological applications. During self-assembly, the favorable binding energy competes with the energetic cost of the growing edge and the elastic stresses generated due to the curvature of the capsid. As a result, incomplete structures such as open caps, cylindrical or ribbon-shaped shells may emerge, preventing the successful replication of viruses. Using elasticity theory and coarse-grained simulations, we analyze the conditions required for these processes to occur and their significance for empty virus self-assembly. We find that the outcome of the assembly can be recast into a universal phase diagram showing that viruses with high mechanical resistance cannot be self- assembled directly as spherical structures. The results of our study justify the need of a maturation step and suggest promising routes to hinder viral infections by inducing mis-assembly. Introduction Viruses are fascinating biological and nanoscale systems (Douglas and Young, 2006; Wen and Steinmetz, 2016). In the simplest cases, these tiny pathogens are formed by a chain of RNA or DNA *For correspondence: cmendoza@iim.unam.mx encased in a protein shell, also known as capsid, made from multiple copies of a single protein (Flint et al., 2004). Despite this apparent simplicity, viruses are able to perform many complex func- Competing interests: The tions which are essential in their replication cycle. One of the most amazing one is their ability to authors declare that no self-assemble with an unparalleled efficiency and precision. competing interests exist. In vivo, the capsid of most viruses assembles from its basic building blocks, which could be indi- Funding: See page 13 vidual capsid proteins, dimers, trimers or capsomers (i.e. clusters of five or six proteins, which consti- Received: 07 October 2019 tute the structural and morphological units of the shell). The resulting structure has a precise Accepted: 06 April 2020 architecture, which in most cases is spherical with icosahedral symmetry (Roos et al., 2010). Several Published: 21 April 2020 viruses assemble their capsid before packaging the genetic material. In addition, the proteins of many viruses have the ability to self-assemble in vitro, even in the absence of genetic material, form- Reviewing editor: Pierre Sens, ing empty capsids. Institut Curie, PSL Research University, CNRS, France The mechanisms of viral assembly have been the subject of recent and interesting investigations (Hagan and Chandler, 2006; Elrad and Hagan, 2008; Nguyen et al., 2009; Johnston et al., 2010; Copyright Mendoza and Perlmutter and Hagan, 2015; Hagan and Zandi, 2016). The assembly of a curved empty shell with Reguera. This article is a well-defined geometry and precise arrangement of the building blocks is a non-trivial process that distributed under the terms of resembles 2D crystallization on a curved space (Meng et al., 2014; Gómez et al., 2015). It also the Creative Commons Attribution License, which shares similarities with the formation of other related structures such as colloidosomes (Meng et al., permits unrestricted use and 2014; Dinsmore et al., 2002; Manoharan et al., 2003), carboxysomes (Perlmutter et al., 2016; redistribution provided that the Garcia-Alles et al., 2017) or clathrin-coated pits (Mashl and Bruinsma, 1998; Kohyama et al., original author and source are 2003; Giani et al., 2016). Capsid formation occurs via a nucleation process driven by the favorable credited. binding energy between capsid proteins (Zandi et al., 2006). At the right assembly conditions, Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 1 of 22 Research article Physics of Living Systems thermal fluctuations induce the formation of small partial shells that tend to redissolve unless they reach a minimum critical size. Beyond this size, the shell grows by the progressive binding of subu- nits. As growth continues, the energy penalty of the naturally curved structure, due to the inescap- able presence of the rim and the accumulation of elastic energy, can be larger than the favorable binding energy. This generates a natural self-limiting mechanism for the formation of partial shells of a finite size that do not grow until closing (Grason, 2016). In fact, there are in vitro experimental evi- dences of apparently stable partial capsids (Law-Hine et al., 2016), that seems to contradict the instability of intermediates that follows from Classical Nucleation Theory (CNT) (Zandi et al., 2006). Recently, there has been a lot of interest in geometric frustration and crystal growth on spherical templates (Zandi et al., 2004; Luque et al., 2012; Meng et al., 2014; Gómez et al., 2015; Gra- son, 2016; Azadi and Grason, 2016; Paquay et al., 2017; Li et al., 2018; Panahandeh et al., 2018). Most previous works have focused either on templated growth on the surface of a sphere (Zandi et al., 2004; Luque et al., 2012; Meng et al., 2014; Gómez et al., 2015; Grason, 2016; Li et al., 2018; Panahandeh et al., 2018) or on analyzing the optimal shape of the resulting shell from pure elastic considerations (Paquay et al., 2017; Lidmar et al., 2003; Schneider and Gomp- per, 2007; Morozov and Bruinsma, 2010; Castelnovo, 2017), ignoring the importance of the deli- cate interplay of other ingredients such as the line tension, the chemical potential or the preferred curvature on their global stability and their process of formation. Here, we analyze the conditions and mechanisms leading to mis-assembly of empty viral capsids by elastic frustration, taking into account all these ingredients. We find that the outcome of the assembly depends on three scaled parameters that can be properly tuned to trigger the formation of non-spherical and open shells. Theoretical predictions obtained with the use of Classical Nucle- ation Theory including elastic contributions are confirmed qualitatively using Brownian Dynamics sim- ulations of a simple coarse-grained model. The results of this work help to better understand viral assembly and might have important implications in: envisaging novel routes to stop viral infections by interfering with their proper assembly; determining the optimal conditions for the assembly of protein cages with the desired geometry and properties for nanotechnological applications (Douglas and Young, 2006); and justifying the potential presence of seemingly stable intermediates that have been observed in recent experiments (Law-Hine et al., 2016). Results Self-assembly of a curved elastic shell The continuous description of the assembly of empty spherical viral capsids is based on Classical Nucleation Theory (CNT) (Zandi et al., 2006). In its standard version, the free energy of formation of a partial shell of area S is seen as the competition of an energy gain driving the assembly, and a rim energy penalty, due to the missing contacts at the edge of the shell. Due to the curvature of the shell and the existence of a preferred angle of interaction between capsid proteins there is another ingredient that has to be considered in the energetics of capsid formation: the elastic energy. Accordingly, the free energy of formation of a partial capsid of area S can be modeled as D DG ¼ S þ LlðSÞ þ Ge : (1) a1 The first term represents the gain in free energy associated with the chemical potential difference D between subunits in solution and in the capsid, being a1 the area per subunit. (With this defini- tion, a positive D is required to promote assembly). The second term is the total line energy of the rim, given by the product of the line tension L times its length lðSÞ. Finally, the third term Ge ¼ Gs þ Gb is the elastic energy associated with the in-plane stress, Gs , and the bending, Gb , ener- gies introduced by the curvature of the shell. Both elastic terms will be modeled using continuum elasticity theory. For the bending energy we will use the generalization of Helfrich’s model for sys- tems with non-zero spontaneous curvature introduced recently by Castelnovo (2017) (see the Appendix). For the in-plane elastic energy, we will use results from continuum elasticity theory for small deformations of thin plates, building up on recent work on the formation and growth of crystal domains of different shapes on curved surfaces (Lidmar et al., 2003; Meng et al., 2014; Seung and Nelson, 1988; Morozov and Bruinsma, 2010; Grason, 2016; Paquay et al., 2017; Köhler et al., Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 2 of 22 Research article Physics of Living Systems 2016; Schneider and Gompper, 2007; Majidi and Fearing, 2008; Castelnovo, 2017). Both stretch- ing and bending terms depend on the particular structure of the growing shell. Four different cases will be analyzed: hexagonally-ordered spherical cap without defects; spherical cap with one or many defects; ribbon and cylinder (see Figure 1). The reason to consider these particular structures is that they represent the most advantageous shapes to release the unfavorable elastic energy. In addition, cylindrical shells also appear frequently as outcome of in vitro assembly experiments. However, it is important to stress that the considered structures do not form a complete set of deformations. The relative importance of stretching versus bending contributions is controlled by a single dimensionless parameter: the Föppl-von Kárman number (FvK) defined here as g  YR20 =k, where Y is the two-dimensional Young’s modulus, R0 is the spontaneous radius of curvature and k is the bending modulus. Most previous studies have focused on the elastic energy and growth of crystals on top of a spherical template of fixed radius R. This case resembles the bending-dominated regime discussed below. Bending-dominated regime In the limit g ¼ YR20 =k  1, the bending energy dominates over the stretching energy and thus, all structures will adopt their spontaneous curvature, R ¼ R0 . The situation will be similar to the growth of a crystal on a template of fixed curvature. In the bending-dominated regime, the free energy of formation of all these structures, when properly scaled by the characteristic elastic energy 4pR20 Y, only depends on two parameters: the scaled chemical potential D~  D=ðYa1 Þ and the scaled line tension l  L=ðR0 YÞ. Thus, it is possible to compare them and determine the most stable structure for a given set of conditions. The comparison is performed for different shapes having the same area S , that is having the same number of subunits. The scaled free energy of formation of a hexagonally-ordered spherical cap of radius R0 without defects made of a circular patch of radius 0 (see Figure 1a) is D~ 2 l 1 6 Dgcap ¼ x þ xþ x (2) 4 2 1536 where x  0 =R0 is the scaled patch size, and the third term is the in-plane elastic energy of a circular domain on a curved spherical surface (Schneider and Gompper, 2007; Meng et al., 2014; Morozov and Bruinsma, 2010). (Equation 2 is an approximation strictly valid for small circular patches with an aperture angle   p, since it is assumed that the perimeter of the shell is approxi- mately the same as that of a circular disk, and a flat metric has been used to compute the in-plane elastic energy. However, we have found that a more accurate evaluation of the second and third terms in this equation [Li et al., 2018] does not alter significantly the main results.) The stretching energy stored in the spherical shell grows fast with the area of the patch, and can be partially released by two different mechanisms: by the introduction of pentagonal defects (see Figure 1b), or by growing anisotropically forming curved ribbon-like crystalline domains (see Figure 1c). The free energy of formation for a spherical cap with one defect is (Morozov and Bruinsma, 2010; Castelnovo, 2017) x2   3 2 Dgd1 ¼ Dgcap þ 1 x ; (3) 1152 2 where the last term is the stretching energy due to a pentagonal disclination at the center of the cap. (The energy of an incomplete cap with one defect placed at an arbitrary location is calculated in Li et al., 2018. It is found that the Gaussian curvature attracts the disclination to the center of the cap while the defect self-energy pushes it towards the boundary. The net result is that the minimum energy corresponds to the defect located off the center of the cap. However, we have numerically verified that this approximation introduces only a very small error in our calculations for the scaled energy. This means that not noticeable effect is observed when the exact expression with the off- center defect is considered.) Such mechanism is energetically favorable only if the second term of pffiffiffiffiffiffiffiffi Equation 3 is negative, that is, if x  2=3. Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 3 of 22 Research article Physics of Living Systems Figure 1. Sketch of the different structures considered in this study. (a) A hexagonally-ordered spherical cap of radius R and geodesic radius 0 without defects; (b) a spherical cap with a single disclination at the center (as shown) or multiple disclinations; (c) a rectangular ribbon with length L and width W, that it is called belt when the length becomes L ¼ 2pR; and (d) a cylindrical patch with size L, that eventually becomes a cylinder of radius R. In the bending-dominated regime, R ¼ R0 . For larger shells, the elastic strain is further released by the introduction of additional disclina- tions. The free energy of formation of a spherical shell with n 5-fold disclinations is (Grason, 2012; Grason, 2016; Castelnovo, 2017) Dgdn ¼ Dgcap þ gs1 þ gs2 (4) where gs1 is the self-energy of the isolated disclinations, and gs2 is their pairwise interaction, whose specific expressions are provided in the Appendix. When more than one defect appears, the mini- mum of the free energy typically occurs for a closed shell. An alternative mechanism to alleviate stretching is the anisotropic growth of the originally spheri- cal cap to adopt the shape of a defect-free rectangular curved stripe or ribbon. The free energy of formation of a ribbon of scaled length l  L=R0 , width w  W=R0 , and area s ¼ lw ¼ px2 growing on the surface of a sphere of radius R0 is (Schneider and Gompper, 2007; Majidi and Fearing, 2008) Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 4 of 22 Research article Physics of Living Systems D~ 2 l x2 l 9 Dgrib ¼ x þ þ wþ x2 w4 : (5) 4 2 w 2p 20480 Unlike the spherical cap, as the area of the patch increases, the ribbon grows longitudinally with- out limitation at a nearly fixed optimal width up to the point where l ¼ 2p, where it forms a closed belt with energy D~ 2 9 Dgbelt ¼ x þlþ x10 : (6) 4 327680 The ribbon-like structure with the lowest energy is always a closed belt rather than the open rib- bon, so we will focus our comparison with this structure. Finally, an alternative to the curved belt could be a cylinder with one principal radius of curvature infinitely large and the other R0 (see Figure 1d). The cylinder has the advantage of not having any in-plane stretching cost, but it has a bending energy penalty that prevents its formation in the bend- ing-dominated limit (see the Appendix). Figure 2 shows a comparison of the energy landscape for the different structures for fixed values of D~ and l. The competition between the bulk energy gain, the line tension penalty and the stretch- ing and bending costs will give rise, at the proper conditions, to a barrier that has to be overcome for triggering the formation of these structures. The height of this nucleation barrier and its location, corresponding to the critical cluster size, are mostly controlled by the bulk and line energy contribu- tions, since the critical size typically occurs at small values of x. In terms of shell nucleation, the bar- rier for the formation of a spherical cap is always the smallest, since the line energy of a circular edge is always smaller than for a rectangular stripe of the same area. Accordingly, the initial embryo of all these structures will be a small spherical cap (Paquay et al., 2017). Neglecting the elastic terms, the critical size for the formation of a spherical shell will be x » l=D~, corresponding to a bar- rier height for nucleation of Dgcap » l2 =ð4DÞ. ~ But rather than on the critical cluster for shell formation, we will be mostly interested in what is the most stable final structure for a given set of conditions. Since the free energies of formation only depend on l and D~ we can draw a universal phase dia- gram describing what is the structure (i.e. cap with or without defects, ribbon, or belt) with the low- est free energy in its stable size in terms of these two parameters. The term universal is intended to mean that the phase diagram is independent of the details of the capsomer-capsomer interactions such as range, preferred angle between capsomeres, bending rigidity, etc, as we corroborate with a coarse-grained simulation in the next section. Figure 3a shows the phase diagram in the bending- dominated limit, corresponding to g ¼ 0. As can be seen, belts are the most stable structure at low line tension l and chemical potential differences D. ~ Closed shells with disclinations are the pre- ferred structure for large values of D~ or l. The frontier between the belt zone and the cap with dis- clinations is approximately independent of l and located at D~ ’ 0:0020. Additionally, a small triangular region where the most stable structure is a frustrated cap with only one disclination is also apparent. As shown in the Appendix, a stable defectless cap only appears as metastable structure, since it has always a larger energy than a belt, and it is thus non competitive as stable structure, even though it may have lower energies as intermediate in the assembly process. General case of arbitrary FvK number Most small viral shells form without any underlying spherical template fixing their curvature. There- fore, it is very interesting to analyze shell formation at arbitrary FvK number, beyond the bending- dominated limit, and without the aid of an auxiliar template. In this general case, we have to con- sider the bending energy and the fact that the radius of the structures, R, may deviate from the spontaneous one, R0 , since it would be now dictated by the competition between stretching, bend- ing, and rim energies. Using the expressions for the bending energy of a sphere and a cylinder of radius R (see the Appendix), the free energy of formation of all structures analyzed in the previous section can be derived. Explicitly, the free energy of formation of a defectless spherical cap of radius R becomes Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 5 of 22 Research article Physics of Living Systems Figure 2. Comparison of free energy landscapes for different structures. Free energy of formation Dg versus the radius of the patch x in the bending-dominated regime for a defectless spherical shell (blue line), a spherical shell with defects (green), and a ribbon/belt (orange) for l ¼ 0:0001 and D~ ¼ 0:001. The optimal structure is the one with the minimum energy, which is the belt in this case. The dashed line represents the unfrustrated decrease of energy expected by the classical nucleation picture for the defectless spherical cap in the absence of elastic stresses. The inset zooms the nucleation barrier located at small patch sizes. 2 D~ 2 l 1 x6 1 2 1  Dgcap ðgÞ ¼ x þ xþ þ x 1 ; (7) 4 2 1536 r 4 4g r where r  R=R0 , and the optimal radius of the shell is given by g 4 r 2 ðr 1Þ x ¼ 0: (8) 192 Deviations from the spontaneous radius (i.e r ¼ 1) are only expected for large domain sizes or large FvK numbers. As the domain size increases, it becomes more favorable to release the elastic stress by the intro- duction of one or many 5-fold disclinations. The free energy of formation of a spherical shell with one central defect is x2 3 x2   Dgd1 ðgÞ ¼ Dgcap ðgÞ þ 1 ; (9) 1152 2 r2 pffiffiffiffiffiffiffiffi which becomes favorable over the defectless case when x=r  2=3. The formation energy of a spherical shell with n-defects is Dgdn ðgÞ ¼ Dgcap ðgÞ þ gs1 ðrÞ þ gs2 (10) where the specific expressions for gs1 ðrÞ and gs2 are written in the Appendix. Finally, the free ener- gies of a closed belt and a cylinder (which are the ribbon-like and cylindrical-patch-like structures with the lowest energy) are Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 6 of 22 Research article Physics of Living Systems Figure 3. Assembly phase diagrams. Phase diagrams of the most stable structures in terms of the scaled chemical potential D~  and the scaled line tension l for different values of the FvK number: a) g ¼ 0, corresponding to the bending-dominated regime, (b) g ¼ 80, (c) g ¼ 200, and d) g ¼ 250. Three possible equilibrium regions are present: belts (i.e. closed ribbons, in orange), frustrated capsids with one defect (red), closed shells with defects (green), and cylinders (purple). Additionally, a region corresponding to a metastable spherical cap without defects (blue) is shown only in (a). In the white region, the equilibrium state corresponds to disaggregated individual capsomers. 2 D~ 2 9 x10 1 2 1  Dgbelt ðgÞ ¼ x þ lr þ þ x 1 ; (11) 4 327680 r8 4g r and Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 7 of 22 Research article Physics of Living Systems 2 ! D~ 2  1 1 Dgcyl ðgÞ ¼ x þ lr þ x2 1 þ 1 ; (12) 4 8g r respectively. Remarkably, the free energy of formation of all these structures only depends on three scaled parameters: the chemical potential D, ~ the line tension l, and the FvK number g. Thus, it is possible to compare them and draw a universal phase diagram for the most stable structure in terms of these three parameters, to contrast with the scenario for the bending-dominated limit. Figure 3 shows phase diagrams for different values of the FvK number, showing the structure with the lowest free energy as a function of the normalized line tension l and chemical potential D~. For small values of FvK, that is g < ~ 100, the phase diagram is essentially the same as in the bending-dominated case. As the FvK number increases, the region where belts are formed occupy a larger domain, while the region with closed caps with disclinations reduces its size. However, the most relevant change is the appearance of a zone at D~>1=ð2gÞ, where the cylinder is the optimal structure. This region progres- sively invades the other structures as the FvK number is increased. Roughly for g ’ 250 only cylinders and belts are expected to be stable structures. This is a very important result since it shows that spherical capsids cannot be self-assembled directly as stable structures at large FvK numbers. The reason why cylinders dominate at large FvK numbers, corresponding to the regime where stretching dominates over bending, is because they have the advantage of not having any stretching energy cost (i.e. a flat sheet of hexamers can be bent into a cylinder without any stretching). A cylin- drical structure having a radius equal to the spontaneous radius R0 , that is r ¼ 1, will minimize the bending penalty and will have a free energy of formation, according to Equation 12, that decreases unboundedly with size when D~>1=ð2gÞ. In other words, once the formation of a cylinder becomes more favorable than free capsomers, it will continue growing without limit decreasing indefinitely its free energy of formation without paying any stretching cost, thus overcoming the energetic gain of any finite sized structure. This will be the case when D>1=ð2gÞ. ~ The larger the g (FvK), the smaller the D~ required for this to occur and therefore, regions where finite sized structures where preferred start to be devoured by the region where cylinders dominate (purple regions in Figure 3). Making use of the definition of the scaled variables, the condition for the appearance of the cylindrical phase can be recast as D  a1 k=ð2R20 Þ. In other words, cylinders appear more easily (smaller D required) for larger values of R0 , in agreement with previous results by Castelnovo (2017) predicting that cyl- inders should dominate for small spontaneous curvatures (large R0 ). Simulation A minimal coarse-grained model has been recently proposed to analyze the assembly of empty viral shells (Aznar and Reguera, 2016; Aznar et al., 2018) and other protein cages (Garcia-Alles et al., 2017). The model can successfully reproduce the assembly of the lowest spherical shell structures using capsomers, that is, pentamers and hexamers, as basic assembly units. Capsomers are coarse- grained at low resolution as effective spheres and their interaction is modelled using three contribu- tions capturing the essential ingredients (see Materials and methods): a Mie-like potential describing the attraction driving the assembly and the excluded volume interaction between a pair of capsom- ers; an angular term accounting for the preferred orientation of the interaction between proteins; and a torsion term, included to distinguish the inner and outer surfaces of the capsomers, and to favor the formation of closed shells. The model has been implemented in a Brownian Dynamics simu- lation as described in Materials and methods. One of the advantages of this simple model is that the parameters of the interaction can be related to the elastic constants (Aznar and Reguera, 2016) (see Materials and methods). In terms of YR20 2 these, the three relevant parameters controlling the assembly become g  k ¼ 94nma cos2 n, pffiffi pffiffi 2 3 kB T ln c1 =c l  YRL0 ¼ 2 cos n ~  aD nm , and D 1 Y ¼ 2 3 D pnm 0 ¼ pnm 0 : Thus, by changing the parameters of the model (mainly the exponents n and m controlling the range of the interaction, the preferred angle of interaction between capsomers n, the local bending rigitidy a, and the concentration c1 which controls the effective chemical potential D) we can explore the universality and the different scenarios of assembly discussed in the previous section. Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 8 of 22 Research article Physics of Living Systems Figure 4 shows the results of simulations using different sets of parameters represented in scaled units and contrasted with the theoretical phase diagram for g ¼ 80. For l ¼ 0:00084, that is a rela- tively large line tension, at low concentration of capsomers, the seed dissolves and no nucleation or growth occurs. As the capsomer concentration is progressively increased, a metastable defectless shell and a closed spherical shell with typically 12 defects form, as expected by the theory. At very high concentrations, nucleation occurs simultaneously at many sites, and the final outcome of the Figure 4. Comparison of simulation results with the theoretical phase diagram. (a) Phase diagram of the most stable structure in terms of the scaled chemical potential D~  and the scaled line tension l for g ¼ 80. Snapshots of the final outcome of the simulation for different initial capsomer concentrations for: b) l ¼ 0:00037, obtained with n ¼ 1:45 and a m ¼ 36, n ¼ 18 potential. By increasing the concentration of capsomers one goes from a dissassembled state to the formation of a ribbon to a spherical shell with defects. (The ribbon and spherical shell are not closed in the snapshots due to the limited number of capsomers and finite simulation time). (c) l ¼ 0:000525, obtained with n ¼ 1:38 and a m ¼ 36, n ¼ 18 potential. As concentration increases, one goes from a shell with one defect to a complete shell with many defects. By increasing the FvK number to g ¼ 1280 (by setting a ¼ 0:4), the simulations clearly forms cylindrical tubes. (d) l ¼ 0:00084, obtained with n ¼ 1:45 and a m ¼ 24, n ¼ 12 potential. In this case, the sequence is: disagreggated, metastable defectless shell, and closed spherical shell with many defects (a partial shell is shown). The simulation results agree qualitatively with the predictions from the theoretical phase diagram. Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 9 of 22 Research article Physics of Living Systems simulation are many fragments of spherical capsids that cannot grow any further due to the deple- tion of free capsomers in solution. That would correspond to kinetic trapping, which is an interesting alternative mechanism to prevent the correct capsid assembly, that will be analyzed in a future work. For l ¼ 0:000585, as the concentration is increased we obtained the expected sequence of: seed dis- olution; formation of a stable cap with a single 5-fold defect in a very narrow range of concentra- tions; and the formation of closed caps with many defects. Finally for l ¼ 0:000372, that is a relatively small line tension, as the capsomer concentration is increased we go from no assembly, to the for- mation of ribbon-like stripes, to the growth of spherical shells with many defects. As naively expected, as l increases, higher scaled chemical potentials are needed to nucleate the structures. Finally, by increasing the parameter a, the bending rigidity is reduced and assembly at higher FvK can be analyzed. The results of the simulations show that as the FvK number is increased, the forma- tion of spherical shells is overriden by the formation of cylindrical bodies, as shown in Figure 4c, that also competes with other elongated structures such as spherocylinders or even conical shapes (see Appendix 1—figure 3). Remarkably, simulations that have been performed for widely different val- ues of the interaction parameters, when properly scaled, all fall into the predicted picture. Therefore, the simulation results nicely confirm almost quantitatively the universality of the fate of the assembly and the potential scenarios discussed in the theory. A precise quantitative comparison between the theory and the simulations has not been performed, since they are done at slightly different condi- tions. While theory assumes a reservoir of capsomers, the simulations are done at fixed total number of subunits. This implies that, as the assembly proceeds, the concentration of the remaining free par- ticles, and consequently the chemical potential, decreases. For this reason, we have not intended to reproduce with precision the borders of the phase diagram using the simulations. The fact that the chemical potential is not strictly constant in the simulations due to the depletion of free subunits may cause quantitative discrepancies when comparing with theory, but does not alter the relative stability of the different shapes analyzed. Discussion We have provided a comprehensive analysis of non-templated assembly of curved elastic shells, tak- ing into account all relevant ingredients (i.e. chemical potential, line tension, spontaneous curvature, and elastic contributions) and the potential formation of non-spherical shapes. The importance of accounting for all these ingredients becomes evident, for instance, in the study of the stability of the defectless spherical cap, which turns out to be always metastable, its global stability hindered by the introduction of defects (at high line tensions) or the formation of ribbons (at low line tensions). Our analysis also shows that the outcome of the assembly not only depends on elastic considerations, but also on the assembly conditions, represented here by the scaled chemical potential. Hence, either belts or closed spherical shells or cylinders may be obtained as the most stable structure for fixed interaction parameters, depending on the concentration of assembly units. When assembly takes place at conditions near the vicinity of a phase boundary, a mixture of the two phases, or a structure resulting from their combination (e.g. a spherocylinder) may form. This may justify the observation of coexisting tubes and spherical capsids in the in vitro assembly of viruses such as SV40 (Kanesashi et al., 2003). Although, for the sake of simplicity, our theoretical analysis has been performed using the contin- uous and small curvature approximations, we have verified that releasing these approximations does not alter significantly the results. The exact expression of the perimeter of the growing edge (Zandi et al., 2006; Gómez et al., 2015) influences the height and location of the nucleation barrier, but has a minor impact on the properties of the final stable structure. The accurate evaluation of the in-plane elastic cost of defects taking into account their spatial distribution (Li et al., 2018), modifies the energies of the growing shell, but does not modify significantly the stability of the final structure. Simulations of a coarse-grained model made using widely different values for the parameters and interaction range confirm that the outcome of the assembly only depends on three scaled parame- ters: the scaled chemical potential D,~ line tension l, and FvK number g. Thus, the assembly phase diagram is universal, and different protein shells, interaction potentials and coarse-grained models can be recast into a unifying picture of assembly, that could guide the efficient production of artifi- cial viral cages. For instance, our analysis indicates that relatively long-range interactions are desir- able to increase the line tension, decrease the FvK number and facilitate the assembly of closed Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 10 of 22 Research article Physics of Living Systems spherical shells. In fact, spherical shells with icosahedral symmetry and triangulation number T > 7 could be successfully assembled in simulations without any template or scaffolding protein, provided that the line tension and FvK number are adequate. Alternatively, chemical or physical modifications that increase the FvK number or reduce the line tension or the effective concentration may become a potential therapeutic target to prevent viral replication by inducing the formation of open, and presumably non-infective, cylindrical or belt-like structures. Experimentally, the chemical potential can be tuned by the total protein concentration or by the addition of crowding agents. The line ten- sion (which depends on the strength of the binding interaction), could be modified by the tempera- ture, the pH and the salt concentration. The bending rigidity and spontaneous radius of curvature are also presumably controlled by pH and the presence, concentration and nature of ions or auxiliary proteins in solution. Further experimental and theoretical investigations are required to make a pre- cise quantitative connection between the physical parameters controlling the assembly and experiments. Triggering the formation of closed spherical shells with an incorrect radius, triangulation number (Caspar and Klug, 1962), or arrangement of proteins could also be an alternative to interfere with the assembly of the right viral capsid. But in our study, we have focused on mechanisms interferring with the closing of the shell by elastic frustration, rather than classifying the specific radius and trian- gulation number of the resulting spherical structure. In addition, we do not consider the situation in which the capsomers interact with cargo. Such interactions are crucial for viruses that co-assemble with their genetic material or a cargo, but this is beyond the scope of the present study. A very important conclusion of our analysis is that spherical capsids cannot be self-assembled directly as stable structures at large FvK numbers. This may explain why some viruses that require high mechanical resistance, for instance many dsDNA bacteriophages such as lambda, HK97 and P22, first assemble a relatively soft spherical procapsid before suffering a maturation transition (Roos et al., 2012; Johnson, 2010) that flattens out their faces, which is a clear signature of a high FvK number (Lidmar et al., 2003). The results of our work indicate that a one-step assembly of a spherical shell with the high elastic resistance and Fvk number of the final structure is not viable. Table 1 compares the estimated elastic properties of different empty capsids of real viruses. The table clearly shows that viruses like CCMV or SV40 that assemble easily in vitro as spherical shells, have estimated values of the scaled line tension and FvK in the region where these structures are expected to be stable outcomes of the assembly. Contrarily, the high FvK number of the mature bacteriophage lambda will prevent its direct assembly. However, its procapsid, which is the first structure that is assembled, has a larger scaled line tension and smaller FvK that would facilitate a successful assembly. (The FvK number of lambda procapsid listed in Table 1 is probably overesti- mated, given its noticeable spherical shell. In addition, we have found in our simulations that even though the theoretical threshold for the disappearance of spherical shells as stable structures is around g ¼ 250, in practice larger FvK numbers are typically required to obtain cylindrical structures since the nucleation barrier for their formation is larger than for the metastable spherical shell). Table 1. Estimates of the main geometric and elastic properties of different non-enveloped empty viral capsids. The Young’s modulus E has been evaluated from AFM nanoindentation experiments (Mateu, 2012; Michel et al., 2006; Ivanovska et al., 2007; Sae-Ueng et al., 2014) and, for SV40, from the experimental spring constant (van Rosmalen et al., 2018) using the standard thin shell formula k ¼ 2:25Eh2 =R (Ivanovska et al., 2004). The 2D Young’s Modulus was calculated as Y ¼ Eh; the pffiffiffi qffiffiffiffiffiffi  effective diameter of the capsomers as (Santolaria, 2011) s ¼ R= 5 p 3 T þ p1ffiffi3 cot p5  1 , where T is the triangulation number; the line tension as (Luque et al., 2012) l ¼ p2ffiffi30s considering a typical binding energy 0 » 10kB T; and the FvK number as g ¼ 12ð1 n2p ÞðR=hÞ2 , with np ¼ 0:3 (Buenemann and Lenz, 2008). Virus T-number Diameter (nm) Thickness h (nm) E (Gpa) Y (N/m) s (nm) Scaled line tension l Föppl-von Karman g CCMV 3 28 3.8 0.14 0.53 5.9 0.00107 148 l Procapsid 7 50 4.0 0.16 0.64 6.8 0.000436 427 l Capsid 7 63 1.8 1.0 1.8 8.6 0.0000976 3344 SV40 7 45 6.0 0.033 0.2 6.1 0.00174 152 Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 11 of 22 Research article Physics of Living Systems In summary, we have seen that the fate of the assembly is controlled by a universal phase diagram in terms of three scaled parameters: line tension, chemical potential and FvK number. The phase dia- grams shed light on the physics controlling the assembly of curved shells, and could guide assembly experiments to achieve either an efficient assembly of artificial viral shells of desired geometry and mechanical properties or, alternatively, to envisage the conditions needed to impede viral infections by arresting viral assembly or inducing missasembly into a non-infective structure. Materials and methods Coarse-grained model and simulation details The simulation model, introduced in Aznar and Reguera (2016); Aznar et al. (2018), is coarse- grained at the level of capsomers which are represented as effective spheres of two different diame- ters: sh and sp , reflecting the fact that hexamers and pentamers are made of a different number of proteins (six and five, respectively). The interaction between capsomers, V ¼ VLJ  Va  Vtor , is mod- eled using three contributions: a Mie-like, an angular, and a torsion potential. The Mie-like potential n hsij m m sij n i VLJ ðrij Þ ¼ ij ; (13) m n r n r describes the binding and the excluded volume interaction between a pair of capsomers in terms of their relative distance, rij is the equilibrium distance corresponding to the minimum of the potential, r is the distance between capsomers centers, ij is the binding energy between capsomers, and m and n represent the power of the repulsive and attractive interaction terms, respectively, which set the range of the interaction potential. The angular contribution is given by ! ! ðij nÞ2 ðji nÞ2 Va ðrij ; Wi ; Wj Þ ¼ exp exp ; (14) 2a2 2a2 where ij is the angle between the vector Wi , describing the spatial orientation of the capsomer, and the unit vector rij . The parameter n is the preferred angle of interaction between proteins of differ- ent capsomers, and the parameter a controls the local bending stiffness, that is, the energy cost required to bend two capsomers out of their preferred angle of interaction. The torsion term is given by   ð1 cos Þ Vtor ðWi ; Wj Þ ¼ exp kt ; (15) 2 where kt is the torsion constant and  is the angle between the planes defined by the unit vector rij and both orientation vectors. The elastic properties of a shell can be related to the main parameters of the interaction. In par- pffiffi ticular the Young’s modulus is approximately given by Y ¼ 2nm pffiffi 02 ; the bending rigidity is k ¼ 3 3 02 ; s 3 8 a s and the preferred radius of curvature is R0 ¼ 2 cos n : The line tension of a partially formed cap can be approximated by Luque et al. (2012) l ¼ p2ffiffi30s ; and the chemical potential difference that controls the assembly is given by D ¼ kB T lnðc1 =c Þ; where c1 is the concentration of free capsomers and c is the critical concentration (Zandi et al., 2006). The model has been implemented in a Brownian Dynamics simulation code using a simple sto- chastic Euler’s integration algorithm, as described in Aznar et al. (2018). Simulations were made with only one type of capsomers. We worked using reduced units in terms of the diameter of the basic building blocks s, their diffusion coefficient D, and the binding energy 0 . In these reduced units, the parameters used in the simulation are: torsion constant kt ¼ 1:5, reduced temperature T ¼ 0:1, corresponding to a binding energy between capsomers of 10 kB T, representing the typical order of magnitude of the strength of interactions between viral capsid proteins. Since, in all cases, the critical nucleus is a partial spherical cap, all simulations were started using a small spherical cap of 19 units with the spontaneous curvature as initial seed. The remaining capsom- ers up to a total of N ¼ 200 400 were initially placed randomly inside a cubic box with periodic boundary conditions. The simulations run for a total of 2  109 steps and the final structures were Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 12 of 22 Research article Physics of Living Systems analyzed. To verify the universality of the phase diagram, we performed an extensive set of simula- tions with different interaction parameters. More specifically, the interaction range was varied from m ¼ 12; n ¼ 6 to m ¼ 48; n ¼ 24; the spontaneous angle in the range 1:24 D~ms , the free energy landscape has a maximum, signaling the nucleation barrier, but also a local minimum, xmin, corresponding to a locally-stable spherical cap. Thus, unlike in the standard case where beyond the critical size the free energy goes steadily down and the shell can grow until closing, the high elastic cost associated with the curvature of the shell will prevent its further growth and force it to reach an equilibrium size (Grason, 2016). The condition D~  D~ms marks the onset of a metastable region where the value of free energy Dgcap ðxmin Þ of the locally-stable shell is larger than its value for the dissembled state at x = 0. A fully stable partial shell is obtained if Dgcap ðxmin Þ  0 at the minimum. The onset of the stable region can be obtained from the conditions Dgcap ¼ 0 and dDgcap =dx ¼ 0, which allows to find the critical value  1=5 5 l4 D~scap ¼ ; (A7) 4 6 beyond which Dgcap ðxmin Þ  0. The corresponding value of x is xs ¼ 2ð6lÞ1=5 ; (A8) and represents the minimum size of a partial spherical shell without defects to be stable. The size of stable partial shells grows as D~ increases. For large patches or when g  1, the rim energy can be neglected and the location of the stable cap size is described approximately ~ 1=4 . by xeq » ð128DÞ Appendix 1—figure 1. Free energy of formation Dgcap of a spherical cap without defects, Equation A4, versus the radius of the patch x for l ¼ 0:0001 and different values of the scaled chemical potential D, ~ illustrating the situations in which no assembly is possible (short dashed line, D~ ¼ 0:0002), a geometrically frustrated metastastable cap (dashed line, D~ ¼ 0:0005) or stable (solid line, D~ ¼ 0:001) finite shell. Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 17 of 22 Research article Physics of Living Systems When the area of the spherical crystalline patch gets large, it becomes favorable to introduce defects to release the elastic stress. The in-plane elastic energy of a spherical cap with one 5-fold disclination is (Morozov and Bruinsma, 2010; Grason, 2016) pY20 3 20   DGsd1 ¼ 1 : (A9) 288 2 R20 Thus, the free energy of formation, in scaled units, for a spherical cap with one defect at the center becomes x2   3 2 Dgd1 ¼ Dgcap þ 1 x : (A10) 1152 2 Such mechanism is energetically favorable only if the second term of Equation A10 is pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi negative, that is, if x  2=3. The fact that for x  2=3 a shell with a defect becomes energetically favorable imposes an important restriction on the values that the scaled chemical potential and line tension may have in order to allow the existence of a stable defectless pffiffiffiffiffiffiffiffi spherical cap. Imposing xs  2=3 one gets the requirements l  67=2 and D~  5=864. Using the conditions Dgd1 ¼ 0 and dDgd1 =dx ¼ 0, we can find the onset of stability for a circular cap with one central disclination, described by 2l 4 6x21 þ 3x41 D~sd1 ¼ þ ; (A11) x1 1152 with  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1=2 1 þ 1 þ 768l x1 ¼ : (A12) 2 For larger shells, the elastic strain is further released by the introduction of additional disclinations. The resulting free energy of formation of a spherical shell with n 5-fold disclinations in scaled units is (Grason, 2012; Grason, 2016; Castelnovo, 2017) Dgdn ¼ Dgcap þ gs1 þ gs2 (A13) where  n  2 x2 x2a  3 2 X gs1 ¼ 1 x 1 ; (A14) 1152 2 a¼1 x2 is the self-energy of the isolated disclinations, and x2 X  gs2 ¼ Vint xa ;xb ; (A15) 1152 b6¼a is the pairwise interaction of disclinations, with 8 2 39 2 2 ! x2a x2b    < xa xb xa xb = Vint xa ;xb ¼ 1 1 þ ln4   5 ; (A16) : x2 x2 x2 x2 x2a x2 2 xb =x2 þ xa xb 2 ; being xa  a =0 , a the geodesical position of the disclinations and xa xb is the (normalized) geodesic distance between disclination a and b. When more than one defect appear, the minimum of the free energy typically occurs for a closed shell, corresponding in the small curvature approximation to x ’ 2 in which case, the stability region appears for D~ larger than D~sdn ’ l þ 0:00146: (A17) and the cap consists of a fully closed shell with defects. Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 18 of 22 Research article Physics of Living Systems The free energy of formation of a ribbon of length L and width W, growing on the surface of a sphere of radius R0 is (Schneider and Gompper, 2007; Majidi and Fearing, 2008) LW 9Y W 5 L DGrib ¼ D þ 2LðL þ WÞ þ : (A18) a1 5120 R40 The energetic advantage of this configuration over the spherical cap is that, for a fixed width W, the in-plane elastic energy only grows linearly with length. In order to compare the energy of ribbons made with similar number of units as the spherical cap, we consider that both structures have the same area. That is, S ¼ p20 ¼ LW. Thus, this energy can be rewritten in dimensionless terms as D~ 2 l x2 l 9 Dgrib ¼ x þ þ wþ x2 w4 (A19) 4 2 w 2p 20480 where w  W=R0 . The optimal width of the ribbon is obtained by minimization that is, from qDgrib =qw ¼ 0, yielding x2 l 9 2 3 l þ þ x w ¼ 0: (A20) w2 p 2560 As the area of the patch increases, the ribbon grows longitudinally at a nearly fixed optimal width up to the point where l ¼ 2p, where it forms a closed belt with energy D~ 2 9 Dgbelt ¼ x þlþ x10 : (A21) 4 327680 The elastic energy of the belt grows very steeply with the patch size, since after closing the ribbon, it can only grow by increasing its width at a large stretching cost. It can be shown that the equilibrium ribbon-like structure with the lowest energy is always a closed belt rather than the open ribbon. So the competing structures are the spherical cap with or without defects and the belt. The minimum in the free energy, corresponding to a stable belt, is located at 64D~ 1=8   xsbelt ¼ 2 ; (A22) 9 and the stability region for the belt occurs when D~ is larger than rffiffiffi !4=5 5 3 D~sbelt ¼ l : (A23) 8 2 Since D~sbelt 1=ð2gÞ for its formation, meaning that cylindrical shells cannot be formed in the bending-dominated limit, corresponding to g ! 0. Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 19 of 22 Research article Physics of Living Systems General case of arbitrary FvK The description of the free energy of formation of shells at arbitrary FvK numbers involves the consideration of the bending energy and of the radius of the structures as an additional free parameter that may now deviate from the spontatneous radius R0 . Particularizing Equation A1, the bending energy of a sphere of radius R in scaled units is  2 sph 1 2 1 Dgb ¼ x 1 ; (A26) 4g r while for a cylinder of radius R reads  2 ! 1 1 Dgcyl b ¼ x2 1 þ 1 (A27) 8g r where r  R=R0 and the Föppl-von Kárman number (FvK), quantifying the ratio of stretching and bending energies, is still defined as g  YR20 =k. Using these expressions we can generalize the free energy of formation of all structures analyzed in the previous section. Explicitly, in reduced units, the free energy of formation of a defectless spherical cap of radius r becomes 2 D~ 2 l 1 x6 1 2 1  Dgcap ðgÞ ¼ x þ xþ þ x 1 : (A28) 4 2 1536 r 4 4g r In the general case, the optimal radius of the shell is obtained from the condition qDgcap =qr ¼ 0, yielding g 4 r 2 ðr 1Þ x ¼ 0: (A29) 192 This equation shows that deviations from the preferred radius (i.e r ¼ 1) are only expected for large domain sizes or large FvK numbers. The corresponding solution with positive second derivative can be obtained analytically, although its expression is a bit cumbersome. For small shells or FvK numbers, it can be well approximated by its two leading terms, yielding g 4 r ’ 1þ x : (A30) 192 In the limit of large shells or Fvk numbers, the radius goes as 1 1 g 4 1=3 r’ þ x : (A31) 3 4 3 Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 20 of 22 Research article Physics of Living Systems Appendix 1—figure 2. Scaled optimal radius r of a spherical cap without defects as a function of the patch size x for a FvK number g ¼ 80. The solid line is the exact solution, Equation A29, the short dashed line is the approximation for large shells or FvK, Equation A31, and the dashed line is the approximation for small shells or FvK, Equation A30. In Appendix 1—figure 2, the curve r vs x is sketched as obtained from the complete solution of Equation A29 (solid line). As can be seen, small caps adopt a curvature close to the one of the bending dominated case, r ’ 1. Larger caps tend to flatten out. The dashed line shows the small shell approximation, Equation A30, while the short dashed line shows the assymptotic approximation given by Equation A31. As the domain size increases, it becomes more favorable to release the elastic stress by the introduction of one or many 5-fold disclinations. The free enery of formation of a spherical shell with one central defect is x2 3 x2   Dgd1 ðgÞ ¼ Dgcap ðgÞ þ 1 : (A32) 1152 2 r2 pffiffiffiffiffiffiffiffi which becomes favorable over the defectless case when x=r  2=3. The formation energy of a spherical shell with n-defects is Dgdn ðgÞ ¼ Dgcap ðgÞ þ gs1 ðrÞ þ gs2 (A33) where  n  2 x2 3 x2 X x2a  gs1 ðrÞ ¼ 1 1 ; (A34) 1152 2 r2 a¼1 x2 and gs2 is given by Equation A15. The free energy of a ribbon of scaled width w and length l in the general case becomes 2 D~ 2 l x2 l 9 w4 1 1 Dgrib ðgÞ ¼ x þ þ wþ x2 þ x2 1 : (A35) 4 2 w 2p 20480 r 4g r The optimal width and radius of curvature of the ribbon, w and r, are obtained by minimization that is, from qDgrib =qw ¼ 0 and qDgrib =qr ¼ 0, yielding Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 21 of 22 Research article Physics of Living Systems x2 l 9 2 w3 l þ þ x ¼0 (A36) w p 2560 r 4 2 and 9gw4 r 2 ð1 rÞ ¼ 0: (A37) 2560 Finally, the free energies of a belt and a cylinder are 2 D~ 2 9 x10 1 2 1  Dgbelt ðgÞ ¼ x þ lr þ þ x 1 (A38) 4 327680 r 8 4g r and 2 ! D~ 2  1 1 Dgcyl ðgÞ ¼ x þ lr þ x2 1 þ 1 ; (A39) 4 8g r respectively. Appendix 1—figure 3. Other structures obtained in the simulations. (a) Bullet-shaped shell and (b) conical shell obtained in two different repetitions for g ¼ 224, l ¼ 0:00069, with m ¼ 36, n ¼ 18, n ¼ 1:345, a ¼ 0:4 and  ¼ 0:02; (c) T = 13 icosahedral shell obtained for g ¼ 100, l ¼ 0:0005246, with m ¼ 36, n ¼ 18, n ¼ 1:40, a ¼ 0:1 and  ¼ 0:02; (d) coexistence between a cylinder and a partial spherocylindrical shell obtained for g ¼ 900, l ¼ 0:0005246, with m ¼ 36, n ¼ 18, n ¼ 1:40, a ¼ 0:3 and  ¼ 0:02; (e) branched ribbon-like structure (Köhler et al., 2016) obtained at low values of the scaled line tension for g ¼ 71, l ¼ 0:00031, with m ¼ 36, n ¼ 18, n ¼ 1:47, a ¼ 0:1 and  ¼ 0:0125. Mendoza and Reguera. eLife 2020;9:e52525. DOI: https://doi.org/10.7554/eLife.52525 22 of 22