Interactive molecular dynamics in virtual reality for accurate flexible protein-ligand docking

Simulating drug binding and unbinding is a challenge, as the rugged energy landscapes that separate bound and unbound states require extensive sampling that consumes significant computational resources. Here, we describe the use of interactive molecular dynamics in virtual reality (iMD-VR) as an accurate low-cost strategy for flexible protein-ligand docking. We outline an experimental protocol which enables expert iMD-VR users to guide ligands into and out of the binding pockets of trypsin, neuraminidase, and HIV-1 protease, and recreate their respective crystallographic protein-ligand binding poses within 5–10 minutes. Following a brief training phase, our studies shown that iMD-VR novices were able to generate unbinding and rebinding pathways on similar timescales as iMD-VR experts, with the majority able to recover binding poses within 2.15 Å RMSD of the crystallographic binding pose. These results indicate that iMD-VR affords sufficient control for users to carry out the detailed atomic manipulations required to dock flexible ligands into dynamic enzyme active sites and recover crystallographic poses, offering an interesting new approach for simulating drug docking and generating binding hypotheses.

For performing iMD-VR runs, all proteins were parameterized with the Amber 14ffSB force field [2]. Benzamidine, indoleamidine, zanamivir, and amprenavir were parameterized with GAFF using antechamber. If needed, hydrogens were added to the systems using reduce [3]. Parameters for oseltamivir were taken from Woods et al [1]. All three systems used the OBC2 implicit water model. All systems were energy minimized before being used as the starting coordinates for each bound complex. All iMD-VR simulations were run at 300K with a timestep of 0.5 fs. When interacting with atoms, we utilized a gaussian force, where the amount of force applied depends on the distance between the user's controller and the atom they were selecting. Once users stopped interacting with a selection, the velocities of all atoms in that selection were reinitialized according to a Boltzmann distribution based on the temperature in order to prevent any remaining interaction energy in molecules propelling them past where the user intends to place them. Further details of how the velocities are reset can be found in ref [4].
For the expert-level tests, an 800 kJ/mol/nm 2 restraint was added to all backbone atoms in both trypsin and neuraminidase; additionally, this restraint was applied to any Ca2+ cations embedded within both systems. An 800 kJ/mol/nm 2 was applied to all backbone atoms in HIV-1 protease, excluding those which make up the flaps which gate the active site (defined as residues 49 to 55 in chain A and residues 48 to 54 in chain B). A separate force was applied to the HIV-1 flap backbone atoms which, unlike the other backbone restraint, the user could turn on and off during the interactive simulation. The software default value of 2000 kJ/mol/nm2 force constant was used for these restraints. A higher force constant is used compared to the backbone restraints due to their interactive nature. When a user wishes atoms to be held in place, it is expected that they want to be held in place firmly. All training and testing non-expert user tasks for trypsin and hiv-1 protease had the same backbone restraints as the expert tests; these restraints were removed for the neuraminidase blind task for one two novice cohorts that were recruited. The other cohort used the same backbone restraints as the expertlevel test. In future work, we intend to carry out a more detailed sensitivity analysis of different restraint regimes.

NarupaXR input files & settings
Input files for every system discussed in this work can be downloaded from the ESI. The input files consist of a PDB structure showing the starting bound coordinates, an .xml file which contains the simulation information to be loaded into NarupaXR, and a .json file which will automatically render the molecules into the same representations used in this work. Files are compatible with NarupaXR version 2019.0.0; details on how to load these files into NarupaXR can be found in the version 2019.0.0 documentation.

Equilibration protocols
For three of the protein systems, a single frame (the minimum RMSD along the iMD-VR trajectory) was extracted from the simulation trajectories for production run molecular dynamics, whose results are shown in Fig 4. The extracted frame from trypsin had an RMSD of 0.041 for benzamidine compared to the starting coordinates. For neuraminidase the extracted frame had an RMSD of 0.27 for oseltamivir and 0.302 for the 150-loop compared to the starting coordinates. For HIV-1 protease, the extracted frame had an RMSD of 0.218 for amprenavir and 1.383 for the active site flaps compared to the starting coordinates.
Extracted frames were solvated and neutralized: 10515 water molecules and 9 Cl-ions were added to trypsin and benzamidine, 17956 water molecules and 1 Cl-ion were added to neuraminidase and oseltamivir, 10696 water molecules and 5 Cl-ions were added to HIV-1 protease and amprenavir. All water was treated with the TIP3P water model. Forcefields for the protein and ligand were kept consistent between the interactive simulation and production run.
The solvated frames were minimized and equilibrated using the following protocol. Post minimization, the system was linearly heated over 20 ps, starting at 30K and increasing the temperature by 30K every 2 ps, until a temperature of 298K was reached. Next, a 10 kcal/mol/Å 2 restraint was added to the protein backbone atoms and the simulation was run for five sets of 10ps under the NVT ensemble, reinitializing the velocities of atoms in the system before each stage. Next, any restraints on the protein were removed and the box volume was allowed to relax to a pressure of 1 bar by running under the NPT ensemble for 9 nanoseconds. Finally, the system was switched back to the NVT ensemble, before being run for a final nanosecond. We note that there are a number of ways which we could have chosen in order to initiate the subsequent 200 ns MD simulations, and it is unclear a priori which choice is optimal. We chose the minimum simply because it is uniquely defined and unambiguous. Given that the snapshots which we extracted for production molecular dynamics were first solvated, minimized and then equilibrated for 10ns, the minimum RMSD structure would have been immediately relaxed to other nearby fluctuating structures; it is therefore unlikely that choosing the minimum as the starting point gives a markedly different result than choosing some other average structure nearby.

RMSD calculation
Narupa requires a set of atom coordinates from which to begin an iMD-VR simulation. These starting coordinates were used as the reference frame for the RMSD calculation. As such, the RMSD was calculated as a measure of the distance the ligand has moved compared to the starting frame (i.e. when it is in a bound position). To calculate the RMSD for the ligand, the protein backbone atoms in each step of the trajectory were aligned to those in the starting coordinates. Once the trajectory was aligned, the RMSD was taken for the ligand atoms compared to the starting coordinates (excluding hydrogens). However, an additional consideration is that all molecules tested display symmetry in some of their functional groups. Narupa allows users to quickly explore rotational orientations, which means that users can flip the orientation of certain rotationally invariant groups, and inadvertently introduce a degree of variation into RMSD calculations. Therefore, when calculating the RMSD of the ligand, symmetrical groups were identified (see figure  S1) and the swap plugin of VMD was used to calculate a more optimal RMSD value.
Additionally, to give an idea of the RMSD when docked, a 'baseline' RMSD was calculated by allowing each system to run for 50 ps in Narupa without any perturbation from the user, using a timestep of 0.5 fs and recording interval of 0.25 ps, totalling 200 RMSD data points. In order to gain an idea of the native behaviour of each ligand within a fully flexible protein, the backbone restraints were relaxed to 5 kJ/mol/nm 2 . In total, 60 ps of simulation time was recorded, however, the first 10 ps was discarded to give the simulation the opportunity to relax itself. Any ions in the system used backbone restraints of 800 kJ/mol/nm 2 . As with above, the RMSD was taken excluding the position of any hydrogens in the ligand. The RMSD value was either taken as (a) a single mean value over this 'reference' trajectory or (b) a range, calculated as plus or minus two standard deviations from the mean. These values can be found in Table S2. A dataset containing the RMSD values for each reference trajectory of each system can be found in the supplementary files of the ESI.

NarupaXR technical specifications
Narupa consists of two components: A molecular dynamics server engine that propagates a simulation forward and transmits the atomistic data to a client application, which renders the simulation in VR, and provides the user interface. For small to moderate sized protein simulations, these two components can be run locally on the same machine. Alternatively, a server machine can run the simulation engine and transmit the atomistic data to other VR clients, allowing multiple people to interact with the same iMD-VR simulation. For the expert-level user tests, both the molecular dynamics engine and VR front end were run on a high-end Alienware 15 R4 gaming laptop with an Intel i7-8750H 2.20GHz CPU and a NVIDIA GeForce GTX 1070 graphics card. For non-expert user tests, both the molecular dynamics engine and VR front end were run on a high-end desktop PC with an AMD Ryzen 31200 3.10GHz CPU and a NVIDIA GeForce GTX 1060 GPU. During the novice user tests, an expert explained the protein system to the user in a multiplayer simulation, which required two client machines, with the data from the simulation engine transmitted to a computer with identical specifications over a local area network.

7.
Participant-generated bound poses Figure 6 in the main text shows an overlay of the minimum RMSD bound pose achieved by each of the five novice participants, showing the increased variation when participants had ghost atoms (top three images) and when they did not (bottom three images). Additionally, for the final 5 ps of the participant-generated trajectories in cohort 2, the distance between atoms in the ligand molecule with key atoms in the protein residue are shown in Figure S4. Table S1 shows the minimum RMSD of each participant for each task.

Docking and redocking without backbone restraints
After establishing the usability of our framework, the experts repeated the tests with the 800 kJ/mol/nm 2 backbone restraints removed from all simulation files. Experts were able to manipulate all three systems with minimal interference to the protein tertiary structure, indicating that backbone restraints can potentially be removed. However, an added benefit of using relatively high backbone restraints is that paths generated are more straightforward to interpret because the motion from thermal fluctuations is curtailed. As a result, PathReducer provides principle components which better capture the key features of user-sampled iMD-VR pathways (i.e. the movement of a ligand from a protein binding site rather than side chain fluctuations). For cases where loop displacement is important (e.g., the motion of residues in the 150-loop during the novice neuraminidase testing tasks), this can be a potential contributor to poor task accomplishment. For example, comparing the re-establishment of the 150-loop for the neuraminidase testing task where both backbone restraints are on and off, figure S5 shows that the RMSD of these key residues is poorer and more varied when backbone restraints are not present. Furthermore, Figure S4 shows that contact between the zanamavir carboxylate group and Arg371 (one of the trio of arginines which is important for binding) is consistently poor, indicating non-optimal placement of the ligand by participants -which is further supported by figure S3, which shows placement of the carboxylate group as better where backbone restraints are present. This would suggest that where placement of key residues is important, such as the 150-loop, backbone restraints have the benefit of restricting the degrees of freedom in a protein. For novice users who are unfamiliar with both iMD-VR controls and the protein system itself, a high backbone restraint gives users the ability to familiarize themselves with the process of undocking and redocking without fear of accidentally pushing the protein system into a state where binding is unfavorable.

PathReducer analysis
PathReducer takes as input an xyz file containing a series of molecular structures (in this case, an iMD-VR-generated trajectory) and uses principle component analysis to determine the optimal rotation and scaling of the original coordinate system that captures the most structural variance in the fewest coordinates. This simple rotation and scaling means that the resultant coordinate system is comprised of coordinates that are linear combinations of the input coordinates (or input 'features'). By looking at the coefficient (or 'weight') of each input coordinate in the new coordinate system, one can see the relative amount that input coordinate is contributing to the direction of the new coordinate. The projections of the original data onto the resultant coordinates are referred to as principal components (PCs). Plots of iMD-VR pathways in PC space provides an idea of the similarity of structures sampled along the pathway. These plots can also be used to understand the 'route' taken by a trajectory. Figure 3 of the main text shows the iMD-VR-generated binding and unbinding pathways in the space defined by their top two and top three PCs. The structures input into PathReducer were represented as squared interatomic distance matrices, owing to their rotational and translational invariance compared to Cartesian coordinates.
For computational tractability, the PCA considered only those 70,000 interatomic distances that varied the most along a trajectory (only taking into account non-hydrogen atoms). In the plots in Figure 3 of the main text, purple corresponds to the beginning of the trajectory and yellow corresponds to the end, with blue and green passed through in between. Each scatter point is a frame from the iMD-VR trajectory, and a linear interpolation was used to connect each point.

Demographic Information
Novice users were recruited in two cohorts of five, totalling ten participants total. Cohort one completed the testing and training tasks for trypsin and neuraminidase, where all protein systems had a restraint placed on their backbone atoms. Cohort two completed the testing and training tasks for trypsin, neuraminidase and HIV-1 protease, however, for the neuraminidase testing task, all restraints were turned off. For cohort one, of the five participants, three participants were female and one was male. One participant opted not to specify their gender. Four participants were between the ages of 18 and 24 and one was between the age of 25 and 34; one participant was an undergraduate taught student, one participant was a postgraduate taught student, and three were postgraduate research students. Prior experience with VR was low amongst the cohort: three participants rated themselves as having a little experience using VR and one participant rated themselves as somewhat experienced. One member of the cohort rated themselves as very experienced. For cohort two, of the five participants, four were male and one was female; two participants were between the ages of 18 and 24 and three were between the ages of 25 and 34; two participants were postdoctoral researchers, two were postgraduate research students and one was a postgraduate taught student. Prior experience with VR was low amongst the cohort: two participants rated themselves as having no experience using virtual reality, two participants rated themselves as being slightly experienced, and one participant rated themselves as somewhat experienced. Figure S1 The symmetrical groups in each of the interactively undocked and redocked molecules. In these cases, the swap plugin of VMD was used so they could be treated as rotationally invariant when calculating RMSD. The atoms labelled with * are the key ligand atoms which participate in the distances shown in Fig S4.   Figure S4: atom distances between key ligand atoms (labelled with * in Fig S1) and key protein residues (shown in figure 1 and 2 of the main text) through the final 5 ps of trajectories generated by 5 user study participants from cohort two, where each color represents the results generated by a different participant. The top row shows the results for training tasks, and the bottom row shows the results for testing tasks. The solid black line shows the system-specific average (with black dotted lines representing ±2s) for a non-interactive Narupa simulation run for 50 ps in the bound pose and the calculated average atom distance is shown by the blue line. For benzamidine and indole-amidine, the distance represents the contact between their amidine group and Asp189 of trypsin. For oseltamivir and zanamavir, the distance represents the contact between their carboxylic acid group and Arg371. For amprenavir, the distance represents the contact between the hydroxyl oxygen and the protonated catalytic residue of HIV-1 protease, Asp25B.

Additional figures and tables
final 5 ps of simulation time pose which represents the minimum RMSD achieved by each of the five participants across the six tasks described in figure 2 of the main text. For each task, each of the five achieved poses are overlaid on top of one another. Figure S5 shows the RMSD of the two 150-loop residues, Asp151 and Arg152, of neuraminidase for the testing task, both with and without backbone restraints, showing that the RMSD of the 150-loop without restraints was (a) higher and (b) exhibited a greater degree of variation compared to with restraints. Each participant is assigned their own color. n/a n/a n/a C1 -2 0.38 3.29 1.22 1.12 n/a n/a n/a C1 -3 1.16 2.70 0.81 3.47 n/a n/a n/a C1 -4 0.50 1.02 0.52 2.47 n/a n/a n/a C1 -5 0.39 3.23 0.57 6.34 n/a n/a n/a C2 -  Table S1 The reported minimum ligand RMSD values for each participant for each task. Protein systems are categorized by letter, where A is trypsin, B is neuraminidase and C is HIV-1 protease. Participants are sorted by cohort, where C1 denotes the first cohort and C2 denotes the second cohort. Where a participant did not complete a given task, the RMSD value is n/a.

Table S2
The reported final ligand RMSD values for each participant for each task. Protein systems are categorized by letter, where A is trypsin, B is neuraminidase and C is HIV-1 protease. Participants are sorted by cohort, where C1 denotes the first cohort and C2 denotes the second cohort. Where a participant did not complete a given task, the RMSD value is n/a.