Physical properties of the HIV-1 capsid from all-atom molecular dynamics simulations

1The family of retroviruses is characterized by their ability to incorporate viral DNA into a host cell’s genome. Most retroviruses infect cells during mitosis when the chromatin is exposed to the cytoplasm. Conversely, the genus of lentiviruses, like the human immunodeficiency virus type 1 (HIV-1), have evolved to infect nondividing cells. Since the host cell’s chromatin is protected by the nucleus, the HIV-1 replication process requires coordination between reverse transcription of viral RNA and nuclear import. Viral RNA is encased in a shell made of the capsid protein CA. Originally thought to play a trivial role in the infection process, it is now well established that the viral capsid fulfills several essential functions. In particular, capsid involvement in the prevention of innate sensor triggering, regulation of reverse transcription and regulation of the nuclear import pathway is of central importance to the successful infection of a host cell.

2Structurally, CA consists of a single monomer composed of two independent domains, the C-terminal domain (CTD) and the N-terminal domain (Fig. 1a). In solution, CA readily forms dimers connected by their CTDs. At higher concentrations of CA and at high salt concentrations (>1 M NaCl), the protein assembles into elongated tubes. Such tubes are made of the hexameric form of CA, arranged in a honeycomb lattice. Assembled HIV-1 capsids are highly asymmetric cone-shaped structures, with sizes ranging from 100 to 200 nm long and 45 to 50 nm wide (Fig. 1b), that are made of pentamers and hexamers following Eberhard’s theorem. The location of the pentamers induces high curvature and permits the closure of the capsid. Indeed, in the HIV-1 capsid seven pentamers are located at the base and five pentamers at the tip (Fig. 1b).

3The function of a CA monomer is to assemble into a complete capsid, while the function of the assembled capsid is to infect a living cell; the chemical–physical properties of the assembled capsid are essential to elucidating its biological function. In the present study, the investigation of the chemical–physical properties of the HIV-1 capsid without genome require the simulation of an entire virus-like particle (VLP), namely a 1.2 μs, 64 million atom molecular dynamics (MD) simulation (Fig. 1b). Analyses of the trajectory reveal atomic details as well as global emergent properties of an entire HIV-1 VLP. In particular, the simulations reveal spontaneous translocation of water and ions through the capsid core, a molecular process related to reverse transcription. Furthermore, analysis of the electrostatic potential produced by the empty capsid reveals electrostatic equivalence between the interior and exterior of the capsid. In addition, analysis of the distribution of ions around the capsid reveals a binding pattern and specific binding sites on the capsid for chloride and sodium. Finally, the mechanical vibrations and normal modes of the empty capsid reveal possible pathways related to capsid uncoating. Altogether, our findings set the stage for a better understanding—from an atomic-resolution perspective—of the emergent properties of the HIV-1 capsid and their relationships to the infective cycle.

4We began our analysis by assessing the stability of the HIV-1 capsid model, formed by 1,300 copies of the CA, shown in Fig. 1. Stability of the individual constituents of the capsid, namely CA pentamers and hexamers, was evaluated by means of root mean squared deviation (r.m.s.d.) from the starting structure over the course of the simulation (Fig. 2a). The r.m.s.d. values were calculated excluding parts of the protein known to be disordered or mobile, such as the cyclophilin-A (cypA) binding loop and residues 210 to 220 directly after α-helix 11. The average r.m.s.d. for CA hexamers and pentamers of 3.5±1.1 Å shows the local stability of the capsid. Furthermore, the measured r.m.s.d. over the trajectory for CA pentamers and hexamers of 2.7±0.6 and 3.8±1.0 Å, respectively, shows that the pentamers are on average more rigid than the hexamers. Similarly, the variability within the ensemble of hexamers in the capsid was 3.7±0.9 Å as compared with 2.6±0.5 Å for the 12 pentamers. In addition, the stability of the capsid as a whole was evaluated by calculating the r.m.s.d. matrix between pairwise structures over the course of the trajectory (Supplementary Fig. 1 and Supplementary Movie 1); the distance matrix shows that after 400 ns all capsid structures are within 3.5±1.0 Å of each other.

5Changes in the cross-sectional area and height of the capsid are related to its global stability (Fig. 2b–d). Therefore, both the height and cross-sectional area were calculated along the three principal axes of inertia of the capsid (Fig. 2b). Interestingly, during the first 400 ns of simulation, a shrinking of the capsid was observed, evidenced by a reduction in both the height and cross-sectional area at a rate of 0.025 nm ns−1 and 0.113 nm2 ns−1, respectively (Fig. 2b,c). After 400 ns, the capsid height and cross-section reach a plateau for the remaining 800 ns of the simulation that is used in the subsequent analyses.

6Examination of the water density in the interior and exterior of the capsid reveals the origin of capsid shrinkage (Supplementary Fig. 2a). During the first 100 ns of simulation, a decay in water density was observed inside the capsid, followed by stabilization around 1.000±0.002 gm cm−3. Importantly, the rates of water moving in and out of the capsid are 20.1 × 103±1.6 × 103 and 19.7 × 103±1.5 × 103 water molecules per ns, respectively (Supplementary Fig. 2b). The transfer rates of water indicate thermal equilibrium between the interior and exterior of the capsid and imply that the HIV-1 capsid is capable of replacing all of its contained water in ∼1 μs. Remarkably, the transfer rates observed for HIV-1 are in stark contrast with order of magnitude smaller transfer rates observed for poliovirus capsids.

7The capsid protein is negatively charged; in fact, at pH 7, the net electric charge of the whole capsid is −3,528 e. The electrostratic potential due to the distribution of electric charge of a fully solvated capsid, that is, explicitly including water and ions, was calculated at 20 ps intervals over the course of the simulation, over a grid containing 350 × 353 × 600 bins along the, and direction, respectively. The resulting electrostatics maps were averaged over the last 400 ns of simulation. Snapshots of the electrostatic maps at regular intervals over the trajectory, averaged over 10 ns of simulation, are presented in Supplementary Fig. 3. The time-averaged electrostatic potential is shown in Fig. 3a,d. Notably, a gradient of up to 7 V is observed between the CypA binding loop and the inner core of CA (Fig. 3c); in addition, the innermost layer of the VLP constitutes an isopotential volume spanning the entire capsid (Fig. 3b). Intriguingly, the effective electrostatic potential of hexamers and pentamers is remarkably similar. Furthermore, the regions of the solvent inside and outside of the capsid are at the same electrostatic potential as shown in Fig. 3b.

8The presence of water and ions in the simulation box offers a unique view into the interactions between the capsid and its native environment. Taking advantage of the small variations within each of the ensembles of hexamers and pentamers, 186 hexamers/12 pentamers were r.m.s. fitted to a common reference; then, occupancy of chloride and sodium ions was measured over the last 400 ns of simulation, resulting in a total sampling per hexamer or pentamer of 186 × 0.4 μs=74.4 μs and 12 × 0.4 μs=4.8 μs, respectively. The analysis reveals regions of high occupancy for chloride (cyan) and sodium (yellow) (Fig. 4). Comparison of our results with published crystal structures show that our observation of chloride at the centre of the hexameric rings is accurate. However, in our simulations, the binding of chloride to other regions is also detected, in particular assisted by highly conserved residues K70 and K180 (ref. ). In addition, we also observe the presence of sodium near the surface of the capsid at multiple sites, especially near genetically important residues E71, E75, E76, E79, E212 and E213.

9Binding of ions to the whole capsid reveals an interesting overall pattern (Fig. 4): Chloride ions form an inner layer adjacent to the surface of the capsid, while sodium binds to the exterior of the capsid. Examination of the transfer rates for both ionic species present in the simulation reveal that ions are in thermal equilibrium between the interior and exterior of the capsid (Supplementary Fig. 2c,d). Interestingly, over the course of the simulation, we were able to observe numerous translocation events of chloride through the central pores in hexamers and pentamers (Supplementary Movie 2), and sodium through the pores located between adjacent CTD dimers (Supplementary Movie 3). Notably, the inwards and outwards transfer rates for sodium are 9.4±2.4 and 8.5±1.1 molecules per ns, respectively. In contrast, we observed inwards and outwards transfer rates twofold that of sodium for chloride, 22.2±2.9 and 20.8±3.3 molecules per ns, respectively.

10At each binding site, ion occupancies are directly related to the binding free energies of both ionic species; therefore, from the occupancy maps, we derived a three-dimensional (3D) potential of mean force (PMF) for both pentamers and hexamers, as explained in Methods and shown in Fig. 4. The resulting PMFs illustrate that irrespective of the oligomeric (hexamer or pentamer) state of CA, the ion binding pattern is remarkably similar. Nevertheless, the PMFs show that binding of chloride in the central pore, near R18, is weaker for hexamers than it is for pentamers, as demonstrated by their Gibbs free energy =−1.5±0.9 kcal mol−1 and =−2.7±1.1 kcal mol−1, respectively. The difference in binding energies between pentamers and hexamers suggests that chloride ions translocate more easily through hexameric channels than they do through pentameric channels. Similarly, sodium binding at the hinge region (residues 144 to 148) and the threefold interface between hexamers exhibits a lower affinity, =−1.8±0.5 kcal mol−1 compared with chloride in any other binding site. Remarkably, the free energies of binding for chloride and sodium ions are similar to the experimentally measured value for binding of sodium to wild-type thrombin (−2.3 kcal mol−1); nonetheless, we find lower values comparing the free energy of binding of sodium to glutamate receptors (−10.0±5.0 kcal mol−1), and the sodium symporter (−10.0±2.0 kcal mol−1). Overall, the distribution of ions and exchange rates reveal that the HIV-1 capsid is able to translocate ions with remarkable specificity; such specificity may play an important role during infection, as it would assist and filter translocation of molecules that are key for successful infection (for example, DNA nucleotides).

11Structural fluctuations of the capsid during the simulation revealed complex dynamics. Deviations from the equilibrium position of individual Cαs were projected over the capsid for different time points, as described in Methods (Fig. 5a and Supplementary Movie 4). The projection of the variations reveals an oscillatory behaviour of the surface of the capsid with magnitude 3.0±1.0 Å. Interestingly, similar oscillations have been observed for lipid vesicles. These ‘surface waves’ are spread over the capsid and dynamically correlate large regions of the capsid (Fig. 5a and Supplementary Figs 4 and 5).

12The capsid is a closed container that can be characterized based on its acoustic properties. For this purpose, periodograms were calculated for each Cα time series (Fig. 5b). A periodogram is the squared magnitude of the discrete Fourier transform and identifies the most dominant frequencies in a time series. Therefore, the periodograms permit one to identify the most dominant frequencies in the measured time series for every Cα in the HIV-1 capsid. Because of the large number of Cαs (>150,000,) analysing the periodograms by visual inspection is impractical. Instead, by using each periodogram as a descriptor for every fifth Cα and the Manhattan distance between descriptors as a similarity score between Cαs, every fifth Cαs were classified into clusters using a partition around medoids (PAM) algorithm. The silhouette criteria (Supplementary Fig. 6) was used to establish the optimal number of clusters, indicating four of them (Fig. 5b). The same clustering analysis was performed for halved and quartered fragments of the simulation that also resulted in four major clusters (Supplementary Fig. 4).

13Remarkably, each of the four clusters identified by the PAM analysis are dominated by a different fundamental frequency (Fig. 5b and Supplementary Fig. 4). Indeed, the two largest clusters (N=10,743 and N=11,898) are dominated by the two lower fundamental frequencies (2.38 and 4.76 MHz); conversely, the two smallest clusters (N=4,679 and N=2,786) are dominated by two higher fundamental frequencies (7.14 and 11.9 MHz). The clusters were mapped to their locations on the structure of the capsid (Fig. 5c,d). Interestingly, the highest frequency cluster (four) is located at the tip of the capsid and around the circumference of the base. On the other hand, lower frequency clusters are located on the large regions spanning between the base and tip of the capsid.

14As previously mentioned, analysis of halved and quartered fragments of the trajectory reveals that oscillations of the capsid can still be clustered for shorter lengths of the simulation. Nonetheless, the mapping of the frequencies on to the capsid structure identify different frequency regions for different fragments of the simulation (Supplementary Fig. 5). These differences between frequency-to-structure mappings suggests that the waves we observe are not stationary. We therefore conclude that the waves observed behave like capillary waves common in several types of membranes and fluids. A remaining question, beyond the scope of the present study, is the relationship between capillary waves and permeation of water and ions through the capsid.

15Normal modes are informative to elucidate the collective motions of flexible molecules that underlie changes in their conformation. The displacement patterns obtained from normal-mode analysis (NMA) provide insights into the mechanical nature of biological nanomachines and can be used to derive conformational transitions related to biological function. A common practice is to extract slow motions from MD trajectories using principal component analysis (PCA) that reveal principal components (PCs) akin to normal modes. In the context of viral capsids, several studies have been performed to extract such collective motions by using a plethora of methods. Nonetheless, since all NMA studies of entire viral capsids were performed in vacuo, a long-standing question in physical virology is the influence of hydrodynamic effects resulting from the presence of solvent in the vibrational modes of viral capsids.

16The vibrational modes of a fully solvated HIV-1 capsid were obtained in the present study based on all-atom MD simulations. In particular, PCA was employed by performing singular value decomposition of the last 400 ns of the trajectory (described in detail in Methods). The PCs were calculated using the Cα atoms (N=150,528), thus resulting in a total of d=3 × N−6=451,578 degrees of freedom. Convergence of the subspace spanned by the PCs was evaluated for different numbers of modes and lengths of the trajectory using the mutual similarity measurement, as explained in Methods and illustrated in Supplementary Fig. 7. A total of n=20,000 frames were employed for the PCA calculations. The number of PCs that can be extracted from the simulation is given by min(d,n)=20,000 (ref. ). Interestingly, 300 modes are required to account for over 80% of the variance observed in the simulations (Fig. 6a); similarly, only 100 modes are neessary to account for over 60% of the variance (Fig. 6a). Projection of the last 400 ns of the simulation onto the first two PCs further demonstrates the stability of the capsid (Supplementary Fig. 8). In addition, for the lowest 800 modes the percentage of the capsid involved in each mode was quantified using the so-called participation number. The participation number reveals that the modes are not localized, and that the PCs involve large portions of the capsid ranging from 70 to 84% of the entire capsid (Fig. 6b).

17Projection of the displacement patterns of the first and second PCs (which together contribute 25% of the overall variance) onto the structure of the capsid reveal complex motions of the capsid subunits (Fig. 6c,d and Supplementary Fig. 9, respectively). The motions are represented as vectors in the figures, indicating the displacement of a given atom. The pattern of the motions is intriguing, as a belt can be observed that effectively divides the capsid into two hemispheres, corresponding to the tip and the base of the cone (Fig. 6d). Additionally, related to the capillary waves observed in the previous section, vortices are observed in the motions throughout the structure of the capsid (Fig. 6c,d). Remarkably, the first 100 PCs reveal the absence of a distinct radial breathing mode, and show that all modes exhibit a combination of displacements that imply motion both radially and within the surface of the capsid. The first 100 PCs are included in Supplementary Dataset 1 that can be visualized using the NormalModeWizard available in visual molecular dynamics (VMD).

18The HIV-1 capsid is a finely tuned nanodevice that modulates several molecular events during HIV-1 infective cycle by interacting with multiple cellular host factors. The structural stability of the capsid during different stages of the infective cycle remains poorly understood. In the present study, despite observing a global shrinking of the capsid compared with our original model, the overall shape of the capsid and the structure of CA hexamers and pentamers remained largely unaltered over 1.2 μs of simulation time. Such global and local stability further supports the fullerenic model of the HIV-1 capsid, as well as the stability of the capsid interfaces, formed by dimers and hexameric/pentameric rings. Our observation that a shrinkage of the capsid is followed by a stable structure is consistent with coarse-grained simulations of HIV-1-like particle assembly, where a global arrangement of the VLP occurs followed by local re-arrangement of the VLP structure. In addition, our results suggest that a closed capsid, irrespective of its cargo, once assembled, remains stable and closed while still permitting the passage of water and ions.

19So far, acoustic and normal mode analysis for large macromolecular systems have been limited to symmetric particles or coarse-grained representations. In the present study we determined acoustic properties and PCs from all-atom simulations of a solvated capsid. We found that the capsid exhibits an oscillatory behaviour manifested as surface waves. Our analysis indicates that the fundamental frequencies associated with such waves occur in the ultrasound regime (2.38 to 11.9 MHz). To our surprise, the fundamental frequencies were associated with different regions of the capsid. Interestingly, the tip of the capsid, where pentamers are closer to each other, presents the highest fundamental frequency. In addition, our NMA results suggest that the capsid exhibits long-range collective dynamics. Such correlations over long distances may play important roles during infection, as they permit the capsid to transfer information to/from distant regions of its structure. This allosteric communication across the capsid could serve a function during nuclear import, during which interactions between the nuclear pore complex and the capsid modulate infection (Supplementary Fig. 10).

20A significant result from the present study is the discovery of binding of cations and anions in genetically important regions of the assembled capsid. We are not the first group to observe ions bound to CA in high-order assemblies; in addition to water molecules, chloride ions were detected by X-ray crystallography at multiple interfaces of a planar hexameric lattice (including the central pore, the CTD dimer and the threefold symmetry axis). However, whether chloride or sodium ions play a biological role is unclear, as they may be substituting for other ions available in the cytoplasm. Furthermore, not all of the binding sites we observe have also been observed experimentally; particularly, sodium has not yet been observed to interact with the hexameric lattice in the X-ray structures. Moreover, the observed electrostatic signature of the assembled capsid together with the ion binding sites could indicate the existence of protein–protein interfaces for as yet undiscovered cellular factors. Interestingly, the distribution of charge on the surface of the capsid of another retrovirus, the Rous sarcoma virus, is essential for successful infection. Finally, our observation of ion-specific binding sites could also explain the requirement for high-salt concentrations for in vitro assembly of CA tubes and VLPs.

21Remarkably, we also observe that ions translocate through chloride- and sodium-specific channels within the surface of the capsid. While homeostasis could be a simple explanation for ion translocation, there are other biological functions of the capsid that could require the existence of such channels. For instance, reverse transcription (which is coupled to capsid uncoating) requires that DNA nucleotides be able to move from the exterior to the interior of the capsid. Therefore, it is possible that the channels found within the capsid serve to translocate small molecules, including nucleotides, necessary during different stages of the infective cycle.

22Importantly, in authentic virions, deviations should be expected from the distribution of ions and electrostatics presented here due to the presence of HIV-1 single-stranded RNA and other proteins in the interior of the capsid. In particular, the native HIV-1 capsid contains two copies of its negatively charged ∼9.7 kb-long genomic RNA that likely induces osmotic and mechanical pressures between the interior and exterior of the capsid; a similar shift in pressure has also been proposed for poliovirus capsids. Remarkably, changes in pressure of the capsid have been observed during reverse transcription in time-lapse atomic force microscopy experiments that result in rupture of the capsid at the narrow end. In the context of in vivo cargo, the ions channels identified in the present study could help regulate the osmotic shock produced during reverse transcription of the genomic RNA.

23The present work is a comprehensive study of the chemical–physical properties of the HIV-1 capsid—an important therapeutic target. Through the combination of state-of-the-art MD simulations with scalable, robust and model-free statistical analysis, we characterized—at atomic resolution—the stability, electrostatics, water/ions permeability and the dynamic and acoustic properties of the HIV-1 capsid. Our results may provide a new avenue for the development of therapeutic interventions that seek to alter the biophysical properties of the HIV-1 capsid towards the treatment of disease.

24In the present study, computer simulations were employed to investigate the stability, electrostatics, ion permeability, acoustics and mechanical properties of the HIV-1 capsid at atomic resolution. In particular, a 1.2 μs MD trajectory of 64,423,983 atoms probed high spatial and temporal resolution characteristics of an empty HIV-1 capsid embedded in its native environment. This level of detail is currently inaccessible to experimental methods alone.

25The HIV-1 capsid composed of 186 CA hexamers and 12 CA pentamers was embedded in a water box with 150 mM sodium chloride, resulting in a simulation box of dimensions 70 nm × 76 nm × 121 nm and a total of 64,423,983 atoms (Fig. 1b), as described in a previous publication and in Supplementary Note 1. In particular, the sequence of HIV-1 subtype B, NL4-3 with the A92E mutation, was used for all CA monomers. For the present study, MD simulations of the unrestrained, fully solvated HIV-1 empty capsid were performed using NAMD 2.10 (refs, ) for a total of 1.2 μs. MD simulations were performed on the Oak Ridge Leadership Computing Facility (OLCF) TITAN supercomputer (INCITE allocation BIO024) using 3880 GPU accelerated Cray-XK nodes; details regarding the computational challenges of the present simulation are described in Supplementary Note 2. The CHARMM36 (ref. ) force field was employed with the TIP3P water model at 310 K and 1 atm. Simulations carried out in the present study used the r-RESPA integrator available in NAMD. Long-range electrostatic force calculations employed the particle mesh Ewald method, utilizing a grid spacing of 2.1 Å and eighth-order interpolation with a 1.2 nm cutoff. The simulations employed an integration time step of 2 fs, with nonbonded interactions evaluated every 2 fs and electrostatics updated every 4 fs; all hydrogen bonds were constrained with the SHAKE algorithm.

26Analysis of the capsid electrostatics and ion permeability from the MD trajectories took advantage of OLCF and National Center for Supercomputing Applications (NCSA) Blue Waters high-performance Lustre filesystem, with a stripe count of 16; details regarding the performance of the Lustre filesystem for analysis of biomolecular simulations have been studied by Stone et al.. Using OLCF Rhea cluster, electrostatic maps were calculated for the entire system, including protein, water and ions, at intervals of 20 ps over the trajectories, resulting in 20,000 maps. Long-range electrostatic calculations were performed using the particle-mesh-Ewald method in VMD. A 3D grid with voxel size of 2 Å was employed for the electrostatics calculations, resulting in 350 × 343 × 600 bins along the, and direction, respectively. The 20,000 maps were then averaged using volutil in VMD. Ion occupancies were calculated at 1 Å resolution (700 × 646 × 1,200 bins) using the volmap plugin in VMD. The occupancy maps were then processed with volutil in VMD to calculate the Gibbs free energy of binding using, where ρ is the probability of finding an ion species at a particular site (voxel) as compared with the probability of finding the same ion species in the bulk ρ0=0.000151 (explained in detail by Cohen et al.).

27To distinguish between the interior versus exterior of the capsid, a 3D ray-tracing method based on the digital differential analysis algorithm was written in C++ and implemented in VMD. A grid with the same dimensions as the simulation box was created for each frame of the trajectory. To identify the regions of the grid occupied by the protein, the molecular surface of the capsid was calculated using the QuickSurf algorithm in VMD, yielding a continuous surface without holes. Subsequentely, for each voxel in the grid, six rays were cast in the, and direction, respetively. If any of the rays hit the external wall of the grid, the voxel was considered to be outside of the capsid; conversely, if none of the rays hit an external wall, the voxel was considered to be inside of the capsid. Voxels assigned to the protein were not evaluated using the digital differential analysis algorithm and were classified as protein voxels. Using the classification of each voxel, atoms within a voxel were labelled as belonging to the inside, outside or protein bound. A grid length of 2.0 Å was employed for all calculations. Water density was calculated using the relationship,, where Mw is the molecular weight of TIP3P water (18.016 g mol−1), V is the volume of the region of interest in (Å3) and N is the number of particles in contained the given region.

28Exchange rates for water and ions were calculated from the number of molecules that were located on one side of the capsid surface at a reference frame t0 and found on the opposite side of the surface at a later frame t>t0 (ref. ). The reference frame was selected at different intervals to minimize the effects of molecule recrossing through the boundary of the capsid (Supplementary Fig. 2). Analysis of the densities and exchange rates from the MD trajectories took advantage of NCSA Blue Waters high-performance Lustre filesystem, with a stripe count of 120. Trajectories were analysed in parallel with a load of 12 frames per compute node.

29First, a reference structure for the capsid was determined using the statistical package R as follows. The similarity matrix containing pairwise r.m.s.d. between every time frame of the capsid simulation for the last 400 ns of simulation was generated (using Cαs); r.m.s.d. was calculated using the Kabsh algorithm. The parallel implementation of the partition around medoids PAM algorithm available in the SPRINT package was utilized to determine the structural medoid of the capsid. Medoids are representative structures observed during the simulation, whose average similarity to all other frames in a trajectory is minimal. Medoids are similar in concept to the mean structure, yet have the advantage of being physical structures, as they are always members of the trajectory. With the medoid as the reference structure, the normalized radial fluctuations for every Cα atom was calculated as, where i=1…NCα. Subsequentely, periodograms were calculated using R for every fifth Cα (Fig. 5b). A periodogram is the squared magnitude of the discrete Fourier transform and identifies the most dominant frequencies in a time series. Periodograms were then clustered using the Manhattan distance as a similarity score and employing the parallel PAM algorithm in SPRINT. The silhouette criteria were evaluated for 2 to 200 clusters, indicating an optimal partioning of the data into four clusters (Supplementary Fig. 6).

30The MD trajectory can be written as a matrix, of dimensions d × n, where d=3NCα=451,578 is the number of cartesian coordinates for all Cαs, and n=20,000 is the total number of frames in the trajectory. The mean position of each atom is then given by the vector. Similarly, the covariance matrix of is given by, with entries, sjk and. Therefore, the scaled-centred trajectory can be calculated using, where Sdiag is the diagonal matrix with entries obtained from the covariance matrix S. Subsequentely, PCA was performed by calculating the singular value decomposition of the scaled-centred trajectory by means of the thick-restart Lanczos algorithm in R. In particular, singular value decomposition calculates the r=min(d, n) eigenvalues λ1…r and eigenvectors/PCs ηi(λ) of the covariance matrix S.

31The resulting ηi(λ) PCs and their associated eigenvalues λ were analysed as follows. The cummulative contribution to variance from the first k PCs was calculated using. The participation number of each PC was calculated using, where is the entropy (information) of the PC. The probabilites pi(λ) are the squared component of the normalized relative displacements for each atom/direction i, pi(λ)=|ηi(λ)|2 (ref. ). A PC equally distributed over every atom results in a participation number equal to the number of atoms in the system. Conversely, small values of the participation number indicate that the PC is not global.

32Convergence of the subspace spanned by PCs was evaluated using the mutual similarity measurement. Mutual similarity was evaluated using the relationship,, where 〈·〉 denotes an ensemble average over the configurations observed in the trajectory, x denotes a protein configuration, denotes the norm and denotes the projection to the m-dimensional PCA subspace. Importantly, the mutual similariy between PCA subspaces is a number between 0 and 1. To evaluate the convergence of the PCA subspaces, PCs were calculated for trajectory fragments of different size, ranging from 50 to 200 ns. Subsequently, mutual similarities, using 10 to 1,000 PCs, were calculated between fragments using R and plotted as a function of fragment size (Supplementary Fig. 7).

33The data that support the findings of this study are available from the corresponding author on request.

34How to cite this article: Perilla, J. R. & Schulten, K. Physical properties of the HIV-1 capsid from all-atom molecular dynamics simulations. Nat. Commun. 8, 15959 doi: 10.1038/ncomms15959 (2017).

35Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.