Virus Capsid Dissolution Studied by Microsecond Molecular Dynamics Simulations

Dissolution of many plant viruses is thought to start with swelling of the capsid caused by calcium removal following infection, but no high-resolution structures of swollen capsids exist. Here we have used microsecond all-atom molecular simulations to describe the dynamics of the capsid of satellite tobacco necrosis virus with and without the 92 structural calcium ions. The capsid expanded 2.5% upon removal of the calcium, in good agreement with experimental estimates. The water permeability of the native capsid was similar to that of a phospholipid membrane, but the permeability increased 10-fold after removing the calcium, predominantly between the 2-fold and 3-fold related subunits. The two calcium binding sites close to the icosahedral 3-fold symmetry axis were pivotal in the expansion and capsid-opening process, while the binding site on the 5-fold axis changed little structurally. These findings suggest that the dissociation of the capsid is initiated at the 3-fold axis.


Introduction
Non-enveloped icosahedral viruses often contain binding sites for divalent cations, usually Ca 2z . The ions are typically bound between coat proteins or on the icosahedral symmetry axes. This is broadly observed in three plant virus taxa: the family Tombusviridae (and an associate satellite virus), the genus Sobemoviruses and the family Bromoviridae [1][2][3][4]. Binding sites for calcium ions have also been found in bacteriophages of the Leviviridae family [5], fish and insect viruses of the Nodaviridae family [6] and in the Picornaviridae family, e.g. several human rhinoviruses [7].
In many of the plant viruses it is possible to induce a conformational change in vitro by removing the ions, either by a chelating agent such as ethylenediaminetetraacetic acid (EDTA) or by exhaustive dialysis against deionized water. Ion-deprived virions reversibly expand on the order of 5-10% at neutral or slightly alkaline pH. In the swollen state internal parts of the virion as well as the RNA molecule may become susceptible to degrading enzymes [8,9]. Chelation of the metal ions is also required for synthesis of virus proteins in cell-free translation systems [9]. Only two low-resolution crystal structures of expanded virons are available: tobacco bushy stunt virus (TBSV) at 8 Å [10] and satellite tobacco necrosis virus (STNV) at 7.5 Å [11]. The radial increases are about 11% and 4%, respectively. In addition, an expanded cowpea chlorotic mottle virus (CCMV) virion was imaged with cryo-electron microscopy at 29 Å and interpreted using rigid body fitting of the high-resolution structures of the native proteins [4]. The dynamic nature of the swelling process as well as the limited resolution of swollen virus particles structures prompted us to perform a simulation study of the capsid of STNV, with and without bound Ca 2z , over one microsecond. The simulations allowed us to reproduce the swelling behavior upon removal of the calcium in silico and develop an atomistic description of the process.
The T = 1 capsid of STNV consists of 60 identical coat proteins with one protein per icosahedral asymmetric unit. The coat protein is 195 amino acid residues long where residues  make up the main domain that constitutes the capsid shell. The virions readily crystallize and the major part of the coat protein has been resolved by X-ray crystallography [2,[12][13][14]. The shell domain at the C-terminus folds as a b-jelly roll similar to many other single-stranded RNA plant viruses. Residues 12-24 form a helical structure that together with the helices of two neighboring subunits form a short stalk that projects inwards into the central cavity around the icosahedral 3-fold axis. The first 11 residues at the N-terminus are disordered and cannot be detected in the electron density maps -in the simulations these residues were modeled as a helix as well. This N-terminal arm and the interior surface of the capsid are lined with positively charged residues that presumably interact with the single-stranded positive-sense RNA molecule [14]. The 1239 nucleotide long genome encompasses only one open reading frame that encodes the coat protein and hence STNV is dependent on the co-infection of a helper virus (tobacco necrosis virus) for copying its RNA genome.
The capsid has three different types of Ca 2z binding sites ( Figure 1). Type I is between two subunits close to the 3-fold symmetry axis. The protein ligands are the carboxyl groups of Asp194 and Glu25 as well as the main chain carbonyl oxygens of Ser61 and Gln64. Type II is on the 3-fold symmetry axis 8.05 nm from the center of the virion. It is coordinated by the carboxyl groups of three Asp55 residues. Type III is on the 5-fold symmetry axis 9.04 nm from the center. This Ca 2z is coordinated by the main chain carbonyl oxygen of five Thr138 residues. In total the capsid can accommodate 92 Ca 2z ions (60 at type I sites, 20 at type II sites and 12 at type III sites).
Simulations were performed of the capsid with and without Ca 2z at two different salt concentrations for one microsecond each ( Table 1). The carboxyl groups of one of the three Asp55 residues at each of the type II calcium binding sites were protonated in the two simulations without Ca 2z , effectively simulating the capsid at a slightly acidic pH to mimic the conditions of the expanded capsid in the 7.5 Å crystal structure [11]. The RNA molecule was not included in our simulations since it cannot be modeled completely in the electron density maps [14,15]. The aim of this work was to probe the dynamic behavior of a virus capsid over timescales that are more than an order of magnitude longer than what has been reported from simulations of viruses previously [16][17][18], and therewith to investigate the role of the structural calcium ions in initiating the dissolution of the satellite tobacco necrosis virus. We particularly looked into the structural features facilitating breaking up of the capsid associated with virus infection.

The capsid expands upon calcium removal
The capsids were stable and remained intact throughout all four trajectories. Removing the Ca 2z had a pronouced effect on the capsid radius. An increase in the radius of gyration (R G ) of 2.6% in Sim NoCaDneut and of 2.4% in Sim NoCaDphys was found over the course of 1 ms. The two trajectories with bound Ca 2z (Sim CaDneut and Sim CaDphys ) showed a weak tendency to increase in size, 0.78% and 0.58% respectively (Table 1 and Figure 2).
The capsids retained their overall spherical shape with just some minor degree of elongation. Figures 3A and 3C emphasize the local anisotropy in the structural changes. The largest changes occurred in the parts of the shell close to the icosahedral 3-fold axes. This area accommodates four calcium binding sites and all of them have charged carboxyl groups as ligands. The charge repulsion induced the formation of small water-filled cavities between the proteins at this protein/protein interface ( Figure 4D). An analysis where the root mean square deviation (RMSd) from the crystal structure for protein dimers, trimers and pentamers was computed is presented in Table 1. The trimers have clearly larger than average RMSd, whereas dimers and pentamers form very stable complexes. This effect is more pronounced in the simulations without Ca 2z .

Dynamics of the coat protein
The secondary structure elements and the overall fold of all the coat proteins were stable throughout the simulations. The Nterminal arm was the most flexible element (Figures 1, 3 and 5A). The a-helix in the N-terminal arm observed in the crystal structure was stable: at the end of the simulation the number of a-helical  DR shell : Increase in the radius of gyration (C a of res. . Crossings: Water permeation events during 1 ms. P f : Osmotic water permeability coefficient. p f : Single pore osmotic water permeability coefficient for the 3-fold axis. a nterm : Number of a-helix residues in the N-terminus. HB Nx : Number of intermolecular hydrogen bonds between pairs of subunits related by 2-, 3-or 5-fold symmetry (disregarding interactions between the N-terminal 11 residues). RMSd Nx : RMSd from the crystal structure for 2-, 3-or 5-protein oligomers calculated for the main domain (C a of residue

Author Summary
We have studied the capsid of satellite tobacco necrosis virus using large scale molecular dynamics simulations, where the atomic motions of 1,2 million particles were tracked over one microsecond. We find that the capsid swells in the simulations, and that the permeability for water increases 10-fold upon removal of the structural calcium ions. The water leaks in predominantly near the three-fold symmetry axis, suggesting that this is the spot where capsid dissociation is initiated following infection.
residues in the N-terminal arm was close to 11, slightly lower in the two trajectories without bound Ca 2z than the simulations with bound Ca 2z (Table 1). Residues 1-11 did not show any propensity to stay in a helical conformation. These residues were modeled as a helix in the starting structure, but they progressively lost that structure. This might be a result of the absence of the RNA molecule since the addition of molecules that mimic the phosphate backbone of RNA has been shown to promote formation of helices in the positively charged N-terminus of CCMV, that presumably plays a similar RNA-binding role [19,20]. The number of intermolecular protein-protein hydrogen bonds was slightly lower in the STNV capsids without Ca 2z , in particular the number of hydrogen bonds between pairs of 2-fold and 3-fold related subunits decreased ( Table 1).
The C a atoms in the shell domain moved on average 0.27 nm in Sim NoCaDneut , 0.26 nm in Sim NoCaDphys , 0.15 in Sim CaDneut and 0.14 in Sim CaDphys predominantly due to the overall radial expansion of the capsid ( Figure 5A). The RMSd after fitting each protein to the crystal structure individually was small (ƒ0.1 nm) apart from the termini and some of the loop regions ( Figure 5B). The flexibility was highest in the two termini, in the loops between secondary structure elements and in the short helix centered at residue Thr119 ( Figure 5C). The bD-bE loop (using the same nomenclature as in [2]) centered at residue Thr80 and the bH-bI loop centered at residue Leu180 show high flexibility both with and without bound Ca 2z ( Figure 5C). Both of these loops are located close to the 5-fold axis, facing the exterior of the capsid ( Figure 1). The bC-bD loop at residue Asp55, the bG-bH loop at residue Thr160 and the C-terminus on the other hand show an  increased flexibility mainly in the two capsids that had no Ca 2z ( Figure 5C). These elements are all located close to the 3-fold axis and Asp55 is the Ca 2z ligand for the type II sites.

Calcium binding sites
When the Ca 2z ions were removed from the capsid, the type III calcium binding sites on the icosahedral 5-fold symmetry axis were populated by either Na z (Sim NoCaDphys only) or water. Sim NoCaDphys had 7 out of 12 type III binding sites occupied by sodium ions throughout the simulation. The remaining type III sites as well as all the type III sites in the capsid of Sim NoCaDneut were occupied by one or two water molecules. The sodium ions were bound at the same position as the calcium ions -at the geometric center of the five Thr138 carbonyl groups -but half of the time, one of the five carbonyl groups pointed away from the 5fold axis and engaged in hydrogen bonds with other water molecules. Water molecules bound to the type III sites were located on the 5-fold axis slightly exterior or interior to the binding site for cations such that they could make hydrogen bonding interactions with two or three carbonyl groups. In these cases the 5-fold symmetry of the protein subunits was broken with one (when two water molecules were bound) or two (when one water molecule was bound) carbonyl groups turned away and engaged in hydrogen bonds with water molecules from the bulk.
The two types of calcium binding sites close to the icosahedral 3-fold symmetry axis were less stable than the type III sites upon Ca 2z removal as reflected in the higher RMSd and RMSf of residues close to the 3-fold axis ( Figure 5). In Sim NoCaDphys all type II sites were occupied by one or two sodium ions subsequent to the equilibration simulation (in which the protein was restrained from moving away from the crystal structure). However, at the end of the production simulation only 1 of the 20 sites was intact with all three carboxyl groups of Asp55 coordinating a sodium ion. The rest of the type II sites were broken up with the carboxyl groups facing other directions( Figure 4D). The type I site also lost its structure in the calcium-free capsids. The four residues that contributed oxygen atoms to coordinate the calcium ion had high RMSf ( Figure 5C) and did not retain their relative positions. The negatively charged carboxyl side chains of all the type I sites together bound approximately 30 sodium ions in an irregular fashion.

Water permeability
The capsids were permeable to water in all trajectories. Assuming a homogeneous permeability across the entire surface of the capsid, an osmotic water permeability coefficient, P f , was calculated [21]. With Ca 2z bound there were about 10,000 permeation events in either direction and the capsid had an average permeability coefficient P f = 2:6|10 {3 cm/s. After removing the Ca 2z P f increased to about 30|10 {3 cm/s ( Table 1). The water transport resulted in a net inflow of water in conjunction with the swelling.
In order to detect structural features associated with water permeability, the crossing events were mapped onto the virus particle. For each water molecule traversing the width of the capsid shell successfully, the closest protein residue (C a ) was determined for the water molecule half-way (in time) through, and statistics of permeation per residue were gathered. In the simulations with bound Ca 2z the permeation mainly occurred between the icosahedral 2-fold and 3-fold symmetry axis, at the junction between three subunits ( Figures 6A-C and 4). The protein-water contacts suggest that this potential water pore is lined by five motifs: the short helix centered at residues 118-119, residues 156-158 in the bG-bH loop, the end of the bC-bD loop, the beginning of the bD-strand and to some degree residues 25-30 close to the flexible loop connecting the N-terminal arm to the shell domain ( Figure 5D). In the simulations without bound Ca 2z the permeability increased at this site, as well as at the protein/protein interface at the 3-fold symmetry axis ( Figures 5D,  6D-F).
If we consider the region around each icosahedral 3-fold axis to be a water pore, we can estimate the single pore permeability coefficient, p f , to be 1:1|10 {14 cm 3 =s in the capsid with bound Ca 2z and 13|10 {14 cm 3 =s without these ions (Table 1).

Discussion
Many plant viruses employ the extremely low concentration of Ca 2z in the cytoplasm of their host cells (the homeostatic concentration of free Ca 2z is between 100-200 nM [22]) as a cue to initiate the replication stage of their life-cycle. The capsid will be effectively depleted of calcium as the equilibrium between occupied and unoccupied sites shifts to the latter.

Capsid swelling
Four different crystals of STNV treated with EDTA were investigated by Montelius et al., but the capsid expanded in only one of those [11] and that particular crystal diffracted to a resolution of 7.5 Å only. In our simulations we could observe a higher degree of fluctuations as well as loss of the icosahedral symmetry in the calcium-deprived capsids compared to the ones with bound calcium. Both higher flexibilty and a lower degree of symmetry would contribute to a lower crystal quality and could explain why there have not been any successful attempts to solve a high-resolution structure of an expanded capsid so far. The trajectories did not quite reach an equilibrium swollen state. The radius of the calcium-depleted capsid can be extrapolated to an ''equilibrium radius'' R g (?) by fitting the R g (t) curve ( Figure 2) to R g (t) = R g (?)2(R g (?)2R g (0))exp({t=t). With this fit an equilibrium radius of &8:22nm is predicted, about a quarter of a percent higher than the value at the end of the Sim NoCaDneut simulation. The weak increase in the radius of the calciumcontaining capsids is probably an artifact due to the lack of the RNA molecule and confirms the long-term unstable nature of the RNA-free capsid.
Swollen STNV capsids can be returned to their native radius by lowering the pH [8], suggesting that the mechanism behind the swelling is electrostatic repulsion of charged aspartate and glutamate ligands. A similar effect was deduced from electrostatics calculations of the native and expanded structures of the CCMV capsid [23]. The protonation state of the carboxyl groups is probably coupled to the magnitude of the expansion. The crystals of the expanded STNV formed at pH 6.5 and the expansion was only moderate, so the capsid may have been protonated at one of the calcium binding carboxyl groups. Since the type II site has the highest local density of negative charge, we decided to protonate one of the three carboxyl groups at each type II binding site in our simulations of the calcium-free capsid.
Analytic ultracentrifugation measurements estimate that the STNV particle can expand up to 7% when treated with EDTA [8]. The crystal structure of the expanded virion showed a radial expansion of between 0.1-0.4 nm, equivalent to 1%-5% (higher closer to the icosahedral 3-fold axis). This agrees well with what we observed in the simulations: an average increase in the R G of 2.5% ( Table 1) and peak RMSd values of the C a of the shell domain above 0.35 nm at the 2-fold subunit interface, e.g. residues 96 and 122, and at the 3-fold axis, e.g. residues 25 and 195 ( Figure 5A).

Permeability coefficients
The osmotic water permeability P f of the capsid is comparable to lipid membranes, which usually have P f = 2{16|10 {3 cm/s [24]. The protein shell is thinnest at the 5-fold axis, but the calcium ion at the type III site is the one that is most difficult to chelate with EDTA [25], something which was corroborated by spectroscopic measurements showing that the type III site has a remarkably high binding affinity to the Ca 2z analog Eu 3z (1:1+0:3 nM) [26]. In the simulations of the capsid without calcium, the coordinating cage of oxygen ligands was best preserved at the type III site. Sodium ions or water molecules replaced the divalent ion on the icosahedral 5-fold axis and obstructed the opening, resulting in minimal permeability even without Ca 2z . Instead, we observed extensive water permeation on the 3-fold axis and in a region close to the 3-fold axis between the 3-fold and the 2-fold axis (Figure 4). The cluster of four calcium binding sites around the 3-fold axis contains several carboxyl groups from aspartate and glutamate residues. Removing the calcium ions introduces a large amount of net negative charge that caused the subunits to move apart, creating water pockets at the 3-fold axes causing increased water permeability. If the entire region around the 3-fold axis is considered a water pore, the permeability of it is comparable to that of membrane proteins that function as water pores, e.g. mammalian aquaporins that have reported single pore permeability of 2{25|10 {14 cm 3 =s [27].

Comparison to the expanded structure of TBSV
The difference in size and triangulation number makes it difficult to compare the expanded structure of TBSV and STNV. The capsid of TBSV consists of 180 identical subunits in a T = 3 arrangement where each of the 60 icosahedral asymmetric units consists of three proteins in slightly different configuration. At the center of these three proteins, there is a so called quasi-3-fold axis that relates three approximately equivalent protein positions [28]. The six calcium binding sites of TBSV are located pairwise between pairs of subunits in the asymmetric unit and each site has five acidic residues from both proteins. The swollen TBSV was crystalized at pH 7.5 and the structure is about 7% larger than the native one. The most predominant structural change is that large openings appear between the quasi-3-fold and 2-fold related subunits [10]. The 3-fold axis of STNV resembles the quasi-3-fold axis of TBSV. The interfaces between the proteins around these axes contain six (TBSV) respectively four (STNV) Ca 2z binding sites (Figure 1), coordinated by carboxylic groups. In both capsids the largest structural differences between the expanded and native structures can be found here.

Factors affecting capsid stability
Not all viruses that have calcium binding sites in the capsid show a swelling behavior. The rhinoviruses have a calcium binding site on the 5-fold axis [7], but do not swell upon removal of these ions. Interestingly, this binding site show striking similarities to the type III calcium binding site in STNV. In both cases five backbone carbonyl oxygens point at the 5-fold axis where the ions are bound. The 5-fold axis is a region with high degree of symmetry constraints and putting an ion there solves the problem of fulfilling the symmetry at a very congested interface. We therefore propose that the carbonyl type of binding sites have a more structural role, while the carboxyl type of binding sites may play a role in the dissolution of the capsid upon infection. This would explain the low degree of structural change at the 5-five fold site in our simulation and their high propensity to bind other cations or water molecules, and this is in line with the relatively higher affinity for ions at these sites [25,26]. There are no atomic structures of the STNV genome available, although a recently published high-resolution structure shows short fragments of RNA [14]. RNA-free virus-like particles of the STNV coat protein have not been observed so far, but the coat proteins of similar viruses readily form empty shells [29,30]. Rather than inserting a modeled RNA structure into the capsid we decided to perform the simulations without RNA. The fact that RNA nevertheless may contribute to the stability of the capsid is evidenced from our simulations at physiological salt concentration, providing more shielding of electrostatic interactions. The radial expansion is less pronounced in this case ( Figure 2) and the number of protein-protein hydrogen bonds seems to increase slightly (Table 1), a typical ''salting out'' effect.
Previous all-atom molecular dynamics simulations of capsids and virus-like particles focused on specific properties like the mechanical strength using force probes [17,18] or the effect of hydration on coherent diffraction from single virus particles [31]. It has become clear though that virus simulations are sensitive to the starting conditions, like a balanced amount of water on the inside, and they require long equilibration times [18,32]. Therefore we paid special attention to the preparation of the starting structures with an extensive equilibration protocol with part of the counter ions specifically added to the interior cavity to avoid a sudden influx of solvent ions that could disrupt the proteinprotein interactions and by carefully balancing of the hydrostatic pressure inside and outside of the capsid. By doing control simulations with Ca 2z and be varying the ionic strength we can draw firm conclusions on the effect of Ca 2z removal on STNV.

Conclusions
The simulations presented here illustrate the mechanism by which an entire virus capsid can transform from a closed configuration into an open one with significantly increased water permeability. The magnitude of the expansion in the simulations is in good agreement with experiments. The higher flexibility and the degradation of the icosahedral symmetry in combination explain why it has been difficult to crystallize expanded capsids. Our work strongly suggests that there are two types of calcium binding sites, playing different roles in the virus lifecycle. The binding site on the 5-fold axis has a more structural role and is less involved in the capsid expansion. The binding sites between 2-and 3-fold related subunits contain many charged carboxylic side chains. Dissolution of the capsid is initiated here due to the electrostatic repulsion between these residues, if the Ca 2z ions are removed.

Materials and Methods
All preparations, simulations and analyses were performed with the GROMACS simulation package version 4.5.3 [33] unless otherwise stated.

System preparation
The starting structures were prepared from the X-ray crystal structure of the coat protein of STNV (PDB ID: 2buk) [2]. Since residue 1-11 can not be resolved in the crystal structure, these were modeled as an a-helix in a direction that did not cause steric clashes with neighboring proteins. The full capsid was generated by applying the icosahedral rotation-translation matrices in the PDB file. One Ca 2z ion was kept at each of the 92 calcium binding sites. The Amber99sb-ILDN forcefield [34,35] was used for the protein combined with TIP3P water [36] in a rhombic dodecahedron simulation box with a side of about 25 nm. Water molecules on the inside of the capsid were replaced with Cl { ions to obtain a neutral system (Sim CaDneut ). Additional Na z and Cl { ions were added to obtain a system with an approximately physiological ion concentration of 150 mM (Sim CaDphys ). The calcium ions were removed, a proton was added to one Asp55 side chain at each of the type II calcium bindings sites and the number of counter ions was adjusted to obtain two neutral calcium-free systems (Sim NoCaDneut and Sim NoCaDphys respectively). Each system consisted of roughly 1.2 million atoms (Table 1).

Preparatory simulations
The starting structures were energy minimized and subsequently the solvent was equilibrated for 10 ns while restraining the protein atoms and the calcium ions to the crystal coordinates with harmonic potentials with a force constant of 1000 kJ mol {1 nm {2 . During the equilibration a 2 fs integration timestep was used and the neighbor lists were updated every 10 th timestep. Short range nonbonded Van der Waals (Lennard-Jones) and Coulomb interactions were calculated within a cut-off radius of 1.15 nm. The long range electrostatic interactions were calculated with the particle mesh Ewald (PME) method [37] with a grid spacing of 0.133 nm. The long range Lennard-Jones interactions were analytically corrected for in the calculation of the pressure and the energy. The pressure of the simulation box was kept at an average of 1 bar using the isotropic Berendsen barostat [38] with a time constant of 25 ps and a compressibility of 4:5|10 {5 bar {1 . The solvent and the capsid were coupled separately to an external heat bath at 300 K with the velocity-rescaling thermostat [39] using a time constant of 0.5 ps. Water molecules were constrained using the SETTLE algorithm [40] and the covalent bonds in the proteins were constrained using the P-LINCS algorithm [33]. Boundaries were treated periodically.

Production simulations
After the equilibration, the position restraints were removed and an integration time step of 4 fs was used to generate 1 ms trajectories. In addition to constraining bond lengths, virtual hydrogen atoms were used [41] which allows slightly longer time steps. The isotropic Parrinello-Rahman barostat [42,43] was used to keep the average pressure at 1 bar with a time constant of 1 ps. The calcium ions were tethered to the oxygen ligands using harmonic potentials with force constants of 5000 kJ mol {1 nm {2 . All other simulation parameters were the same as during the equilibration. The trajectories were sampled every 50 ps for analysis. The production simulations were calculated in parallel on a Cray XE6 system over 2016 cores (612 of which for PME calculations) at a speed of 30 ns/day.

Analysis
Unless otherwise stated, all trajectories were analyzed at 1 ns intervals and final values were calculated as the average over the last 100 ns. Root-mean-square displacement (RMSd) and fluctuation (RMSf) was calculated as: The number of a-helical residues was calculated using the g_helix program of the GROMACS package. In the water permeability analysis water molecules were counted when they had traversed the entire width of the capsid shell. Visual inspection did not imply a single-file type of permeation mechanism, which justifies that the distinction between a diffusive and osmotic permeation coefficients was not required [21]. The osmotic water permeability coefficient, P f , and the single pore osmotic water permeability coefficient, p f , was calculated as the average permeability in both directions using these formulas: Where N is the number of permeation events, t is the duration, A is the area and DC is the concentration gradient, i.e. 55 mol/L for pure water. The capsid was approximated to have the same area as a sphere with radius 8 nm.