# Locating binding poses in protein-ligand systems using reconnaissance metadynamics

See allHide authors and affiliations

Contributed by Michele Parrinello, February 6, 2012 (sent for review December 15, 2011)

## Abstract

A molecular dynamics-based protocol is proposed for finding and scoring protein-ligand binding poses. This protocol uses the recently developed reconnaissance metadynamics method, which employs a self-learning algorithm to construct a bias that pushes the system away from the kinetic traps where it would otherwise remain. The exploration of phase space with this algorithm is shown to be roughly six to eight times faster than unbiased molecular dynamics and is only limited by the time taken to diffuse about the surface of the protein. We apply this method to the well-studied trypsin–benzamidine system and show that we are able to refind all the poses obtained from a reference EADock blind docking calculation. These poses can be scored based on the length of time the system remains trapped in the pose. Alternatively, one can perform dimensionality reduction on the output trajectory and obtain a map of phase space that can be used in more expensive free-energy calculations.

Understanding how proteins interact with other molecules (ligands) is crucial when examining enzymatic catalysis, protein signaling, and a variety of other biological processes. It is also the basis for rational drug design and is thus an important technological problem. Ligand binding is primarily examined using X-ray crystallography experiments together with measurements of the binding free energies. Additionally, numerous computational methods have been applied to this problem so as to extract more detailed information. The fastest of these approaches are based on an extensive configurational search of the protein surface (docking), in which the various candidate poses found are scored in accordance with some approximate function that treats solvation, protein flexibility, and entropic effects in some approximate manner.

Free-energy methods, based on either molecular dynamics (MD) or Monte Carlo simulations, can be used to calculate accurate binding free energies (1⇓–3). However, it is far more difficult to use these methods to search for candidate poses as the timescales involved in ligand binding are typically much longer than those that are accessible in MD. Thus, one often finds that the ligand becomes trapped in a kinetic basin on the surface of the protein and does not escape during the remainder of the calculation.

We recently developed a method, reconnaissance metadynamics, for increasing the rate at which high-dimensional configurational spaces are explored in MD simulations (4). This enhanced sampling is obtained by using a Gaussian mixture model to identify clusters in the stored trajectory. The positions of these clusters correspond to the kinetic basins in which the system would otherwise be trapped, which means that a history-dependent bias function that uses the information obtained from the clustering can be used to force the system away from the traps and into unexplored portions of phase space. In what follows we demonstrate how this algorithm can be used to examine the binding of benzamidine to trypsin, and perform a blind docking simulation based entirely on enhanced sampling.

## Background

Extensive conformational search procedures combined with fast and simple scoring functions give a surprisingly good description of protein-ligand docking in a variety of systems. In fact, for a number of systems so-called blind docking calculations can be performed in which the binding pose is found without using any experimental insight (5). The two greatest, unsolved problems for this field are to find universal scoring functions and to develop protocols for incorporating protein flexibility (6). These two problems are interlinked as an accurate scoring function must take the energetic cost of the conformational changes into account. Standard biomolecular force fields, together with implicit solvent models, provide the best approach for balancing these contributions. However, empirical and knowledge-based scoring functions often perform better for certain classes of problems and are thus frequently employed (7⇓–9).

Simulations based on MD force fields provide an alternative to simple docking calculations and both MD and Monte Carlo simulations have been used to locate sites with favorable interaction energy (10, 11). Furthermore, recent studies have exploited the power of modern computers to examine the process of ligand binding directly (12⇓–14). In these calculations the ligand is initially placed outside the protein and MD is used to find favorable binding sites. In the limit of long simulation time, the sites are visited according to the Boltzmann distribution and thus can be scored based on the amount of time the ligand spends at each site. This approach allows one to incorporate the protein flexibility, to treat the water explicitly, and to use established techniques for improving the force fields. In addition, one can obtain dynamical information on the binding process as well as structural information. However, these calculations still use an enormous amount of computational time and produce so much data that specialist tools are required for analysis. For example, the recent paper on the binding of benzamidine to trypsin by Buch et al. (12) used 500 unbiased simulations at a length of 100 ns.

Using plain MD simulations for locating binding poses is expensive because kinetic traps prevent the ligand from diffusing freely over the whole protein surface during short simulations. This problem is encountered frequently in MD and can be resolved by using enhanced sampling methods. A number of such methods have been applied to ligand binding (15⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–30). Typically these methods accelerate sampling by either increasing the temperature or by introducing a bias that prevents the system from becoming trapped in a basin. The bias is often constructed in terms of a small number of collective variables (CVs) that are selected by the user based on what is known about the location of the binding site, the binding pathway, and the conformational changes in the protein that occur during binding (31, 32). Using these methods one can calculate binding free energies for a small number of putative poses (33). Alternatively, one can find new, poorly characterized binding sites by using them in tandem with docking calculations (34).

The reconnaissance metadynamics method (RMD) (4) inserts the rich data that can be obtained from short MD simulations into a self-learning algorithm and thereby generates local collective coordinates that can push the system away from the kinetic traps it encounters. This procedure saves one from selecting a small number of appropriate CVs at the outset and thus provides a way to perform simulations when the reaction mechanism is uncertain. Thus far we have applied this method to model systems for polypeptide folding (4) and to small clusters of water and argon (35). These studies have demonstrated that RMD performs an extensive exploration of the energetically accessible portions of phase space and that this method can be used to locate global minima in energy landscapes. However, in the problems that we have examined the free-energy landscape is dominated by energetic contributions so these systems could alternatively be studied through a combination of optimization and transition state searches (36). Applying the RMD algorithm to the blind docking problem, as we do in this paper, represents a far greater challenge to the methodology because ligand binding involves a delicate balance between enthalpic and entropic contributions.

## Results and Discussion

We chose to examine the well-studied trypsin–benzamidine system, which has been extensively examined using free-energy perturbation (37). Both the benzamidine ligand and the trypsin protein are relatively rigid (38), so the binding site can be found using conventional blind docking (39).

One can use a large set of CVs in a reconnaissance metadynamics calculation and thus avoid many of the problems associated with choosing a small number of CVs for conventional metadynamics or umbrella sampling. However, it is important to realize that the CVs selected will influence the scope of the sampling. Thus for trypsin–benzamidine, where we know that the binding is not accompanied by large configurational changes in the protein, we selected CVs that describe only the position and orientation of the ligand relative to the protein and assume that MD alone will account for any protein flexibility. The CVs we chose are based on the distances between the C_{4}, N_{1}, and N_{2} atoms of the ligand (Fig. 1) and 16 uniformly spaced points on the protein surface (*Materials and Methods*). These distances are then transformed by a switching function so that whenever the ligand is far from the protein the collective variables have essentially the same values. The switching function is given by[1]where *r*_{i} is the *i* th distance and *r*_{0} = 13 *Å*. This set of 48 coordinates contains redundancy. However, because the self-learning algorithm at the heart of the RMD algorithm selects the most appropriate linear combinations of these to push, this redundancy does not present a particular problem. What is important to stress is that the parameters in this function and the points on the surface are chosen without using any information on the location of the binding site. As such this approach is general enough that it could be used for any globular protein. In addition, this description of the ligand’s position, orientation, and conformation can be systematically refined by either increasing the number of points on the surface of the protein or by increasing the number of points in the ligand. However, the cost of the calculations will increase as the number of CVs is increased (*Materials and Methods*).

### Extent of Exploration.

To test whether or not RMD is doing a good job of exploring phase space we generated a set of putative binding poses using conventional blind docking. These calculations were done using EADock (40), which is known to reproduce the correct binding pose for a range of systems (41) and which generated a large number of structurally diverse poses (Table S1). We then ran 10, 200-ns reconnaissance metadynamics simulations and calculated the rmsd distance between snapshots taken every 10 ps from our trajectory and the 27 poses found in our EADock calculations. Fig. 2 shows that during our simulations we come close to every single one of these putative poses. More importantly, in five out of the ten simulations we were able to find the binding site. These results are in stark contrast to the results we obtain from MD simulations of similar length. During the course of these calculations we were only able to find a subset of the poses and the binding site was never visited. This observation appears, at first glance, to be at variance with the results of Buch et al. (12) who found that in 37% of their 100-ns, unbiased MD simulations on this system the experimental binding pose was found. However, in their simulations some information on the location of the binding site was employed, as constraints were applied on the relative position of the protein and ligand to ensure that the ligand only explored one side of the protein.

Fig. 3*A* provides an alternative representation of the data on the extent to which phase space is explored during the RMD and MD simulations. This figure shows the fraction of reference poses found as a function of time and suggests that RMD is on average six to eight times faster at finding poses than MD [fitting the curves in Fig. 3*A* to the function 1 - exp(-*t*/*t*_{0}) we find that the ratio of *t*_{0} values for MD to RMD are 5.6, 7.5, and 8.3 for the three rmsd cutoffs we tested]. This increased speed does not appear particularly dramatic, but it is important to remember that if in any of the MD simulations the ligand had found the binding site it would have almost certainly remained there for the remainder of the simulation time, which is not what happens in the RMD simulations that find the binding site. In other words, the exploration in RMD is only slowed down because diffusion of the ligand about the protein is relatively slow—a fact of life that will be present in any method based on molecular dynamics.

### Generating Candidate Poses.

To generate meaningful output from any ligand-binding trajectory, it is necessary to predict which poses have high binding affinities, much as one scores poses in traditional docking calculations. Making such predictions from MD simulations is in principle straightforward, because the time spent in a given configuration is connected to its free energy. The only caveat is that one must see multiple transitions between states. If one appropriately accounts for the bias, similar strategies can be used in methods involving a bias potential. The problem with RMD is that multiple transitions between states are seldom observed because of the high-dimensionality space of collective variables. In contrast, when using methods like metadynamics, the small number of collective coordinates forces these transitions to occur.

In an RMD simulation it will take some time to generate sufficient bias to push the system out of a basin. The specific amount of time will depend on the basin’s depth and hence its kinetic stability. Low free-energy poses are usually narrow minima in the potential energy surface. These states will be both thermodynamically and kinetically stable. It may, therefore, be possible to find low free-energy poses by extracting the most populated clusters from an RMD trajectory. To further explore this idea, we analyzed the RMD trajectory frames using the method of Daura et al. (42) that is implemented in the GROMACS `g_cluster` utility. This procedure ranks each trajectory frame based on the number of neighboring frames that are within 1 Å rmsd. The top-ranked frame, together with all its neighbors, is then removed and the ranking process is repeated.

Fig. 3*B* shows that the clusters generated from the analysis of the RMD trajectories are much smaller than those generated from an analysis of the MD trajectories. This result confirms that the MD simulations are spending a great deal of time (up to 30 ns) trapped at a small number of sites on the protein surface. In contrast RMD spends at most 0.4 ns in any given pose and is thus able to explore more of the protein surface. In addition, this analysis of the RMD simulations identifies the binding pose as important. In three of the five RMD simulations that found the binding site, the cluster corresponding to the binding site is the most populated, whereas in the remaining two the binding site is ranked second and third.

Fig. 3*C* provides further evidence that clustering of the RMD trajectory gives reasonable binding poses. In this figure we show the vacuum interaction energy between the protein and the ligand for the top 50 clusters (i.e., the most populated ones) from each simulation. This interaction energy neglects solvent and entropic effects but is still often correlated with the binding free energy (43). Hence, the fact that the clusters found in RMD have consistently lower energies than those found in MD suggests that they correspond to more strongly bound conformations. Furthermore, if we examine all the frames in the trajectory we find that, in contrast with MD, the top clusters in RMD correspond to the structures with the lowest energies. There is no such shift in MD, which suggests that in these simulations the ligand becomes trapped in many basins that do not have particularly low interaction energies. As such, the MD simulations are too short to express the relationship between the residence time in a given structure and its free energy.

The clustering procedure does not take into account the bias, and thus some of the well-populated clusters might not correspond to minima on the unbiased free-energy surface. Hence, to probe the kinetic stabilities of the poses from one of the RMD simulations, we ran unbiased MD trajectories starting from the 136 most populated clusters. During these simulations we took the time spent within 2.5 Å rmsd of the initial configuration as a measure of the stability of the pose and found that 89 poses were stable for more than 100 ps, 25 were stable for more than 1 ns, and 7 of them were stable for more than 5 ns. Out of these seven poses, one was the crystallographic pose and one was a similar pose in which the ligand was separated from the Asp-189 residue by a water molecule. In addition, this set of poses contained the S2 and S3 states that were identified as stable in the MD studies of Buch et al. (12) (Table S2). Intriguingly, a stable pose (Fig. S1) was found in a part of the protein surface that was deliberately not explored in the MD investigation in reference (12). This configuration remains unchanged for 60 ns of unbiased MD and we predict that it is one of most stable interaction sites outside of the binding pocket. It is possible that, like the S2 site, it acts as a secondary binding site (44). For the EADock calculations a similar analysis showed that only eight of the poses generated were stable for more than 100 ps (Table S1) and that this set included the binding site, the S2 site, and a similar pose to that shown in Fig. S1.

### Dimensionality Reduction.

Clustering is one way of examining the data from an extensive sampling of a high-dimensional phase space, such as that obtained from docking, MD, or an enhanced sampling calculation. An alternative is to perform dimensionality reduction (45, 46). This way of examining ligand binding is appealing, because the largest changes in the position of the ligand are those corresponding to motion across the two-dimensional protein surface so the data should lie on a low-dimensionality manifold. Furthermore, a low-dimensional representation of the protein surface is a useful tool for visualizing the kinetic information that can be extracted from MD-based approaches.

Many dimensionality reduction algorithms work by endeavoring to reproduce the rmsd distances between the trajectory frames in a lower-dimensionality space (47). Clearly the rmsd distances between the ligand in the various trajectory frames can be approximately reproduced in a three-dimensional space as they will be dominated by differences in position of the center of mass of the ligand. To further lower the dimensionality of the projection requires one either to incorporate periodicity in the low-dimensionality projection or to make less of the global features in phase space and, instead, focus on the local connectivity between basins. We recently developed the sketch-map algorithm (48) as a tool for analyzing trajectory data. This algorithm uses the second of these approaches to the problem as it endeavors to reproduce the immediate connectivity between states rather than the full set of distances between frames. The algorithm’s focus is controlled by transforming the distances in the high-dimensionality and low-dimensionality spaces using a sigmoid function. This procedure also ensures that close-together points are projected close together, whereas far apart points are projected far apart but not necessarily at the same distance.

We used the RMD trajectories to produce the sketch-map projections because, unlike our MD simulations which didn’t visit the binding site, we have sufficient sampling in the RMD to build a reliable map. To record the high-dimensional positions, we used the coordinates of the ligand’s C_{4}, N_{1}, and N_{2} atoms in a protein-centered frame of reference. Fig. 4*A* shows that the resulting two-dimensional map clearly separates the poses around the binding site from other low energy poses on the protein surface and that there are specific pathways and channels that connect the various clusters. Moreover, Fig. 4 *B* and *C* and Fig. S2 show that, in the area around the binding site, we are able to separate the metastable sites described by Buch et al. (12), in spite of the fact that some of them are rather close in space (center of mass separation of ∼4 *Å*). This result suggests that sketch-map is also able to describe the orientation of the ligand and that using multiple atoms to define the ligand’s position is worthwhile. The resolution can be further improved by constructing a map using only points that are close to the binding site (Fig. S3).

We can use the projection shown in Fig. 4 to do a qualitative comparison between the results of our RMD simulation and the results of the extensive MD simulations by Buch et al. (12). In agreement with the previous study there is a significant population in the S3 state and a pathway from this state to the binding pose that passes through the TS1, TS2, and TS3 transition states. There are also other pathways between the bulk solvent and the binding site that pass through TS2 and TS3. In particular, during six of the ten binding or unbinding events that we observed, the ligand passed through the TS2 state on its way to or from the binding site, which suggests that this state is on the main binding pathway.

## Conclusions

Molecular dynamics with explicit solvent has enormous potential for predicting protein-ligand interactions because it is based on a physically motivated and systematically improvable potential energy surface and because it incorporates conformational, solvent, and entropic effects in a physically consistent manner. Its one major drawback is that it is considerably more computationally expensive than using docking calculations based on a configurational search with approximate scoring functions. One reason for this expense is that there are many energetic basins on the surface of the protein which can kinetically trap the ligand and slow down diffusion. This problem can be resolved by using a simulation bias to force the system away from kinetic traps and to flatten the energy surface. However, the requirement to find a small set of CVs that describes all the potential traps makes it difficult to apply a suitable bias using many established methods. In contrast, in reconnaissance metadynamics we can use large numbers of collective variables and let the algorithm work out which linear combination best describes each trap. The procedure outlined in this paper can thus be used to tackle problems where conformational and solvent effects play a large role, which would be difficult to examine using standard docking. Furthermore, the method is considerably cheaper than unbiased MD.

Reconnaissance metadynamics simulations provide an extensive exploration of the low-energy portions of phase space. One can use this data to find the approximate locations for the various basins in the free-energy surface or alternatively use dimensionality reduction techniques to create low-dimensionality maps of phase space. The fact that these maps are low-dimensional allows one to reexplore the interesting parts of phase space using other, more quantitative, enhanced sampling algorithms. In future, we will use this idea to extract accurate free energies for the various binding poses found during the RMD simulations.

## Materials and Methods

### System Setup and Computational Details.

The simulations were performed using GROMACS 4.5 (49) and the PLUMED plug-in (50). We used the Amber ff99 force field (51) for the protein and TIP3P for the water molecules. For the ligand, van der Waals parameters were taken from the corresponding amino acids (phenylalanine and arginine), and appropriate charges were calculated using a RESP fit (52) to a Hartree–Fock calculation with the 6-31G* basis set—a procedure identical to that described in ref. 27. Long-range electrostatics was treated using the particle mesh Ewald approach with a grid spacing of 1.2 Å. A cutoff of 10 Å was used for all van der Waals and the direct electrostatic interactions and the neighbor list was updated every 10 steps. All production simulations were performed in the canonical ensemble at 300 K and this temperature was maintained using the stochastic velocity rescaling thermostat (53). To prevent the system from sampling fully solvated configurations we used a restraining wall that limited the exploration to configurations where the sum of all the switching functions between the C_{7} carbon and the points on the surface was greater than 1. This wall only has any effect when the minimum distance between the protein and the ligand is greater than 12 Å and represents a relatively small perturbation of the underlying energy surface.

The trypsin–benzamidine complex [Protein Data Bank (PDB) ID code 1J8A] (54) was used as the starting structure in this study. All histidines were protonated on the *N*_{ϵ} site other than the catalytic H57, which was doubly protonated. This protein was then placed in a truncated octahedral simulation box that extended at least 7 Å from any protein atom. Prior to production a 10 ns constant pressure simulation, in which the protein atoms were initially restrained, was performed to equilibrate the system. Ten RMD production simulations were performed together with 10 MD simulations. These calculations were started from ten statistically inequivalent configurations, where the ligand was outside the protein. For each calculation we ran one RMD and one MD simulation. The initial starting configuration was generated by displacing the ligand from the binding site by 20 Å and running a short equilibration run. The remaining nine starting points were selected from the MD trajectory launched from the first point. In all these initial configurations the protein-ligand distance was greater than 10 Å. Furthermore, we visually inspected the starting configurations to ensure the widest possible spread of initial configurations.

### RMD Setup.

Relevant points on the surface of the protein were selected by constructing a graph which had all the C_{α} atoms at its vertices and connections between any pair of vertices closer than 14 Å. A heuristic algorithm was then used to find the maximum independent set of this graph (55). This procedure produces a uniformly distributed set of C_{α} atoms on the surface. For trypsin these were the C_{α} atoms of residues 23, 47, 60, 74, 92, 97, 109, 127, 147, 159, 164, 173, 186, 193, 229, and 244. The switching function was set up so that its value for a test point moving along the protein surface (5 Å above it) changed smoothly from approximately 1 when it was immediately above one of the surface points to approximately 0.4 once it was above the neighboring surface point 14 Å away. For the reconnaissance metadynamics, data was collected every 0.5 ps, which was then clustered every 100 ps. The bias was constructed from the clusters that had a weight greater than 0.2 in these fits and by endeavoring to add hills of width 1.5 and height 1 kJ mol^{-1} every 2 ps. Hills were only added when the distance from one of the cluster centers (in the metric of that particular cluster) was less than 8.356—a distance that, at variance with previous applications of RMD, was kept constant for the entirety of the simulation.

As discussed in the main text we can easily create a more fine grained representation of the space by increasing the number of CVs and thus increasing the cost of the calculation. It is not straightforward to quantify the scaling with the number of CVs because it is unclear how much longer it will take to sample these higher dimensionality spaces. What we can say with certainty is that calculating the distance between a basin center and the instantaneous position scales with the square of the number of CVs. However, the cost of calculating the force because of the bias is for the most part small when compared to the cost of a single MD step.

### Docking Calculations.

The docking calculations presented in this paper were used to provide a set of interesting poses that we could refind using our RMD simulations. We thus chose not to dwell on these calculations and just used the default (fast) protocol for EADock, which is provided on the Swissdock web server (56). The crystallographic structure of the protein (with the ligand removed) was used directly and 256 binding poses were obtained. These poses were then clustered using an rmsd cutoff of 2 Å and only clusters with at least five members were used. More details on these structures can be found in Table S1, which also shows that the crystallographic pose has an energy that is considerably lower than that of the other poses.

### Sketch-Map Calculations.

The distances, *d*, between frames in the nine-dimensional space were transformed using 1 - [1 + (2^{a/b} - 1)(*d*/*σ*)^{a}]^{-b/a} with *σ*, *a*, and *b* taking values of 20 Å, 1, and 3, respectively. The projection was then generated by minimizing the discrepancies between these transformed distances and the set of distances between the frames’ projections. These distances in the low-dimensionality space were once again transformed by the sigmoid function above, but in this case the *a* and *b* parameters were set to 2 and 3, respectively. The data from the 10 RMD trajectories was fitted by first projecting a set of 500 landmark points, 100 of which were selected at random and 400 of which were selected using farthest point sampling. Each point in this fit was weighted based on the number of unselected frames that fell within its voronoi polyhedra. Once this fitting was completed the unselected trajectory frames were mapped using the out-of-sample projection technique detailed in ref. 48.

## Acknowledgments

We thank Dr. M. Ceriotti for help with the sketch-map calculations and Dr. V. Limongelli for fruitful discussions. We also acknowledge computational resources from the central high-performance cluster of Eidgenössische Technische Hochschule Zurich (Brutus). Financial support for this work was obtained from European Union Grant ERC-2009-AdG-247075 and from The Swedish Research Council Grant 623-2009-821.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. E-mail: par.soderhjelm{at}phys.chem.ethz.ch or parrinello{at}phys.chem.ethz.ch.

Author contributions: P.S., G.A.T., and M.P. designed research; P.S., G.A.T., and M.P. performed research; P.S., G.A.T., and M.P. analyzed data; and P.S., G.A.T., and M.P. wrote the paper.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1201940109/-/DCSupplemental.

## References

- ↵
- ↵
- Singh N,
- Warshel A

- ↵
- Essex JW,
- Severance DL,
- Tirado-Rives J,
- Jorgensen WL

- ↵
- Tribello G,
- Ceriotti M,
- Parrinello M

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Buch I,
- Giorgino T,
- De Fabritiis G

- ↵
- Dror R,
- et al.

- ↵
- ↵
- ↵
- Woods C,
- Jonathan W,
- King M

- ↵
- ↵
- ↵
- ↵
- Torrie GM,
- Valleau JP

- ↵
- ↵
- ↵
- Laio A,
- Parrinello M

- ↵
- Woo HJ,
- Roux B

- ↵
- Jensen M,
- Park S,
- Tajkhorshid E,
- Schulten K

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Limongelli V,
- et al.

- ↵
- ↵
- ↵
- ↵
- Wales DJ

- ↵
- ↵
- Guvench O,
- Price D,
- Brooks C III.

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Oliveira M,
- et al.

- ↵
- Das P,
- Moll M,
- Stamati H,
- Kavraki L,
- Clementi C

- ↵
- Ferguson A,
- Panagiotopoulos A,
- Debenedetti P,
- Kevrekidis I

- ↵
- Cox T,
- Cox M

- ↵
- Ceriotti M,
- Tribello G,
- Parrinello M

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Cuesta-Seijo JA,
- García-Granda S

- ↵
- Balaji S,
- Swaminathan V,
- Kannan K

- ↵
- Grosdidier A,
- Zoete V,
- Michielin O

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Chemistry

- Biological Sciences
- Biophysics and Computational Biology