All-atom molecular dynamics simulations and coarse-grained pulling simulations snapshots and contact maps analysis scripts for SARS-CoV-2 RBD variants in complex with H11-H4 nanobody

The first dataset contains 10,000 snapshots per trajectory, each representing a 1-microsecond-long simulation per replica, for all-atom molecular simulations of the RBD/H11-H4 complexes. The variants included are WT, Alpha, Delta, Omicron XBB.1.5, Omicron BA.2.86, and Omicron JN1.

1. SARS-CoV-2 RBD WT H11-H4_AA.tar.xz (one replica)
2. SARS-CoV-2 RBD Alpha H11-H4_AA.tar.xz (one replica)
3. SARS-CoV-2 RBD Delta H11-H4_AA.tar.xz (one replica)
4. SARS-CoV-2 RBD XBB.1.5 H11-H4_AA.tar.xz (one replica)
5. SARS-CoV-2 RBD BA.2.86 H11-H4_AA.tar.xz (five replicas)
6. SARS-CoV-2 RBD JN.1 H11-H4_AA.tar.xz (five replicas)

The simulations were conducted using AMBER22 and the FF19SB force fields (Salomon-Ferrer et al. 2013) with an enabled pmemd.cuda module for high performance. The initial structure of the RBD/H11-H4 complex was modeled from the ternary complex CR3022/H11-H4 and the WT RBD (PDB ID: 6ZH9 (Huo et al. 2020)). The coordinates of the CR3022 antibody were removed. To avoid artificial charges, the N- and C-termini residues of the RBD were capped with ACE and NME groups, respectively. Using the WT RBD as the structure template, we modeled several SARS-CoV-2 variants: i) Alpha, ii) Delta, iii) XBB.1.5, iv) BA.2.86, and v) JN.1 with UCSF Chimera v1.17.3 (Pettersen et al. 2004) using the Dunbrack rotamer library (Shapovalov and Dunbrack 2011). The corresponding protonation state corresponded to a pH value of 7.4 and was fixed with PDBfixer (Eastman et al. 2017). Each protein complex was solvated in a dodecahedral box initially extending 10 Å further from the solute in each direction, using the four-site OPC water model (Izadi, Anandakrishnan, and Onufriev 2014), and total charges in the system were neutralized with the appropriate number of counterions (9 Cl- ions for both the WT and Alpha variants, 9 Cl- ions for the Delta variant, 10 Cl- ions for the XBB.1.5 variant, 12 Cl- ions for the BA.2.86 variant, and 11 Cl- ions for the JN.1 variant). Afterward, the systems underwent geometric optimization using the steepest descent algorithm for 5,000 cycles to adjust the solvent orientation and remove local clashes. MD equilibration was carried out in multiple steps. The initial temperature equilibration of the complexes in the NVT ensemble involved a staged process in which the temperature was incrementally raised through four steps: 150, 200, 250, 300, and finally to 310 K. Each step was 200 ps long. In this process, position restraints were applied to the heavy atoms of each protein, using decreasing spring constants of 5.0, 4.0, 3.0, and 1.0 kcal/mol/Å2, respectively to allow for gradual relaxation of the complexes. Subsequently, a 1 ns equilibration at 310 K in the NPT ensemble without position restraints was conducted. Production MD trajectories were performed with the NPT ensemble using periodic boundary conditions and PME (Simmonett and Brooks 2021; Cerutti et al. 2009) with a grid spacing of 1.0 Å for treating long-range electrostatic interactions. The non-bonded interactions were described by the Lennard-Jones potential with a cutoff distance of 9 Å. The simulations employed Langevin dynamics (Sindhikara et al. 2009) for temperature control with a collision frequency of 4.0 ps-1, and the Berendsen barostat for pressure control (Berendsen et al. 1984), with a relaxation time of 2.0 ps and 1 bar pressure. Bond constraints involving hydrogen atoms were maintained using the SHAKE algorithm (Ryckaert, Ciccotti, and Berendsen 1977). The hydrogen mass repartition scheme was applied using ParmEd (Eastman et al. 2017), which allowed the use of a 4-fs integration time step (Hopkins et al. 2015). Each protein complex was simulated for 1 μs, except for the BA.2.86 and JN.1 variants, where five replicas of 1.5 μs long were run for each system. Most of the RBD/H11-H4 complexes remained stable during the AA-MD simulation, except for the Omicron variants BA.2.86 and JN.1, which dissociated upon equilibration. These variant have a lower sensitivity to antibody neutralization, enabling it to bypass both therapeutic and vaccine-induced immune responses, as reported in (Planas et al. 2024; Zhang et al. 2024). Hence, we excluded these last variants from the subsequent analysis.

The second dataset consists of 600 snapshots for each replica, totaling 50 coarse-grained pulling trajectories per system, across the four SARS-CoV-2 RBD/H11-H4 complexes (WT, Alpha, Delta, and Omicron XBB.1.5).

1. SARS-CoV-2 RBD WT H11-H4.tar.xz
2. SARS-CoV-2 RBD Alpha H11-H4.tar.xz
3. SARS-CoV-2 RBD Delta H11-H4.tar.xz
4. SARS-CoV-2 RBD XBB.1.5 H11-H4.tar.xz

The CG topology files for each protein complex were created with martinize2 (Kroon et al. 2022) and the Martini 3 force field (Souza et al. 2021). The secondary structure was identified using the DSSP v3.0.0 program (Touw et al. 2015). The GōMartini 3 approach (Poma, Cieplak, and Theodorakis 2017) was employed, substituting conventional harmonic bonds with Lennard-Jones (LJ) interactions based on contact maps obtained from AA-MD simulations. This method utilizes LJ potentials for virtual sites, which allows for the exploration of a broader conformational space, including critical unfolded states, enhancing our understanding of protein dynamics and functionality. A parameter study from standard values of 9.14 kJ/mol to 20 kJ/mol was conducted to determine the optimal value for the LJ potential depth (data not shown). The interaction energy of the contacts in GōMartini was set at 15.0 kJ/mol. This effective value allowed us to recover agreement with all-atom results of previously reported studies (Golcuk et al. 2022). The CG structures were initially minimized in vacuum during 5,000 steps using the steepest descent algorithm. Subsequently, the complexes were solvated in a 10x10x60 nm³ box using the Martini water model. This solvation involved about 38,511 coarse-grained water beads, corresponding to around 154,044 water molecules, with the addition of Na+ and Cl- ions to form a 0.15 M NaCl solution. Systems were then minimized using the same parameters as above. Positional restraints were imposed on the BB beads of each protein to prevent drifting during the equilibration phases. Both the NVT and NPT equilibrations, along with the production phase, utilized the V-rescale thermostat (Bussi, Donadio, and Parrinello 2007). The temperature coupling time constant was set at 1.0 ps for both protein and non-protein parts of the system, keeping the temperature at 300 K. The NVT equilibration was run for 2 ns, with an integration time of 20 fs. During the NPT equilibration and production, an isotropic pressure coupling was used with a compressibility set at 10-4 bar-1 and 1 bar pressure. The NPT equilibration was run for 5 ns, using the C-rescale barostat (Bernetti and Bussi 2020) with a pressure coupling time constant of 18 ps and an integration time of 10 fs. For the production phase, the Parrinello-Rahman barostat (Parrinello and Rahman 1981) was used, with a pressure coupling time constant of 15 ps. The cutoff distances for Coulomb and Van der Waals interactions were set at 1.2 nm across all equilibration and production phases. The pulling simulations (see Figure 1C) were conducted over 1.2 microseconds, with a time step of 20 fs. For the RBD/H11-H4 complexes, specific constraints were applied: the positions of the heavy atoms of the final three residues from the C-terminus of RBD were frozen along the z-axis, and similarly, the coordinates of heavy atoms of residues S126, S127, and K128 in H11-H4 were fixed along the x- and y-axes. The center of mass (COM) of these coordinates was targeted for steered molecular dynamics simulation at a constant speed of 1x10-4 nm/ps and a spring constant of 37.6 kJ/mol nm². A total of 50 independent replicas were conducted for each system using GROMACS 2023 (Abraham et al. 2015).

For further details on the trajectories, please contact Luis F. Cofas-Vargas (fcofas@ippt.pan.pl).

AMBER input files

1. 6ZH9_WT_AMBER_inputs.tar.xz
2. 6ZH9_alpha_AMBER_inputs.tar.xz
3. 6ZH9_delta_AMBER_inputs.tar.xz
4. 6ZH9_XBB1.5_AMBER_inputs.tar.xz
5. 6ZH9_BA.2.86_AMBER_inputs.tar.xz
6. 6ZH9_JN1_AMBER_inputs.tar.xz

*.pdb - Starting coordinates for MD simulations

*.parm7 - Topology file for AMBER

*.rst7 - Coordinates file for AMBER

COMMANDS.sh - This script runs a molecular dynamics simulation workflow consisting of minimization, equilibration, and production stages. It utilizes MPI and CUDA for parallel computing, while periodically checking and correcting the center of mass alignment during equilibration.

check_com.sh - 
This script aligns the center of mass of atoms in a current restart file to a reference structure, ensuring structural consistency by translating the coordinates if necessary.

min.in - Minimization file

eq*.in - NVT and NPT equilibration files

md.in - Production in NPT ensemble

Gromacs input files

1. WT input.tar.xz
2. Alpha input.tar.xz
3. Delta input.tar.xz
4. XBB.1.5 input.tar.xz

6ZH9_WT_pull.pdb - Starting structure for contact map calculation (http://pomalab.ippt.pan.pl/GoContactMap/)

6ZH9_WT_pull.map - Contact map

go_system.gro - Protein complex + Gō contacts molecular structure 

go_martini.itp - GōMartini 3 contacts topology

go_molecule1.itp - Topology file for RDB

go_molecule1.itp - Topology file for H11-H4

*.mdp -  Molecular dynamic parameters files required for running simulations

*.ndx - Index file for pulling simulations.

Gromacs output files

1. WT output.tar.xz
2. Alpha output.tar.xz
3. Delta output.tar.xz
4. XBB.1.5 output.tar.xz

*.gro - Final snapshots for minimization, NVT and NPT equilibrations

*.tpr - Topology and coordinate information for minimization, NVT and NPT equilibrations

pullf*.xvg - Pull-force files for each replica

pullx*.xvg - Pull-coordinate files for each replica

Contact map analysis scripts 

1. Native_contact_analysis.tar.xz

Native contacts

Distances_native_WT.ipynb - A script that calculates the distances and persistence of native contact pairs along their respective trajectories.

npt_dry.pdb - Topology for calculations

distances.py - Script that is designed to identify and calculate the critical breaking points for each contact pair in a molecular structure. This file must be in the "native_distances" folder

average_native_distance.py - Script for statistics calculation of breaking distance of contact pairs. This file must be in the  "distances" folder.

Graph.ipynb - Script to graph the contact persistence along trajectories. This file must be in the "distances" folder

Nonnative (NON) contacts

1. NON_contact_analysis.tar.xz

Distances_nonnative_WT.ipynb - A script that calculates the distances and persistence of nonnative contact pairs along their respective trajectories.

sort_contacts.py - Script that sorts native and nonnative contacts.

distances.py -  Script that is designed to identify and calculate the critical breaking points for each contact pair in a molecular structure. This file must be in the "nonnative" folder

count_nonnative_contacts.py -  Script to count the number of nonnative contacts along the trajectories. This file must be in the "nonnative" folder

average_distances.py - Script for statistics calculation of breaking distance of contact pairs. This file must be in the  "distances" folder.

Graph_nonnative.pynb - Script to graph the contact persistence along trajectories. This file must be in the "distances" folder

RMSD analysis

1. RMSD.tar.xz

RMSD.ipynb - Script to calculate and graph the RMSD of RBD and H11-H4