DOI: 10.5281/zenodo.15608522
title: All-atom molecular dynamics and coarse-grained steered MD simulations for "A Comparative Nanomechanical Study of Antibody and Nanobody Binding to SARS-CoV-2 Variants"

description: The first dataset contains 500 snapshots per system, each representing a 500-nanosecond-long simulation, for all-atom molecular simulations of the SARS-CoV-2 RBD/Abs or Nbs complexes. The variants included are Wuhan (WT), Omicron BA.4, and Omicron JN.1.



RBD-BA4_PDI-231_AA.tar.xz (one replica)

RBD-JN1_PDI-231_AA.tar.xz (one replica)

RBD-WT_PDI-231_AA.tar.xz (one replica)

RBD-BA4_R14_AA.tar.xz (one replica)

RBD-JN1_R14_AA.tar.xz (one replica)

RBD-WT_R14_AA.tar.xz (one replica)

RBD-BA4_S2X259_AA.tar.xz (one replica)

RBD-JN1_S2X259_AA.tar.xz (one replica)

RBD-WT_S2X259_AA.tar.xz (one replica)

RBD-BA4_C1_AA.tar.xz (one replica)

RBD-JN1_C1_AA.tar.xz (one replica)

RBD-WT_C1_AA.tar.xz (one replica)

RBD-BA4_R1-32_AA.tar.xz (one replica)

RBD-JN1_R1-32_AA.tar.xz (one replica)

RBD-WT_R1-32_AA.tar.xz (one replica)

RBD-BA4_n3113.1_AA.tar.xz (one replica)

RBD-JN1_n3113.1_AA.tar.xz (one replica)

RBD-WT_n3113.1_AA.tar.xz (one replica)


Each folder contains the following: 

*.parm7 - dry topology in amber format

*.rst7 - dry coordinates in amber format

*.nc - dry trajectories in netcdf format

*_amber.pdb - starting structure prepared with pdb4amber

frame_*.pdb - frame containing most of the high-frequency contacts

frame_*.map - contact map of the frame containing most of the high-frequency contacts 

 

We first generated complete structural models of antibody and nanobody complexes with the SARS-CoV-2 RBD, which served as an input for subsequent relaxation and nanomechanical analysis. A total of six representative complexes were selected, targeting three distinct epitopic regions of the RBD. These included three antibodies —PDI-231 (PDB ID: 7MZN)(Wheatley et al., 2021), S2X59 (PDB ID: 7M7W)(Starr et al., 2021), and R1-32 (PDB ID: 7YDI)(He et al., 2022), and three nanobodies, R14 (PDB ID: 7WD1)(H. Liu et al., 2023), C1 (PDB ID: 7OAP)(Huo et al., 2021), and n3113.1 (PDB ID: 7VNE)(Yang et al., 2021). Missing residues in 7M7W (residues 141-147) were modeled using Modeller v10.5(Sali & Blundell, 1993). In the 7YDI structure, ACE2 was removed, while for 7VNE, only the RBD-nanobody portion was retained. To avoid artificial termini charges, the N- and C-termini of the RBD, heavy chain, light chain, and nanobody (where needed)  were capped with acetyl and methylamide groups, respectively. Additionally, BA.4 and JN.1 SARS-CoV-2 variants were modeled using the WT RBD as template, incorporating the corresponding mutations via Modeller while preserving the original antibody or nanobody binding geometry. 

All-atom molecular dynamics (AA-MD) simulations were then performed to relax each complex and to obtain equilibrated structures suitable for coarse-grained modeling within the GōMartini 3 framework.  These simulations were carried out using the AMBER 24 software package with the ff19SB force field(Case et al., 2024; Tian et al., 2020). Simulations were run with the GPU-accelerated pmemd.cuda engine(Salomon-Ferrer et al., 2013). Initial energy minimization in vacuum was performed for each complex to eliminate steric clashes, particularly around mutated residues introduced during modeling. The protonation state of titrable residues was assigned at pH 7.4 using PDBFixer(Eastman et al., 2017).

Each minimized complex was solvated in a dodecahedral box of explicit water molecules using the four-site OPC water model(Izadi et al., 2014), with a minimum distance of 10 Å between any protein atom and the box edge. System neutrality was achieved by adding Cl- counterions, with the number of ions adjusted based on the net charge of each system. The number of ions added ranged from 2 to 12, depending on the complex variant.

Following solvation, systems underwent a 5,000 steps of energy minimization using the steepest descent algorithm to resolve unfavorable contacts and to allow for solvent relaxation. Temperature equilibration was performed under constant volume (NVT) conditions by incrementally increasing the temperature from 150 K to 310 K in five 200-ps state. During this phase, harmonic position restraints were applied to protein heavy atoms, with spring constants progressively reduced from 5 to 1 kcal mol-1 Å-2 to permit gradual relaxation of the protein backbone and side chains. This was followed by a 1 ns unrestrained equilibration at 310 K under constant pressure (NPT) conditions.

Production simulations were conducted in the NPT ensemble using periodic boundary conditions. Long rage electrostatics were treated with the particle mesh Ewald method, employing a grid spacing of 1 Å. Short-range interactions were modeled using a Lennard-Jones (LJ) potential with a 9 Å cutoff(Cerutti et al., 2009; Simmonett & Brooks, 2021). Temperature was regulated using Langevin dynamics with a collision frequency of 4 ps-1 (Sindhikara et al., 2009), while pressure was maintained at 1 bar using the Monte Carlo barostat with a relaxation time of 2 ps(Åqvist et al., 2004). Covalent bonds involving hydrogen atoms were constrained using the SHAKE algorithm(Ryckaert et al., 1977), and hydrogen mass repartitioning was applied using ParmEd(Shirts et al., 2017), allowing the use of a 4 fs integration time step(Hopkins et al., 2015). Each system was simulated for 500 ns. All complexes remained structurally stable throughout the simulations.

 

AA contact map determination

To identify high-frequency native contacts at the interface of each complex, we analyzed the atomistic molecular dynamics trajectories of antibody or nanobody–RBD systems. Each system was simulated for 500 nanoseconds, and contacts were calculated every 1 nanosecond using a custom ad hoc Python script (https://github.com/Multiscale-Modelling-of-Complex-Systems). This script uses the executable version of the Go Contact Map server (http://pomalab.ippt.pan.pl/GoContactMap/), available at https://zenodo.org/records/3817447, to apply both the van der Waals (vdW) radii overlap (OV) and repulsive chemical structural units (rCSU) contact map methods (Moreira et al., 2020; Wołek et al., 2015).

The OV contact map method is a purely geometric criterion that identifies interactions based on the spatial overlap of vdW spheres centered on heavy atoms. To incorporate attractive contributions, each vdW radius is scaled by a factor of 1.24. A native contact is defined when spheres from two residues, separated by at least four positions in the sequence, overlap. This approach has been widely used to identify stabilizing interactions relevant for folding and mechanical unfolding of protein domains (Cofas-Vargas, Olivos-Ramirez, et al., 2024; Z. Liu et al., 2022; Moreira et al., 2020; Poma et al., 2017, 2021).

The rCSU method complements the OV approach by integrating chemical specificity and electrostatics. It considers both attractive and repulsive interactions and classifies residues as hydrophobic, hydrophilic, aromatic, or ionic. A contact is considered valid when the number of attractive interactions exceeds the number of repulsive ones. Each atom is represented as a sphere, and its surface is sampled using a Fibonacci grid, which ensures a uniform and unbiased distribution of points.

Only high-frequency contacts, defined as those present in more than 70 percent of the frames, were retained. The frame containing the largest number of contacts was selected for coarse-grained modeling using the GōMartini framework

 

The second dataset comprises coarse-grained pulling simulations, with one snapshot recorded per nanosecond. For each system, a total of 50 trajectories were generated per pulling configuration—corresponding to each chain in the case of antibodies, across the three SARS-CoV-2 RBD complexes: WT, Omicron BA.4, and JN.1.



RBD-BA4_PDI-231_CG.tar.xz (50 replicas per pulling chain)

RBD-JN1_PDI-231_CG.tar.xz (50 replicas per pulling chain)

RBD-WT_PDI-231_CG.tar.xz (50 replicas per pulling chain)

RBD-BA4_R14_CG.tar.xz (50 replicas)

RBD-JN1_R14_CG.tar.xz (50 replicas)

RBD-WT_R14_CG.tar.xz (50 replicas)

RBD-BA4_S2X259_CG.tar.xz (50 replicas per pulling chain)

RBD-JN1_S2X259_CG.tar.xz (50 replicas per pulling chain)

RBD-WT_S2X259_CG.tar.xz (50 replicas per pulling chain)

RBD-BA4_C1_CG.tar.xz (50 replicas)

RBD-JN1_C1_CG.tar.xz (50 replicas)

RBD-WT_C1_CG.tar.xz (50 replicas)

RBD-BA4_R1-32_CG.tar.xz (50 replicas per pulling chain)

RBD-JN1_R1-32_CG.tar.xz (50 replicas per pulling chain)

RBD-WT_R1-32_CG.tar.xz (50 replicas per pulling chain)

RBD-BA4_n3113.1_CG.tar.xz (50 replicas)

RBD-JN1_n3113.1_CG.tar.xz (50 replicas)

RBD-WT_n3113.1_CG.tar.xz (50 replicas)


 

Each folder contains the following: 

*_CG.pdb - Starting structure 

*_CG.gro - Protein complex + Gō contacts molecular structure 

go_atomtypes.itp - GōMartini 3 contacts definitions

go_nbparams.itp - GōMartini 3 contacts topology

molecule_0.itp - Topology file for the complex

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

*.ndx - Index file for pulling simulations.

*.gro - Final snapshots for NVT and NPT equilibrations

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

npt_dry_ok.pdb - PDB with the correct numbering. Useful for visualization

Folder Pulling_H or Pulling_L (in the case of Abs) or Pulling (in the case of Nbs)

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

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

*.xtc- Pulling trajectories

pull_H.mdp - Molecular dynamics parameter file required for pulling from the H chain 

pull_L.mdp - Molecular dynamics parameter file required for pulling from the L chain

 

Coarse-grained (CG) models for all antibody and nanobody complexes were prepared using the Martini 3 force field(Souza et al., 2021) and the martinize2 tool(Kroon et al., 2024). Secondary structure definitions were assigned based on DSSP v3.0. For nanomechanical simulations, we applied the GōMartini 3 approach (Poma et al., 2017; Souza et al., 2025), which incorporates a LJ potential between virtual sites derived from native contact maps (Cofas-Vargas, Moreira, et al., 2024). The depth of the LJ potential wells (epsilon) for native contacts was set to 15 kJ mol-1, consistent with prior parameterizations of GōMartini 3 (Cofas-Vargas, Olivos-Ramirez, et al., 2024).

All CG systems were initially energy-minimized in vacuum for 5,000 steps using the steepest descent algorithm. To accommodate the pulling geometry, the complexes involving conventional Abs were solvated in a 10 x 10 x 100 nm3 water box, while the RBD-Nb complexes were solvated in a 10 x 10 x 80 nm3 box. In both cases, the standard Martini water model was used. Sodium and chloride ions were added to neutralize each system and to achieve a physiological strength of 0.15 M. A second minimization step was performed after solvation using the same settings. 

During the equilibration stages, position restraints were applied to the backbone (BB) beads of each protein to prevent diffusion or undesired rotations. Both NVT and NPT equilibrations, as well as the production phase, were performed using the V-rescale thermostat(Bussi et al., 2007) with a coupling time constant was set at 1 ps, maintaining the temperature at 300 K. NVT equilibration was run for 2 ns using a 20 fs time step. Subsequent NPT equilibration lasted 10 ns and was conducted with isotropic pressure coupling using the C-rescale barostat(Bernetti & Bussi, 2020), a compressibility of 1 x 10-4 bar, and a reference pressure of 1 bar. The pressure coupling time constant was set to 12 ps. Production simulations used the same barostat with a reduced coupling time of 8 ps and a 20 fs integration time step. The cutoff distance for both Coulombic and van der Waals (vdW) interactions was set to 1.2 nm throughout all stages of the CG simulations.

Steered molecular dynamics (SMD) pulling simulations were carried out for 2.5-3.0 μs per replica. For each complex, directional constraints were applied to mimic mechanical dissociation. Specifically, the positions of the heavy atoms of the last three residues at the C-terminus of the RBD were restrained along the z-axis. Simultaneously, the coordinates of the last three residues of the nanobody or of the antibody-s heavy or light chain were fixed in the x- and y-directions. The pulling force was applied to the center of mass (COM) of the selected terminus using a constant velocity of 1 x 10-5 nm ps-1 and a harmonic spring constant of 37.6 kJ mol-1 nm-2. For conventional Abs, two separate pulling sets were defined: in set 1, force was applied through the heavy chain; in set 2, force was applied through the light chain. A total of 50 independent replicas were performed for each configuration using GROMACS 2023.5(Abraham et al., 2015), with a pulling speed of 1x10-5 nm ps-1 and a force constant of 37.6 kJ mol-1 nm-2.

 

The third dataset (found in version 2 of the repository) consists of equilibrium coarse-grained simulations, with one snapshot saved every nanosecond (2 microseconds each). For each RBD/Ab or RBD/Nb complex, five independent trajectories were generated across the three SARS-CoV-2 RBD variants: WT, Omicron BA.4, and JN.1.



RBD-BA4_PDI-231_CG.tar.xz (5 replicas)

RBD-JN1_PDI-231_CG.tar.xz (5 replicas)

RBD-WT_PDI-231_CG.tar.xz (5 replicas)

RBD-BA4_R14_CG.tar.xz (5 replicas)

RBD-JN1_R14_CG.tar.xz (5 replicas)

RBD-WT_R14_CG.tar.xz (5 replicas)

RBD-BA4_S2X259_CG.tar.xz (5 replicas)

RBD-JN1_S2X259_CG.tar.xz (5 replicas)

RBD-WT_S2X259_CG.tar.xz (5 replicas)

RBD-BA4_C1_CG.tar.xz (5 replicas)

RBD-JN1_C1_CG.tar.xz (5 replicas)

RBD-WT_C1_CG.tar.xz (5 replicas)

RBD-BA4_R1-32_CG.tar.xz (5 replicas)

RBD-JN1_R1-32_CG.tar.xz (5 replicas)

RBD-WT_R1-32_CG.tar.xz (5 replicas)

RBD-BA4_n3113.1_CG.tar.xz (5 replicas)

RBD-JN1_n3113.1_CG.tar.xz (2 replicas)

RBD-WT_n3113.1_CG.tar.xz (5 replicas)


 

Each folder contains the following: 

*_CG.pdb - Starting structure 

*_CG.gro - Protein complex + Gō contacts molecular structure 

go_atomtypes.itp - GōMartini 3 contacts definitions

go_nbparams.itp - GōMartini 3 contacts topology

molecule_0.itp - Topology file for the complex

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

*.gro - Final snapshots for NPT equilibrations

*.tpr - Topology and coordinate information for N NPT equilibrations

npt_dry_ok.pdb - PDB with the correct numbering. Useful for visualization

*.xtc- Pulling trajectories

 

The fourth dataset (found in version 2 of the repository) consists of SMD coarse-grained simulations for H/L-only systems, with one snapshot saved every nanosecond (2 microseconds each). For each RBD/Ab complex, five independent trajectories were generated for the SARS-CoV-2 RBD WT.



RBD-WT_PDI-231_CG_onlyH.tar.xz (50 replicas)

RBD-WT_PDI-231_CG_onlyL.tar.xz (50 replicas)

RBD-WT_S2X259_CG_onlyH.tar.xz (50 replicas)

RBD-WT_S2X259_CG_onlyL.tar.xz (50 replicas)

RBD-BA4_R1-32_CG_onlyH.tar.xz (50 replicas)

RBD-BA4_R1-32_CG_onlyL.tar.xz (50 replicas)


Each folder contains the following: 

*_CG.pdb - Starting structure 

*_CG.gro - Protein complex + Gō contacts molecular structure 

go_atomtypes.itp - GōMartini 3 contacts definitions

go_nbparams.itp - GōMartini 3 contacts topology

molecule_0.itp - Topology file for the complex

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

*.ndx - Index file for pulling simulations.

*.gro - Final snapshots for NVT and NPT equilibrations

npt_dry_ok.pdb - PDB with the correct numbering. Useful for visualization

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

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

*.xtc- Pulling trajectories

pull.mdp - Molecular dynamics parameter file required for pulling from the H or L chain.

 

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