RBD ACE2 contact analysis

Code and workflow for running MD simulations on Folding@home to analyze RBD:ACE2 contacts

Publication: https://doi.org/10.1016/j.cell.2021.01.037

License

• This software is licensed under the MIT license - a copy of this license is provided as SOFTWARE_LICENSE
• The data in this repository is made available under the Creative Commons CC0 (“No Rights Reserved”) License - a copy of this license is provided as DATA_LICENSE

Manifest

• 01_system_preparation - Contains all scripts and PDBs before and after structure preparation, system parameterization, and equilibration. Also contains the representative glycan structures added to RBD:ACE2
• 02_analysis - Contains the script used to analyze RBD:ACE2 contacts in the trajectories
• additional_data - Contains all scripts and PDBs before and after structure preparation, system parameterization, and equilibration for systems that were not included in the paper, but are relevant to this study.

Contributors

• Ivy Zhang
• William G. Glass
• Tristan Croll
• Aoife M. Harbison
• Elisa Fadda
• John D. Chodera

Simulation details

Structure preparation The RBD:hACE2 complex was constructed from individual RBD (PDB: 6m0j, Chain E) and hACE2 (PDB: 1r42, Chain A) monomers aligned to the full RBD:hACE2 structure (PDB: 6m0j). The 1r42 structure was used for hACE2 because:

• 1r42 is higher resolution (2.20 Å, whereas 6m0j is 2.45 Å).
• The electron density map of 1r42 clearly reveals the NAG orientation at each glycosylated asparagine residue, providing a reliable building block on which to construct more complex glycan structures.

While the NAGs for each glycan were present in 1r42, we constructed more complex glycans at each NAG due to the important role of glycans in mediating RBD:hACE2 binding, as shown by Casalino et al..

In order to start from the most reliable structural models, we obtained 6m0j and 1r42 from the Coronavirus Structural Taskforce (CST) database, which contains refined structural models based on careful examination of the electron density. In the RBD of the refined 6m0j structure amino acid rotamers and peptide bonds were flipped to increase Ramanchandran favourability, decrease rotamer outliers, and reduce clashes. A more detailed summary of the 6m0j refinement details is available at: https://github.com/thorn-lab/coronavirus_structural_task_force/blob/master/pdb/surface_glycoprotein/SARS-CoV-2/6m0j/isolde/notes.txt. The 1r42 refined structure differs from the PDB-deposited structure in that it includes the missing C-terminal domain of hACE2 (copied from the 6m17 PDB structure). A more detailed summary of the 1r42 refinement details is available here.

The resulting RBD and hACE2 monomers were then aligned in PyMOL 2.3.2 to the CST 6m0j structure to create an initial RBD:hACE2 complex. The overall RMSD was 0.426 Å and the interface RMSD was 0.405 Å, where RMSD was computed for all atoms and the interface residues were defined as all residues within 4 Å of the other binding partner.

Next, the full glycosylation patterns for hACE2 and RBD glycans were determined from Shajahan et al. and Watanabe et al. For the constructed RBD:hACE2 complex, these included sites: N53, N90, N103, N322, N432, N546, and N690 on hACE2 and N343 on the RBD. The glycan structures used for each site (FA2, FA26G1, FA2, FA2, FA2G2, A2, FA2, FA2G2, respectively) correspond to the most stable conformers obtained from multi microsecond MD simulations of cumulative sampling (Harbison et al.). Base NAG residues of each glycan structure were aligned to the corresponding NAG stub in the RBD:hACE2 model in PyMOL 2.3.2 and any resulting clashes were refined in ISOLDE.

System solvation and parametrization The refined glycosylated RBD:hACE2 complex was prepared for simulation using the AmberTools17 tleap suite. All relevant disulfide bridges were specified as well as covalent connectivity within each glycan structure. The glycosylated protein was parameterized with the Amber ff14SB (Maier et al., 2015) and GLYCAM_06j-1 (Kirschner et al., 2008) force fields. The system was solvated using the TIP3P rigid water model (Jorgensen et al., 1983) in a cubic box with 1.5 nm solvent padding on all sides. The solvated system was then minimally neutralized with 0.15 M NaCl using the Li/Merz ion parameters of monovalent ions for the TIP3P water model (12-6 normal usage set) (Li et al., 2015).

System equilibration The system was energy-minimized with an energy tolerance of 10 kJ mol−1and equilibrated using the OpenMM 7.4.2 Langevin integrator for 300 ns in the NPT (p=1 atm, T = 310 K) ensemble with a timestep of 4.0 femtoseconds, a collision rate of 1.0 picoseconds-1, and a constraint tolerance of 1 ✕ 10−5. Hydrogen atom masses were set to 4.0 amu by transferring mass from connected heavy atoms, bonds to hydrogen were constrained, and center of mass motion was not removed. Pressure was controlled by a molecular-scaling Monte Carlo barostat with an update interval of 25 steps. Non-bonded interactions were treated with the Particle Mesh Ewald method (Darden et al., 1993) using a real-space cutoff of 1.0 nm and the OpenMM default relative error tolerance of 0.0005, with grid spacing selected automatically. For improved stability, the structure was then equilibrated using the OpenMMTools 0.20.0 BAOAB Langevin integrator (Leimkuhler and Matthews, 2013) for 10 ns using all of the same simulation parameters described above. This simulation was subsequently packaged to seed for production simulation on Folding@home (Shirts and Pande, 2000, Zimmerman et al., 2020). Default parameters were used unless noted otherwise.

Folding@home simulation The equilibrated structure was then used to initiate parallel distributed MD simulations on Folding@home (Shirts and Pande, 2000, Zimmerman et al., 2020). Simulations were run with OpenMM 7.4.2 (Folding@home core22 0.0.13). Production simulations used the same Langevin integrator as the NpT equilibration described above. In total, 2000 independent MD simulations were generated on Folding@home. Conformational snapshots (frames) were stored at an interval of 0.5 ns/frame for subsequent analysis. The resulting final dataset contained 2000 trajectories, ~184 µs of aggregate simulation time, and 367610 frames. This amount of simulation time corresponds to ~13.7 GPU-years on an NVIDIA GeForce GTX 1080Ti. This trajectory dataset with solvent is available at the MolSSI COVID-19 Molecular Structure and Therapeutics Hub: https://covid.molssi.org//simulations/#foldinghome-simulations-of-the-sars-cov-2-spike-rbd-bound-to-human-ace2