Shape selection and misassembly in viral capsid formation by elastic frustration
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 selfassembly, 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 ribbonshaped shells may emerge, preventing the successful replication of viruses. Using elasticity theory and coarsegrained simulations, we analyze the conditions required for these processes to occur and their significance for empty virus selfassembly. 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 selfassembled 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 misassembly.
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 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 functions which are essential in their replication cycle. One of the most amazing one is their ability to selfassemble with an unparalleled efficiency and precision.
In vivo, the capsid of most viruses assembles from its basic building blocks, which could be individual capsid proteins, dimers, trimers or capsomers (i.e. clusters of five or six proteins, which constitute the structural and morphological units of the shell). The resulting structure has a precise architecture, which in most cases is spherical with icosahedral symmetry (Roos et al., 2010). Several viruses assemble their capsid before packaging the genetic material. In addition, the proteins of many viruses have the ability to selfassemble in vitro, even in the absence of genetic material, forming empty capsids.
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; Perlmutter and Hagan, 2015; Hagan and Zandi, 2016). The assembly of a curved empty shell with a welldefined geometry and precise arrangement of the building blocks is a nontrivial process that resembles 2D crystallization on a curved space (Meng et al., 2014; Gómez et al., 2015). It also shares similarities with the formation of other related structures such as colloidosomes (Meng et al., 2014; Dinsmore et al., 2002; Manoharan et al., 2003), carboxysomes (Perlmutter et al., 2016; GarciaAlles et al., 2017) or clathrincoated pits (Mashl and Bruinsma, 1998; Kohyama et al., 2003; Giani et al., 2016). Capsid formation occurs via a nucleation process driven by the favorable binding energy between capsid proteins (Zandi et al., 2006). At the right assembly conditions, 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 subunits. As growth continues, the energy penalty of the naturally curved structure, due to the inescapable presence of the rim and the accumulation of elastic energy, can be larger than the favorable binding energy. This generates a natural selflimiting 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 evidences of apparently stable partial capsids (LawHine 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; Grason, 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 Gompper, 2007; Morozov and Bruinsma, 2010; Castelnovo, 2017), ignoring the importance of the delicate 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 misassembly 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 nonspherical and open shells. Theoretical predictions obtained with the use of Classical Nucleation Theory including elastic contributions are confirmed qualitatively using Brownian Dynamics simulations of a simple coarsegrained 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 (LawHine et al., 2016).
Results
Selfassembly 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
The first term represents the gain in free energy associated with the chemical potential difference $\mathrm{\Delta}\mu $ between subunits in solution and in the capsid, being ${a}_{1}$ the area per subunit. (With this definition, a positive $\mathrm{\Delta}\mu $ is required to promote assembly). The second term is the total line energy of the rim, given by the product of the line tension $\mathrm{\Lambda}$ times its length $l(S)$. Finally, the third term ${G}_{e}={G}_{s}+{G}_{b}$ is the elastic energy associated with the inplane stress, ${G}_{s}$, and the bending, ${G}_{b}$, energies 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 systems with nonzero spontaneous curvature introduced recently by Castelnovo (2017) (see the Appendix). For the inplane 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., 2016; Schneider and Gompper, 2007; Majidi and Fearing, 2008; Castelnovo, 2017). Both stretching and bending terms depend on the particular structure of the growing shell. Four different cases will be analyzed: hexagonallyordered 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öpplvon Kárman number (FvK) defined here as $\gamma \equiv Y{R}_{0}^{2}/\kappa$, where $Y$ is the twodimensional Young’s modulus, ${R}_{0}$ is the spontaneous radius of curvature and $\kappa $ 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 bendingdominated regime discussed below.
Bendingdominated regime
In the limit $\gamma =Y{R}_{0}^{2}/\kappa \ll 1$, the bending energy dominates over the stretching energy and thus, all structures will adopt their spontaneous curvature, $R={R}_{0}$. The situation will be similar to the growth of a crystal on a template of fixed curvature. In the bendingdominated regime, the free energy of formation of all these structures, when properly scaled by the characteristic elastic energy $4\pi {R}_{0}^{2}Y$, only depends on two parameters: the scaled chemical potential $\mathrm{\Delta}\stackrel{~}{\mu}\equiv \mathrm{\Delta}\mu /(Y{a}_{1})$ and the scaled line tension $\lambda \equiv \mathrm{\Lambda}/({R}_{0}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 hexagonallyordered spherical cap of radius ${R}_{0}$ without defects made of a circular patch of radius ${\rho}_{0}$ (see Figure 1a) is
where $x\equiv {\rho}_{0}/{R}_{0}$ is the scaled patch size, and the third term is the inplane 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 $\theta \ll \pi$, since it is assumed that the perimeter of the shell is approximately the same as that of a circular disk, and a flat metric has been used to compute the inplane 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 ribbonlike crystalline domains (see Figure 1c).
The free energy of formation for a spherical cap with one defect is (Morozov and Bruinsma, 2010; Castelnovo, 2017)
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 selfenergy 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 offcenter defect is considered.) Such mechanism is energetically favorable only if the second term of Equation 3 is negative, that is, if $x\ge \sqrt{2/3}$.
For larger shells, the elastic strain is further released by the introduction of additional disclinations. The free energy of formation of a spherical shell with n 5fold disclinations is (Grason, 2012; Grason, 2016; Castelnovo, 2017)
where ${g}_{{s}_{1}}$ is the selfenergy of the isolated disclinations, and ${g}_{{s}_{2}}$ is their pairwise interaction, whose specific expressions are provided in the Appendix. When more than one defect appears, the minimum of the free energy typically occurs for a closed shell.
An alternative mechanism to alleviate stretching is the anisotropic growth of the originally spherical cap to adopt the shape of a defectfree rectangular curved stripe or ribbon. The free energy of formation of a ribbon of scaled length $l\equiv L/{R}_{0}$, width $w\equiv W/{R}_{0}$, and area $s=lw=\pi {x}^{2}$ growing on the surface of a sphere of radius ${R}_{0}$ is (Schneider and Gompper, 2007; Majidi and Fearing, 2008)
Unlike the spherical cap, as the area of the patch increases, the ribbon grows longitudinally without limitation at a nearly fixed optimal width up to the point where $l=2\pi $, where it forms a closed belt with energy
The ribbonlike structure with the lowest energy is always a closed belt rather than the open ribbon, 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 ${R}_{0}$ (see Figure 1d). The cylinder has the advantage of not having any inplane stretching cost, but it has a bending energy penalty that prevents its formation in the bendingdominated limit (see the Appendix).
Figure 2 shows a comparison of the energy landscape for the different structures for fixed values of $\mathrm{\Delta}\stackrel{~}{\mu}$ and $\lambda $. The competition between the bulk energy gain, the line tension penalty and the stretching 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 contributions, since the critical size typically occurs at small values of x. In terms of shell nucleation, the barrier 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}^{*}\approx \lambda /\mathrm{\Delta}\stackrel{~}{\mu}$, corresponding to a barrier height for nucleation of $\mathrm{\Delta}{g}_{cap}^{*}\approx {\lambda}^{2}/(4\mathrm{\Delta}\stackrel{~}{\mu})$. 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 $\lambda $ and $\mathrm{\Delta}\stackrel{~}{\mu}$ we can draw a universal phase diagram describing what is the structure (i.e. cap with or without defects, ribbon, or belt) with the lowest 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 capsomercapsomer interactions such as range, preferred angle between capsomeres, bending rigidity, etc, as we corroborate with a coarsegrained simulation in the next section. Figure 3a shows the phase diagram in the bendingdominated limit, corresponding to $\gamma =0$. As can be seen, belts are the most stable structure at low line tension $\lambda $ and chemical potential differences $\mathrm{\Delta}\stackrel{~}{\mu}$. Closed shells with disclinations are the preferred structure for large values of $\mathrm{\Delta}\stackrel{~}{\mu}$ or $\lambda $. The frontier between the belt zone and the cap with disclinations is approximately independent of $\lambda $ and located at $\mathrm{\Delta}\stackrel{~}{\mu}\simeq 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. Therefore, it is very interesting to analyze shell formation at arbitrary FvK number, beyond the bendingdominated limit, and without the aid of an auxiliar template. In this general case, we have to consider the bending energy and the fact that the radius of the structures, $R$, may deviate from the spontaneous one, ${R}_{0}$, since it would be now dictated by the competition between stretching, bending, 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
where $r\equiv R/{R}_{0}$, and the optimal radius of the shell is given by
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 introduction of one or many 5fold disclinations. The free energy of formation of a spherical shell with one central defect is
which becomes favorable over the defectless case when $x/r\ge \sqrt{2/3}$. The formation energy of a spherical shell with $n$defects is
where the specific expressions for ${g}_{{s}_{1}}(r)$ and ${g}_{{s}_{2}}$ are written in the Appendix. Finally, the free energies of a closed belt and a cylinder (which are the ribbonlike and cylindricalpatchlike structures with the lowest energy) are
and
respectively.
Remarkably, the free energy of formation of all these structures only depends on three scaled parameters: the chemical potential $\mathrm{\Delta}\stackrel{~}{\mu}$, the line tension $\lambda $, and the FvK number $\gamma $. 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 bendingdominated 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 $\lambda $ and chemical potential $\mathrm{\Delta}\stackrel{~}{\mu}$. For small values of FvK, that is $\gamma \lesssim 100$, the phase diagram is essentially the same as in the bendingdominated 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 $\mathrm{\Delta}\stackrel{~}{\mu}>1/(2\gamma )$, where the cylinder is the optimal structure. This region progressively invades the other structures as the FvK number is increased. Roughly for $\gamma \simeq 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 selfassembled 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 cylindrical structure having a radius equal to the spontaneous radius ${R}_{0}$, 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 $\mathrm{\Delta}\stackrel{~}{\mu}>1/(2\gamma )$. 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 $\mathrm{\Delta}\stackrel{~}{\mu}>1/(2\gamma )$. The larger the $\gamma $ (FvK), the smaller the $\mathrm{\Delta}\stackrel{~}{\mu}$ 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 $\mathrm{\Delta}\mu \ge {a}_{1}\kappa /(2{R}_{0}^{2})$. In other words, cylinders appear more easily (smaller $\mathrm{\Delta}\mu $ required) for larger values of ${R}_{0}$, in agreement with previous results by Castelnovo (2017) predicting that cylinders should dominate for small spontaneous curvatures (large ${R}_{0}$).
Simulation
A minimal coarsegrained 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 (GarciaAlles 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 coarsegrained at low resolution as effective spheres and their interaction is modelled using three contributions capturing the essential ingredients (see Materials and methods): a Mielike potential describing the attraction driving the assembly and the excluded volume interaction between a pair of capsomers; 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 simulation 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 these, the three relevant parameters controlling the assembly become $\gamma \equiv \frac{Y{R}_{0}^{2}}{\kappa}=\frac{4nm{\alpha}^{2}}{9{\mathrm{cos}}^{2}\nu}$, $\lambda \equiv \frac{\mathrm{\Lambda}}{Y{R}_{0}}=\frac{2\mathrm{cos}\nu}{nm}$, and $\mathrm{\Delta}\stackrel{~}{\mu}\equiv \frac{\mathrm{\Delta}\mu}{{a}_{1}Y}=\frac{2\sqrt{3}}{\pi nm}\frac{\mathrm{\Delta}\mu}{{\u03f5}_{0}}=\frac{2\sqrt{3}}{\pi nm}\frac{{k}_{B}T\mathrm{ln}{c}_{1}/{c}^{*}}{{\u03f5}_{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 $\nu $, the local bending rigitidy $\alpha $, and the concentration ${c}_{1}$ which controls the effective chemical potential $\mathrm{\Delta}\mu $) we can explore the universality and the different scenarios of assembly discussed in the previous section.
Figure 4 shows the results of simulations using different sets of parameters represented in scaled units and contrasted with the theoretical phase diagram for $\gamma =80$. For $\lambda =0.00084$, that is a relatively 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 simulation are many fragments of spherical capsids that cannot grow any further due to the depletion 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 $\lambda =0.000585$, as the concentration is increased we obtained the expected sequence of: seed disolution; formation of a stable cap with a single 5fold defect in a very narrow range of concentrations; and the formation of closed caps with many defects. Finally for $\lambda =0.000372$, that is a relatively small line tension, as the capsomer concentration is increased we go from no assembly, to the formation of ribbonlike stripes, to the growth of spherical shells with many defects. As naively expected, as $\lambda $ increases, higher scaled chemical potentials are needed to nucleate the structures. Finally, by increasing the parameter $\alpha $, 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 formation 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 values 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 conditions. 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 particles, 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 nontemplated assembly of curved elastic shells, taking into account all relevant ingredients (i.e. chemical potential, line tension, spontaneous curvature, and elastic contributions) and the potential formation of nonspherical 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 continuous 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 inplane 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 coarsegrained model made using widely different values for the parameters and interaction range confirm that the outcome of the assembly only depends on three scaled parameters: the scaled chemical potential $\mathrm{\Delta}\stackrel{~}{\mu}$, line tension $\lambda $, and FvK number $\gamma $. Thus, the assembly phase diagram is universal, and different protein shells, interaction potentials and coarsegrained models can be recast into a unifying picture of assembly, that could guide the efficient production of artificial viral cages. For instance, our analysis indicates that relatively longrange interactions are desirable to increase the line tension, decrease the FvK number and facilitate the assembly of closed 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 noninfective, cylindrical or beltlike structures. Experimentally, the chemical potential can be tuned by the total protein concentration or by the addition of crowding agents. The line tension (which depends on the strength of the binding interaction), could be modified by the temperature, 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 precise 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 triangulation 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 coassemble 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 selfassembled 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 onestep 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 overestimated, 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 $\gamma =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).
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 diagrams 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 noninfective structure.
Materials and methods
Coarsegrained model and simulation details
Request a detailed protocolThe simulation model, introduced in Aznar and Reguera (2016); Aznar et al. (2018), is coarsegrained at the level of capsomers which are represented as effective spheres of two different diameters: ${\sigma}_{h}$ and ${\sigma}_{p}$, reflecting the fact that hexamers and pentamers are made of a different number of proteins (six and five, respectively). The interaction between capsomers, $V={V}_{LJ}\cdot {V}_{a}\cdot {V}_{tor}$, is modeled using three contributions: a Mielike, an angular, and a torsion potential. The Mielike potential
describes the binding and the excluded volume interaction between a pair of capsomers in terms of their relative distance, $\mathbf{r}}_{ij$ is the equilibrium distance corresponding to the minimum of the potential, r is the distance between capsomers centers, $\u03f5}_{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
where ${\theta}_{ij}$ is the angle between the vector $\mathbf{\Omega}}_{i$, describing the spatial orientation of the capsomer, and the unit vector ${\mathbf{\mathbf{r}}}_{ij}$. The parameter $\nu $ is the preferred angle of interaction between proteins of different capsomers, and the parameter $\alpha $ 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
where ${k}_{t}$ is the torsion constant and $\xi $ is the angle between the planes defined by the unit vector ${\mathbf{\mathbf{r}}}_{ij}$ and both orientation vectors.
The elastic properties of a shell can be related to the main parameters of the interaction. In particular the Young’s modulus is approximately given by $Y=\frac{2nm}{\sqrt{3}}\frac{{\u03f5}_{0}}{{\sigma}^{2}},$ the bending rigidity is $\kappa =\frac{3\sqrt{3}}{8}\frac{{\u03f5}_{0}}{{\alpha}^{2}},$ and the preferred radius of curvature is ${R}_{0}=\frac{\sigma}{2\mathrm{cos}\nu}.$ The line tension of a partially formed cap can be approximated by Luque et al. (2012) $\lambda =\frac{2{\u03f5}_{0}}{\sqrt{3}\sigma},$ and the chemical potential difference that controls the assembly is given by $\mathrm{\Delta}\mu ={k}_{B}T\mathrm{ln}\left({c}_{1}/{c}^{*}\right),$ where ${c}_{1}$ 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 stochastic 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 $\sigma $, their diffusion coefficient $D$, and the binding energy ${\u03f5}_{0}$. In these reduced units, the parameters used in the simulation are: torsion constant ${k}_{t}=1.5$, reduced temperature $T=0.1$, corresponding to a binding energy between capsomers of $10{k}_{B}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 capsomers up to a total of $N=200400$ were initially placed randomly inside a cubic box with periodic boundary conditions. The simulations run for a total of 2 × 10^{9} steps and the final structures were analyzed. To verify the universality of the phase diagram, we performed an extensive set of simulations 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<\nu <1.45$; the bending stiffness in the range $0.05<\alpha <0.4$; and the concentration of capsomers was varied from $\rho =0.005$ to $\rho =0.05$.
Appendix 1
Elastic Model of Shell Assembly
The elastic energy in the formation of a curved shell has two contributions: the bending energy, G_{b}, associated with deviations from the spontaneous curvature, and the inplane elastic energy, G_{s}. The bending energy is described using the expression (Castelnovo, 2017):
where, H is twice the mean curvature of the shell, R_{0} is the spontaneous radius of curvature, K is the Gaussian curvature and κ is the bending modulus.
The expression for the inplane elastic energy of the different analyzed structures is based on previous results from continuum elasticity theory for the deformation of thin plates (Seung and Nelson, 1988; Lidmar et al., 2003; Majidi and Fearing, 2008; Morozov and Bruinsma, 2010; Meng et al., 2014; Grason, 2016; Paquay et al., 2017; Castelnovo, 2017; Schneider and Gompper, 2007). The structures analized are: a hexagonallyordered spherical caps without defects; a spherical cap with one central defect; a spherical cap with n defects; a ribbon; a belt; a cylindrical patch; and a cylinder. The Föpplvon Kárman number (FvK), defined in this work as $\gamma \equiv Y{R}_{0}^{2}/\kappa $, dictates the relative importance of bending and stretching contributions.
Bendingdominated regime
For $\gamma \ll 1$, the bending term dominates and forces all structures to adopt the spontaneous curvature R = R_{0}. We will derive the free energy of formation of the different structures in the small curvature approximation, and compare their relative stability under assembly conditions. The comparison is performed for different shapes having the same area S, that is having the same number of subunits.
The inplane elastic energy of a circular domain of geodesic radius ρ_{0} on a curved spherical surface of radius R is given by Schneider and Gompper (2007); Morozov and Bruinsma, 2010; Meng et al. (2014)
Accordingly, its free energy of formation becomes
or in scaled units
where $\mathrm{\Delta}{g}_{cap}\equiv \frac{\mathrm{\Delta}{G}_{cap}}{4\pi {R}_{0}^{2}Y}$ is the free energy of formation divided by the characteristic elastic energy $4\pi {R}_{0}^{2}Y,\phantom{\rule{thinmathspace}{0ex}}x\phantom{\rule{thinmathspace}{0ex}}\equiv {\rho}_{0}/{R}_{0}$ is the scaled patch radius, $\lambda \equiv \mathrm{\Lambda}/({R}_{0}Y)$ is the scaled line tension, and $\mathrm{\Delta}\stackrel{~}{\mu}\equiv \mathrm{\Delta}\mu /(Y{a}_{1})$ is the scaled chemical potential. The bulk energy grows as x^{2}, the rim energy as x, and the elastic stress as x^{6}. The competition between these three contributions determines the shape of the $\mathrm{\Delta}{g}_{cap}(x)$ curve (see Appendix 1—figure 1). For $\mathrm{\Delta}\stackrel{~}{\mu}$ small, the positive second and third terms of Equation A4 dominate, thus, grows monotonically as shown in Appendix 1—figure 1. However, as $\mathrm{\Delta}\stackrel{~}{\mu}$ increases, there is a particular value
obtained by setting $\frac{d\mathrm{\Delta}{g}_{cap}}{dx}=0$ and $\frac{{d}^{2}\mathrm{\Delta}{g}_{cap}}{d{x}^{2}}=0$, at which an inflection point located at
appears. For $\mathrm{\Delta}\stackrel{~}{\mu}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}\mathrm{\Delta}{\stackrel{~}{\mu}}^{ms}$, the free energy landscape has a maximum, signaling the nucleation barrier, but also a local minimum, x_{min}, corresponding to a locallystable 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 $\mathrm{\Delta}\stackrel{~}{\mu}\phantom{\rule{thinmathspace}{0ex}}\ge \phantom{\rule{thinmathspace}{0ex}}\mathrm{\Delta}{\stackrel{~}{\mu}}^{ms}$ marks the onset of a metastable region where the value of free energy $\mathrm{\Delta}{g}_{cap}({x}_{min})$ of the locallystable shell is larger than its value for the dissembled state at x = 0. A fully stable partial shell is obtained if $\mathrm{\Delta}{g}_{cap}({x}_{min})\le 0$ at the minimum. The onset of the stable region can be obtained from the conditions $\mathrm{\Delta}{g}_{cap}=0$ and $d\mathrm{\Delta}{g}_{cap}/dx=0$, which allows to find the critical value
beyond which $\mathrm{\Delta}{g}_{cap}({x}_{min})\phantom{\rule{thinmathspace}{0ex}}\le \phantom{\rule{thinmathspace}{0ex}}0$. The corresponding value of x is
and represents the minimum size of a partial spherical shell without defects to be stable. The size of stable partial shells grows as $\mathrm{\Delta}\stackrel{~}{\mu}$ increases. For large patches or when $\gamma \ll 1$, the rim energy can be neglected and the location of the stable cap size is described approximately by $x}_{eq}\approx (128\mathrm{\Delta}\stackrel{~}{\mu}{)}^{1/4$.
When the area of the spherical crystalline patch gets large, it becomes favorable to introduce defects to release the elastic stress. The inplane elastic energy of a spherical cap with one 5fold disclination is (Morozov and Bruinsma, 2010; Grason, 2016)
Thus, the free energy of formation, in scaled units, for a spherical cap with one defect at the center becomes
Such mechanism is energetically favorable only if the second term of Equation A10 is negative, that is, if $x\ge \sqrt{2/3}$. The fact that for $x\ge \sqrt{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 spherical cap. Imposing ${x}^{s}\le \sqrt{2/3}$ one gets the requirements $\lambda \le {6}^{7/2}$ and $\mathrm{\Delta}\stackrel{~}{\mu}\le 5/864$.
Using the conditions $\mathrm{\Delta}{g}_{{d}_{1}}=0$ and $d\mathrm{\Delta}{g}_{{d}_{1}}/dx=0$, we can find the onset of stability for a circular cap with one central disclination, described by
with
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$ 5fold disclinations in scaled units is (Grason, 2012; Grason, 2016; Castelnovo, 2017)
where
is the selfenergy of the isolated disclinations, and
is the pairwise interaction of disclinations, with
being ${x}_{\alpha}\equiv {\rho}_{\alpha}/{\rho}_{0}$, ${\rho}_{\alpha}$ the geodesical position of the disclinations and $\left{\mathbf{\mathbf{x}}}_{\alpha}{\mathbf{\mathbf{x}}}_{\beta}\right$ is the (normalized) geodesic distance between disclination $\alpha $ and $\beta $.
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\simeq 2$ in which case, the stability region appears for $\mathrm{\Delta}\stackrel{~}{\mu}$ larger than
and the cap consists of a fully closed shell with defects.
The free energy of formation of a ribbon of length $L$ and width $W$, growing on the surface of a sphere of radius ${R}_{0}$ is (Schneider and Gompper, 2007; Majidi and Fearing, 2008)
The energetic advantage of this configuration over the spherical cap is that, for a fixed width $W$, the inplane 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=\pi {\rho}_{0}^{2}=LW$. Thus, this energy can be rewritten in dimensionless terms as
where $w\equiv W/{R}_{0}$. The optimal width of the ribbon is obtained by minimization that is, from $\partial \mathrm{\Delta}{g}_{rib}/\partial w=0$, yielding
As the area of the patch increases, the ribbon grows longitudinally at a nearly fixed optimal width up to the point where $l=2\pi $, where it forms a closed belt with energy
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 ribbonlike 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
and the stability region for the belt occurs when $\mathrm{\Delta}\stackrel{~}{\mu}$ is larger than
Since $\mathrm{\Delta}{\stackrel{~}{\mu}}_{belt}^{s}<\mathrm{\Delta}{\stackrel{~}{\mu}}_{cap}^{s}$, that leads to the important result that closed belts become stable before defectless spherical caps. Therefore, defectless spherical caps can only be at most metastable structures.
Finally, an alternative to the curved ribbon is a cylindrical stripe of scaled width $w$ and length $l$ having one of his principal radius of curvature zero and a scaled energy
As in the case of the ribbon and the belt, the cylindrical patch eventually closes when $l=2\pi r$ into a cylinder whose energy of formation is
The cylinder has the advantage of not having any inplane stretching cost, but it has a bending energy penalty described by the second term inside the parentesis. An energetically favorable cylinder requires $\mathrm{\Delta}\stackrel{~}{\mu}>1/(2\gamma )$ for its formation, meaning that cylindrical shells cannot be formed in the bendingdominated limit, corresponding to $\gamma \to 0$.
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 ${R}_{0}$. Particularizing Equation A1, the bending energy of a sphere of radius $R$ in scaled units is
while for a cylinder of radius $R$ reads
where $r\equiv R/{R}_{0}$ and the Föpplvon Kárman number (FvK), quantifying the ratio of stretching and bending energies, is still defined as $\gamma \equiv Y{R}_{0}^{2}/\kappa $.
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
In the general case, the optimal radius of the shell is obtained from the condition $\partial \mathrm{\Delta}{g}_{cap}/\partial r=0$, yielding
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
In the limit of large shells or Fvk numbers, the radius goes as
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\simeq 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 5fold disclinations. The free enery of formation of a spherical shell with one central defect is
which becomes favorable over the defectless case when $x/r\ge \sqrt{2/3}$. The formation energy of a spherical shell with $n$defects is
where
and ${g}_{{s}_{2}}$ is given by Equation A15. The free energy of a ribbon of scaled width $w$ and length $l$ in the general case becomes
The optimal width and radius of curvature of the ribbon, $w$ and $r$, are obtained by minimization that is, from $\partial \mathrm{\Delta}{g}_{rib}/\partial w=0$ and $\partial \mathrm{\Delta}{g}_{rib}/\partial r=0$, yielding
and
Finally, the free energies of a belt and a cylinder are
and
respectively.
Data availability
All data generated or analysed during this study are included in the manuscript and supporting files.
References

Neutral versus charged defect patterns in curved crystalsPhysical Review E 94:1–14.https://doi.org/10.1103/PhysRevE.94.013003

Viral nanomechanics with a virtual atomic force microscopeJournal of Physics: Condensed Matter 30:264001.https://doi.org/10.1088/1361648X/aac57a

Physical ingredients controlling stability and structural selection of empty viral capsidsThe Journal of Physical Chemistry B 120:6147–6159.https://doi.org/10.1021/acs.jpcb.6b02150

Physical principles in the construction of regular virusesCold Spring Harbor Symposia on Quantitative Biology 27:1–24.https://doi.org/10.1101/SQB.1962.027.001.005

Viral selfassembly pathway and mechanical stress relaxationPhysical Review E 95:052405.https://doi.org/10.1103/PhysRevE.95.052405

BookPrinciples of Virology: Molecular Biology, Pathogenesis, and Control of Animal Viruses (Second Edition)Washington, DC: ASM Press.

Clathrin assembly regulated by adaptor proteins in CoarseGrained modelsBiophysical Journal 111:222–235.https://doi.org/10.1016/j.bpj.2016.06.003

Perspective: geometrically frustrated assembliesThe Journal of Chemical Physics 145:110901.https://doi.org/10.1063/1.4962629

Dynamic Pathways for Viral Capsid AssemblyBiophysical Journal 91:42–54.https://doi.org/10.1529/biophysj.105.076851

Recent advances in coarsegrained modeling of virus assemblyCurrent Opinion in Virology 18:36–43.https://doi.org/10.1016/j.coviro.2016.02.012

Virus particle maturation: insights into elegantly programmed nanomachinesCurrent Opinion in Structural Biology 20:210–216.https://doi.org/10.1016/j.sbi.2010.01.004

Modelling the selfassembly of virus capsidsJournal of Physics: Condensed Matter 22:104101.https://doi.org/10.1088/09538984/22/10/104101

Simian virus 40 VP1 capsid protein forms polymorphic assemblies in vitroJournal of General Virology 84:1899–1905.https://doi.org/10.1099/vir.0.190670

Stress induced branching of growing crystals on curved surfacesPhysical Review Letters 116:135502.https://doi.org/10.1103/PhysRevLett.116.135502

Budding of crystalline domains in fluid membranesPhysical Review E 68:1–15.https://doi.org/10.1103/PhysRevE.68.061905

Virus shapes and buckling transitions in spherical shellsPhysical Review E 68:051910.https://doi.org/10.1103/PhysRevE.68.051910

Physics of shell assembly: line tension, hole implosion, and closure catastropheThe Journal of Chemical Physics 136:184507.https://doi.org/10.1063/1.4712304

Adhesion of an elastic plate to a sphereProceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 464:1309–1317.https://doi.org/10.1098/rspa.2007.0341

Spontaneouscurvature theory of clathrincoated membranesBiophysical Journal 74:2862–2875.https://doi.org/10.1016/S00063495(98)779937

Invariant polymorphism in virus capsid assemblyJournal of the American Chemical Society 131:2606–2614.https://doi.org/10.1021/ja807730x

The equilibrium structure of selfassembled protein nanocagesNanoscale 10:22802–22809.https://doi.org/10.1039/C8NR07202G

Mechanisms of virus assemblyAnnual Review of Physical Chemistry 66:217–239.https://doi.org/10.1146/annurevphyschem040214121637

Major capsid reinforcement by a minor protein in herpesviruses and phageNucleic Acids Research 42:9096–9107.https://doi.org/10.1093/nar/gku634

ThesisStructure, Mechanical Properties, and SelfAssembly of Viral Capsids PhD thesisUniversitat de Barcelona.

Shapes of crystalline domains on spherical fluid vesiclesEurophysics Letters 70:136–142.https://doi.org/10.1209/epl/i2004104642

Defects in flexible membranes with crystalline orderPhysical Review A 38:1005–1018.https://doi.org/10.1103/PhysRevA.38.1005

Design of virusbased nanomaterials for medicine, biotechnology, and energyChemical Society Reviews 45:4074–4126.https://doi.org/10.1039/C5CS00287G

Origin of icosahedral symmetry in virusesPNAS 101:15556–15560.https://doi.org/10.1073/pnas.0405844101

Classical nucleation theory of virus capsidsBiophysical Journal 90:1939–1948.https://doi.org/10.1529/biophysj.105.072975
Decision letter

Pierre SensReviewing Editor; Institut Curie, PSL Research University, CNRS, France

Aleksandra M WalczakSenior Editor; École Normale Supérieure, France
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
This study reports a theoretical analysis of the mechanical factors influencing the selfassembly (and misassembly) of protein subunits into viral capsids. This subject has been studied for a long time, but the present paper offers new insight regarding the importance of contribution from the protein line tension and preferred curvature. The main result of the study is (i) a universal phase diagram, which shows that viruses with high mechanical resistance cannot be selfassembled directly as spherical structures, and (ii) the ensuing suggestion of the need of a maturation step to stiffen the capsid after assembly. Inducing misassembly could be a promising route to hinder viral infections.
Decision letter after peer review:
Thank you for submitting your article "Triggering misassembly in viral capsid formation by elastic frustration" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Aleksandra Walczak as the Senior Editor. The reviewers have opted to remain anonymous.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
This paper combines a theoretical and a numerical study of the dependence of the selfassembly or misassembly process of viral capsids on parameters such as the chemical potential difference of protein subunits in and out of the capsid, the line tension, and the stretching and bending energies of the planar assembly, characterized by the Fopplvon Karman number. The main result of the study is (i) an universal phase diagram, which shows that viruses with high mechanical resistance cannot be selfassembled directly as spherical structures, and (ii) which justifies the need of a maturation step to stiffen the capsid after assemble. The authors suggest that this might be a promising route to hinder viral infections by inducing misassembly.
The reviewers found your study of the selfassembly (and misassembly) of protein subunits into capsids to be interesting and potentially important. Although the subject of viral capsid assembly has been studied for a long time, investigating the effect of the line tension, chemical potential and the preferred curvature of the protein assembly is a welcome addition to the existing literature.
The reviewers nevertheless have quite a few comments and questions, some of which are major points that should be addressed before a final decision can be made on the suitability of this submission of publication at eLife.
Essential revisions:
The reviewers nevertheless have some questions, some of which are major points that should be addressed before a final decision can be made on the suitability of this submission of publication at eLife.
1) The title of the paper focusses on misassembly, while most of the text discusses focus on assembly rather than misassembly. You should consider modifying the title to account for this.
2) Several important references are not cited. It is important to compare and contrast the results of the current study in some detail with earlier investigations. These include:
Defect in curved crystal:
– Azadi and Grason Phys. Rev. E 2016
Mechanisms of virus assembly:
– M.F. Hagan and D. Chandler, Biophys. J. 91, 42 (2006);
– O.M. Elrad, M.F. Hagan, Nano letters 8, 3850 (2008);
– H. Nguyen et al., J. Am. Chem. Soc. 131, 2606 (2009);
– I.G. Johnston et al., J. Phys. Condens. Matter 22, 104101 (2010);
– J.D. Perlmutter, M.F. Hagan, Ann. Rev. Phys. Chem. 66, 217 (2015).
Formation of clathrincoated pits:
– R.J. Mashl et al., Biophys. J. 74, 2862 (1998);
– T. Kohyama et al., Phys. Rev. E 68, 061905 (2003);
– M. Giani et al., Biophys. J. 111, 222 (2016).
3) A table of numerical values is provided (Table 1). There are some questions about how the Young's modulus E is evaluated. In some references provided, in some of the references, E is not mentioned directly so the authors might have estimated it indirectly. Could the procedure for this evaluation be explained?
4) In order to determine the lowest energy sate, shape with similar number of subunits should be compared. Here, a rescaled size x is used, the definition of which appear to depend upon the shape considered. Hence a given x does not correspond to the same number of subunit for different shape. This should be corrected, and energies of shape with the same number of subunit should be compared.
5) The analytical results are based on a thermodynamic approach where the chemical potential is fixed, which implies a reservoir of particle. It is not clear how this compares with the numerical simulations. Is the chemical potential fixed during the assembly process in the simulations? This point should be made clear, and if the two approaches are based on different thermodynamic assumptions, how meaningfully can they be compared?
6) How does one relate the number of capsomers on the spherical surface, which is controlled by the chemical potential difference, and the number of particles. In the simulation, do we access the number of particles or the number of capsomers?
7) Regarding the bendingdominated regime, it seems a bit strange to rescale the chemical potential and line tension with the Youngs modulus, since this regime should contain the case where the Fopplvon Karman number (and the twodimensional Young modulus Y) vanish, which implies that tilde(Δµ) and λ diverge. What is the optimal structure in this regime, where one does not expect differences been defectfree caps and caps with defects? This regime is not apparent in Figure 3A.
8) For the spherical cap with defect, how is the spatial distribution of the defect should be discussed? Does the energy used to compare with other shapes correspond to the lowest energy when defect localisation is optimised? This should be discussed in some details, since defect distribution is clearly an important factor influencing the energy of the cap.
9) The authors state that spherical capsids cannot selfassemble at large FvK. This statement is probably correct for strictly spherical capsids. However, aggregates could form as planar patches, then bend and form disclinations at sufficient size and line tension – see references on formation of clathrincoated pits. Such structures are not accounted for in the analytical theory. Do they show up in the simulations, and if not, why not? Here again, the distribution of defect is a crucial factor that should be discussed.
10) Figure 3: Similar phase diagrams have been calculated by Schneider and Gompper. In their calculation, also structures consisting of several caps or several belts have been predicted. Are these structures relevant here as well? This should be discussed.
11) In the high concentration regimes "… many fragments of spherical capsids that cannot grow any further.…". Why can't they merge, rearrange their internal structure, and thereby reduce linetension energy?
12) In the general case of non zero FvK numbers, the spontaneous radius of curvature plays an important role in the model, and this is discussed thoroughly by the authors. However, in the phase diagram, it seems that its discussion has been dismissed, or maybe its value has been fixed? In any case, some clarification on how the phase diagram depend on the value of spontaneous curvature is required. For example, it has been shown the reference by Castelnovo in 2017 that for small enough spontaneous curvature, cylinders should dominate. It is not clear how it compares to the many cases discussed here.
13) It is not clear why cylinders should dominate for large FvK numbers and large rescaled spontaneous curvature. Could the author elaborate on this point?
14) The authors should make more effort to recast their results in terms of experimentally useful quantities, so that it becomes apparent how to experimentally explore the phase diagram.
It is proposed that chemical or physical factors that increase the FvK number or reduce the line tension or the effective concentration are potential targets to prevent viral replication. It would be quite insightful to provide a more quantitative and practical version of this statement. Regarding the chemical potential, this statement could be turned into a capsid protein concentration at which viral assembly is expected to fail, given particular values of the other parameters. How does this concentration relates to expected concentration in cells? Regarding the mechanical parameter, in which range can they be expected to vary under the action of which factor, and would this be enough in practice to prevent capsid formation? You should provide more precise statement regarding how misassembly could be induced to hinder viral infections.
https://doi.org/10.7554/eLife.52525.sa1Author response
Essential revisions:
The reviewers nevertheless have some questions, some of which are major points that should be addressed before a final decision can be made on the suitability of this submission of publication at eLife.
1) The title of the paper focusses on misassembly, while most of the text discusses focus on assembly rather than misassembly. You should consider modifying the title to account for this.
We have changed the title of the work now entitled: “Shape selection and misassembly in viral capsid formation by elastic frustration”.
2) Several important references are not cited. It is important to compare and contrast the results of the current study in some detail with earlier investigations. These include:
Defect in curved crystal:
– Azadi and Grason Phys. Rev. E 2016
Mechanisms of virus assembly:
– M.F. Hagan and D. Chandler, Biophys. J. 91, 42 (2006);
– O.M. Elrad, M.F. Hagan, Nano letters 8, 3850 (2008);
– H. Nguyen et al., J. Am. Chem. Soc. 131, 2606 (2009)R;
– I.G. Johnston et al., J. Phys. Condens. Matter 22, 104101 (2010);
– J.D. Perlmutter, M.F. Hagan, Ann. Rev. Phys. Chem. 66, 217 (2015).
Formation of clathrincoated pits:
– R.J. Mashl et al., Biophys. J. 74, 2862 (1998);
– T. Kohyama et al., Phys. Rev. E 68, 061905 (2003);
– M. Giani et al., Biophys. J. 111, 222 (2016).
We thank the reviewers for pointing out those relevant references. We have included all of them in the revised manuscript. More specifically, we have added a sentence at the beginning of the third paragraph to reference previous works on the mechanisms of virus assembly, including an additional reference (Hagan and Zandi, 2016). We have added another sentence to refer to the formation of clathrincoated pits as well as carboxysomes and colloidosomes. Finally, we have included the reference to the work by Azadi and Hagan in the Discussion of crystal growth on spherical templates on the fourth paragraph.
3) A table of numerical values is provided (Table 1). There are some questions about how the Young's modulus E is evaluated. In some references provided, in some of the references, E is not mentioned directly so the authors might have estimated it indirectly. Could the procedure for this evaluation be explained?
The values of Young’s modulus reported on Table 1 are all, except that of SV40, provided by the references included in the caption. More specifically, the Young’s modulus of CCMV and λ capsid are listed in Table 1 of Mateu, 2012 and were evaluated using Finite Elements simulations by Michel et al., 2006 and Ivanoska et al., 2007, respectively. The Young’s modulus of λ procapsid was evaluated and reported by SaeUeng et al., 2014 from the spring constant measured by AFM using the standard thin shell formula. Finally, the Young’s modulus of SV40 was evaluated from the spring constant measured by AFM by van Rosmalen et al., 2018, using the standard thin shell formula k = 2.25 E h^2/R. We have included a sentence in the caption of Table 1 to clarify that. We have also added a reference to Ivanoska et al., 2004 and 2007 in the caption of Table 1 for the sake of completeness.
4) In order to determine the lowest energy sate, shape with similar number of subunits should be compared. Here, a rescaled size x is used, the definition of which appear to depend upon the shape considered. Hence a given x does not correspond to the same number of subunit for different shape. This should be corrected, and energies of shape with the same number of subunit should be compared.
Although, for simplicity and consistency, we write the equations and represent the energies for different shapes as a function of x, we do in fact determine the lowest energy state by comparing the energies for different shapes having the same area π x^{2}, which is equivalent to having the same number of subunits. The variable x=ρ_{0}/R_{0} is always defined as the radius of a circular patch of subunits, ρ_{0}, divided by the spontaneous radius R_{0}, so that its normalized area is π ρ_{0}^{2}/(4π R_{0}^{2}) = x^{2}/4. For other shapes, in order to compare structures with the same number of subunits, what we do is to impose that they have the same area. Indeed, we use this requirement of equal areas to rewrite the scaled free energies of all shapes in terms of x. For example, for a ribbon of length L and width W, the normalized area is A_{rib}/(4π R_{0}^{2})=LW/(4π R_{0}^{2})=lw/(4π). Thus, a ribbon with the same area (and thus, with the same number of subunits) as a circular patch implies lw=π x^{2}, as written before Equation 5.
To clarify this point, we have added a sentence in the manuscript before Equation 2 and we have included this Discussion in the first paragraph of the bending dominated section and before Equation 19 in the Appendix.
5) The analytical results are based on a thermodynamic approach where the chemical potential is fixed, which implies a reservoir of particle. It is not clear how this compares with the numerical simulations. Is the chemical potential fixed during the assembly process in the simulations? This point should be made clear, and if the two approaches are based on different thermodynamic assumptions, how meaningfully can they be compared?
This is a very important point. While theory assumes a controlled chemical potential (equivalent to having a reservoir of particles), the simulations are done at fixed number of particles. This implies that in the simulation as the assembly proceeds, the concentration of the remaining free particles, and consequently the chemical potential, decrease. This may cause discrepancies when comparing both approaches, specifically if a precise quantitative comparison is intended. For this reason, we have not intended to reproduce with precision the values of the chemical potential or the borders of the phase diagram using the simulations. Nonetheless, for simulations with a large number of particles, the change in concentration and chemical potential due to the formation of a shell would be relatively small and can be neglected. But even if the chemical potential is not strictly constant in our simulations, the semiquantitative comparison we were interested in is not altered since in any case we can compare structures with different relative concentrations and we can still say which structure will appear at the larger concentration, for example.
This important discussion and clarification about finite number of particles’ effects is now included at the end of section Simulation.
6) How does one relate the number of capsomers on the spherical surface, which is controlled by the chemical potential difference, and the number of particles. In the simulation, do we access the number of particles or the number of capsomers?
In the simulations we represent each subunit, i.e. (hexameric) capsomer, by a coarsegrained spherical particle. The simulations are performed with a fixed total number of subunits N. The initial configuration consists of a small spherical cap of 19 subunits with the spontaneous radius R_{0} as initial seed, and the remaining N19 subunits are randomly distributed inside a fixedsized simulation box. Thus, the total concentration of subunits (or capsomers) is fixed in a given simulation. Please note that in the simulations we are only inserting an initial seed, but we are not forcing the capsomers to be in any kind of spherical or cylindrical template. Instead, the subunits spontaneously selfassemble into the energeticallypreferred shape and size. Thus, it is the total number of subunits (i.e. capsomers) that is controlled, but it is the time evolution of the seeded system, without any bias or template, which determine how many subunits end up forming the energetically preferred structure.
7) Regarding the bendingdominated regime, it seems a bit strange to rescale the chemical potential and line tension with the Youngs modulus, since this regime should contain the case where the Fopplvon Karman number (and the twodimensional Young modulus Y) vanish, which implies that tilde(Δµ) and λ diverge. What is the optimal structure in this regime, where one does not expect differences been defectfree caps and caps with defects. This regime is not apparent in Figure 3A.
The situation of a zero Young’s modulus is unrealistic, since it would correspond to a mechanical unstable material that could be stretched or compressed without any energy cost. Thus, we are always considering a finite Y to normalize the chemical potential and the line tension. Rather than with a negligible Young modulus, the bendingdominated regime, corresponding to the limit in which the Fopplvon Karman number vanish, is achieved by having a very large bending modulus compared to the Young modulus times R_{0}^{2}. In this case the structure cannot be bend out from its preferred curvature R_{0} since it would cost a lot of bending energy. This is the case depicted in Figure 3A. On the other hand, the limit in which Y vanishes would imply a structure that does not pay any stretching energy. In this case, there is almost no difference for the structure to incorporate defects or not, but also only the first two terms of Equation 1 would remain. Thus, at equal areas, partially formed caps with the preferred curvature would prevail over the ribbon due to the smaller edge energy. In addition, in the limit of vanishingly small Y, the minimum due to elastic frustration will only occur at very large patch areas, exceeding the values corresponding to a closed shell. Euler’s theorem states that at least 12 defects are required to have a closed shell. So, in the bending dominated regime and in the limit Y very small, corresponding to very large values of the scaled ∆µ and λ, the optimal structure would still be a closed spherical cap with defects. In fact, we have shown in the Appendix that the defectless spherical cap is always metastable. Therefore, all realistic situations are already accounted for in the bendingdominated scenario depicted in Figure 3A.
8) For the spherical cap with defect, how is the spatial distribution of the defect should be discussed? Does the energy used to compare with other shapes correspond to the lowest energy when defect localisation is optimised? This should be discussed in some details, since defect distribution is clearly an important factor influencing the energy of the cap.
In Li et al., 2018, the energy of an incomplete cap with one defect placed at an arbitrary location is calculated. It is found that the Gaussian curvature attracts the disclination to the center of the cap while the defect selfenergy pushes it towards the boundary. The net result is that the minimum energy corresponds to the defect located off the center of the cap. In our calculation we have used the approximation in which the defect is located at the center. However, we have numerically verified that this approximation introduces only a very small error, of less than 0.25% , in our calculations for the scaled energy. This means that not noticeable effect is observed when the exact expression with the offcenter defect is considered. On the other hand, for a cap with more than one disclination, the repulsive interaction among them, as well as their repulsion with the free boundary, contribute to distribute them regularly, as considered in our calculation. When the cap finally closes, the maximal number of disclinations is 12, their regular distribution locates them along the vertices of an icosahedron, which corresponds to the lowest energy configuration.
9) The authors state that spherical capsids cannot selfassemble at large FvK. This statement is probably correct for strictly spherical capsids. However, aggregates could form as planar patches, then bend and form disclinations at sufficient size and line tension – see references on formation of clathrincoated pits. Such structures are not accounted for in the analytical theory. Do they show up in the simulations, and if not, why not? Here again, the distribution of defect is a crucial factor that should be discussed.
In our model, the condition for planar patches to assemble is ∆µ(tilde)>=1/γ, but only if the preferred curvature is zero, otherwise, the assembly of cylinders will be preferable due to a smaller bending energy. If the preferred curvature is zero, there would be no reason to bend out of the plane, which would also cost a lot of stretching energy, unless an additional external mechanism forcing the bending is present (which is the case in clathrincoated pits). Thus, in our case, at large FvK the theory does not predict spherical capsids as preferred structure and we indeed have not found them in the simulations performed at large FvK.
10) Figure 3: Similar phase diagrams have been calculated by Schneider and Gompper. In their calculation, also structures consisting of several caps or several belts have been predicted. Are these structures relevant here as well? This should be discussed.
In Schneider and Gompper’s study, the crystalline domains appear on a spherical fluid vesicle acting as a template. Their phase diagrams are constructed in terms of the relative area covered by the crystalline phase, considering the possibility that it may consist of independent belts or caps growing on the same vesicle. In our study, on the other hand, there is no preexisting vesicle or template and the assembled caps or belts form independently of each other depending on the chemical potential and not by imposing a fixed target area. What it is possible, however, is that during the assembly process, shapes similar to a central cap with ribbonlike protrusions (similar to a flower) or branched ribbonlike structures (Appendix 1—figure 3E) may appear. These more complex shapes have not been considered in our study.
11) In the high concentration regimes "… many fragments of spherical capsids that cannot grow any further.…". Why can't they merge, rearrange their internal structure, and thereby reduce linetension energy?
The kinetic trapping mechanism we discuss in that part of the text, is due to the formation of many fragments, that in a simulation with a fixed number of subunits causes that the concentration of the remaining free monomers eventually reduces to a value that is smaller than the critical concentration needed to continue the assembly of the fragments (this is related to point 5 above). Of course, these fragments would be subjected to Brownian motion and may eventually meet at the right orientations and merge, reducing the line tension. But it is very unlikely that many fragments containing a different number of subunits and defects end up forming a perfect closed shell. (It would be similar to the formation of a nice single crystal by the coarsening of multiple small crystalline domains grown independently.) In any case, the kinetic trapping mechanism and the possibility of forming the right structure at longer times by coarsening, go beyond of the present work and will be discussed in a future publication, as mentioned in the manuscript.
12) In the general case of non zero FvK numbers, the spontaneous radius of curvature plays an important role in the model, and this is discussed thoroughly by the authors. However, in the phase diagram, it seems that its discussion has been dismissed, or maybe its value has been fixed? In any case, some clarification on how the phase diagram depend on the value of spontaneous curvature is required. For example, it has been shown the reference by Castelnovo in 2017 that for small enough spontaneous curvature, cylinders should dominate. It is not clear how it compares to the many cases discussed here.
The phase diagrams (Figure 3) have been drawn using scaled variables. The scaled line tension is λ=Λ/(R_{0}Y), where R_{0} is the spontaneous radius of curvature. We can make contact with Castelnovo results noticing that small spontaneous curvature means a large R_{0} which for a given line tension Λ implies a small value of λ which is the region where belts tend to appear. On the other hand, the phase diagram show that cylinders appear at Δμ(tilde)>=1/(2γ), which, after using the definition of the scaled quantities can be rewritten as Δμ>=κ a_{1}/(2R_{0}^{2}). In other words, cylinders appear more easily (smaller Δμ required) for larger values of R_{0}, in complete agreement with Castelnovo results. This is now clarified and discussed in the manuscript.
13) It is not clear why cylinders should dominate for large FvK numbers and large rescaled spontaneous curvature. Could the author elaborate on this point?
At large FvK numbers, corresponding to the regime where stretching dominates over bending, cylinders 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 cylindrical structure having a radius equal to the spontaneous radius R_{0}, i.e. r=1, will minimize the bending penalty and will have a free energy of formation, according to Equation 12, that decreases unboundedly with the size when Δμ(tilde)>=1/(2γ). 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 Δμ(tilde)>=1/(2γ). The larger the γ (FvK), the smaller the Δμ(tilde) 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).
This discussion has been incorporated in the manuscript, just before the simulations section.
14) The authors should make more effort to recast their results in terms of experimentally useful quantities, so that it becomes apparent how to experimentally explore the phase diagram.
It is proposed that chemical or physical factors that increase the FvK number or reduce the line tension or the effective concentration are potential targets to prevent viral replication. It would be quite insightful to provide a more quantitative and practical version of this statement. Regarding the chemical potential, this statement could be turned into a capsid protein concentration at which viral assembly is expected to fail, given particular values of the other parameters. How does this concentration relates to expected concentration in cells? Regarding the mechanical parameter, in which range can they be expected to vary under the action of which factor, and would this be enough in practice to prevent capsid formation? You should provide more precise statement regarding how misassembly could be induced to hinder viral infections.
The main goal of the work was to provide a general understanding of the assembly phase diagram that could guide future experiments. But in order to make precise quantitative predictions that could be experimentally tested, we need a few parameters that, as far as we know, have unfortunately not been measured for real viruses, yet. For instance, the critical concentration c* is needed in order to make a quantitative prediction of the capsid protein concentration at which viral assembly is expected to fall. Without accurate evaluation of these parameters it is too adventurous to make reliable quantitative predictions. We hope that our work will stimulate future experimental work aimed at measuring these important parameters for different viruses of interests. In any case, we have added a paragraph in the Discussion to make a closer connection between the physical parameters controlling the assembly and experiments.
The added paragraph reads:
“Experimentally, the chemical potential can be tuned by the total protein concentration or by the addition of crowding agents. The line tension (which depends on the strength of the binding interaction), could be modified by the temperature, 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 precise quantitative connection between the physical parameters controlling the assembly and experiments.”
https://doi.org/10.7554/eLife.52525.sa2Article and author information
Author details
Funding
Universidad Nacional Autónoma de México (DGAPA IN110516)
 Carlos I Mendoza
Universidad Nacional Autónoma de México (DGAPA IN103419)
 Carlos I Mendoza
Gobierno de Espana (FIS201567837P)
 David Reguera
Ministerio de Economía y Competitividad (PGC2018098373BI00)
 David Reguera
European Regional Development Fund (PGC2018098373BI00)
 David Reguera
Universidad Nacional Autónoma de México (Sabbatical Fellowship)
 Carlos I Mendoza
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
CIM appreciates the hospitality of Dr. David Reguera and the Departament de Física de la Matèria Condensada of the Universitat de Barcelona where this work was carried out during a sabbatical leave and subsequent visits. CIM received financial support provided by DGAPAUNAM through a Sabbatical Fellowship and by grants DGAPA IN110516 and IN103419. DR acknowledges funding from the Spanish government through grants FIS201567837P and PGC2018098373BI00 (MINECO/FEDER, UE).
Senior Editor
 Aleksandra M Walczak, École Normale Supérieure, France
Reviewing Editor
 Pierre Sens, Institut Curie, PSL Research University, CNRS, France
Publication history
 Received: October 7, 2019
 Accepted: April 6, 2020
 Accepted Manuscript published: April 21, 2020 (version 1)
 Version of Record published: April 24, 2020 (version 2)
Copyright
© 2020, Mendoza and Reguera
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,509
 Page views

 192
 Downloads

 10
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Cell Biology
 Physics of Living Systems
Motile cilia are hairlike cell extensions that beat periodically to generate fluid flow along various epithelial tissues within the body. In dense multiciliated carpets, cilia were shown to exhibit a remarkable coordination of their beat in the form of traveling metachronal waves, a phenomenon which supposedly enhances fluid transport. Yet, how cilia coordinate their regular beat in multiciliated epithelia to move fluids remains insufficiently understood, particularly due to lack of rigorous quantification. We combine experiments, novel analysis tools, and theory to address this knowledge gap. To investigate collective dynamics of cilia, we studied zebrafish multiciliated epithelia in the nose and the brain. We focused mainly on the zebrafish nose, due to its conserved properties with other ciliated tissues and its superior accessibility for noninvasive imaging. We revealed that cilia are synchronized only locally and that the size of local synchronization domains increases with the viscosity of the surrounding medium. Even though synchronization is local only, we observed global patterns of traveling metachronal waves across the zebrafish multiciliated epithelium. Intriguingly, these global wave direction patterns are conserved across individual fish, but different for left and right nose, unveiling a chiral asymmetry of metachronal coordination. To understand the implications of synchronization for fluid pumping, we used a computational model of a regular array of cilia. We found that local metachronal synchronization prevents steric collisions, cilia colliding with each other, and improves fluid pumping in dense cilia carpets, but hardly affects the direction of fluid flow. In conclusion, we show that local synchronization together with tissuescale cilia alignment coincide and generate metachronal wave patterns in multiciliated epithelia, which enhance their physiological function of fluid pumping.

 Immunology and Inflammation
 Physics of Living Systems
One of the feats of adaptive immunity is its ability to recognize foreign pathogens while sparing the self. During maturation in the thymus, T cells are selected through the binding properties of their antigenspecific Tcell receptor (TCR), through the elimination of both weakly (positive selection) and strongly (negative selection) selfreactive receptors. However, the impact of thymic selection on the TCR repertoire is poorly understood. Here, we use transgenic Nur77mice expressing a Tcell activation reporter to study the repertoires of thymic T cells at various stages of their development, including cells that do not pass selection. We combine highthroughput repertoire sequencing with statistical inference techniques to characterize the selection of the TCR in these distinct subsets. We find small but significant differences in the TCR repertoire parameters between the maturation stages, which recapitulate known differentiation pathways leading to the CD4+ and CD8+ subtypes. These differences can be simulated by simple models of selection acting linearly on the sequence features. We find no evidence of specific sequences or sequence motifs or features that are suppressed by negative selection. These results favour a collective or statistical model for Tcell self nonself discrimination, where negative selection biases the repertoire away from self recognition, rather than ensuring lack of selfreactivity at the singlecell level.